]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * written by International Business Machines Corp. |
f4cf5f2d | 4 | * Copyright (C) 2001-2012 Free Software Foundation, Inc. |
e4d82761 UD |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify | |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
50944bca | 10 | * |
e4d82761 UD |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
e4d82761 UD |
18 | */ |
19 | /********************************************************************/ | |
50944bca | 20 | /* */ |
e4d82761 UD |
21 | /* MODULE_NAME: dosincos.c */ |
22 | /* */ | |
23 | /* */ | |
24 | /* FUNCTIONS: dubsin */ | |
25 | /* dubcos */ | |
26 | /* docos */ | |
27 | /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */ | |
28 | /* sincos.tbl */ | |
29 | /* */ | |
30 | /* Routines compute sin() and cos() as Double-Length numbers */ | |
31 | /********************************************************************/ | |
32 | ||
33 | ||
34 | ||
35 | #include "endian.h" | |
50944bca | 36 | #include "mydefs.h" |
c8b3296b | 37 | #include <dla.h> |
e4d82761 | 38 | #include "dosincos.h" |
1ed0291c | 39 | #include <math_private.h> |
15b3c029 | 40 | |
31d3cc00 UD |
41 | #ifndef SECTION |
42 | # define SECTION | |
43 | #endif | |
44 | ||
af968f62 UD |
45 | extern const union |
46 | { | |
47 | int4 i[880]; | |
48 | double x[440]; | |
49 | } __sincostab attribute_hidden; | |
50 | ||
e4d82761 UD |
51 | /***********************************************************************/ |
52 | /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ | |
53 | /* as Double-Length number and store it at array v .It computes it by */ | |
54 | /* arithmetic action on Double-Length numbers */ | |
50944bca | 55 | /*(x+dx) between 0 and PI/4 */ |
e4d82761 UD |
56 | /***********************************************************************/ |
57 | ||
31d3cc00 UD |
58 | void |
59 | SECTION | |
60 | __dubsin(double x, double dx, double v[]) { | |
a1a87169 | 61 | double r,s,c,cc,d,dd,d2,dd2,e,ee, |
e4d82761 | 62 | sn,ssn,cs,ccs,ds,dss,dc,dcc; |
58985aa9 | 63 | #ifndef DLA_FMS |
a1a87169 UD |
64 | double p,hx,tx,hy,ty,q; |
65 | #endif | |
50944bca UD |
66 | #if 0 |
67 | double xx,y,yy,z,zz; | |
68 | #endif | |
e4d82761 UD |
69 | mynumber u; |
70 | int4 k; | |
50944bca | 71 | |
e4d82761 UD |
72 | u.x=x+big.x; |
73 | k = u.i[LOW_HALF]<<2; | |
74 | x=x-(u.x-big.x); | |
75 | d=x+dx; | |
76 | dd=(x-d)+dx; | |
a1a87169 | 77 | /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ |
e4d82761 | 78 | MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); |
af968f62 UD |
79 | sn=__sincostab.x[k]; /* */ |
80 | ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */ | |
81 | cs=__sincostab.x[k+2]; /* */ | |
82 | ccs=__sincostab.x[k+3]; /* */ | |
e4d82761 | 83 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */ |
50944bca | 84 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); |
e4d82761 | 85 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */ |
50944bca | 86 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 UD |
87 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */ |
88 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
89 | ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */ | |
50944bca | 90 | |
e4d82761 | 91 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */ |
50944bca | 92 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); |
e4d82761 | 93 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */ |
50944bca | 94 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); |
e4d82761 | 95 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */ |
50944bca | 96 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); |
e4d82761 | 97 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */ |
50944bca | 98 | |
e4d82761 UD |
99 | MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); |
100 | MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
50944bca | 101 | SUB2(e,ee,dc,dcc,e,ee,r,s); |
e4d82761 | 102 | ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */ |
50944bca | 103 | |
e4d82761 UD |
104 | v[0]=e; |
105 | v[1]=ee; | |
106 | } | |
107 | /**********************************************************************/ | |
108 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ | |
109 | /* as Double-Length number and store it in array v .It computes it by */ | |
110 | /* arithmetic action on Double-Length numbers */ | |
111 | /*(x+dx) between 0 and PI/4 */ | |
112 | /**********************************************************************/ | |
113 | ||
31d3cc00 UD |
114 | void |
115 | SECTION | |
116 | __dubcos(double x, double dx, double v[]) { | |
a1a87169 | 117 | double r,s,c,cc,d,dd,d2,dd2,e,ee, |
e4d82761 | 118 | sn,ssn,cs,ccs,ds,dss,dc,dcc; |
58985aa9 | 119 | #ifndef DLA_FMS |
a1a87169 UD |
120 | double p,hx,tx,hy,ty,q; |
121 | #endif | |
50944bca UD |
122 | #if 0 |
123 | double xx,y,yy,z,zz; | |
124 | #endif | |
e4d82761 | 125 | mynumber u; |
50944bca UD |
126 | int4 k; |
127 | u.x=x+big.x; | |
e4d82761 UD |
128 | k = u.i[LOW_HALF]<<2; |
129 | x=x-(u.x-big.x); | |
50944bca | 130 | d=x+dx; |
e4d82761 | 131 | dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ |
50944bca | 132 | MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); |
af968f62 UD |
133 | sn=__sincostab.x[k]; /* */ |
134 | ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */ | |
135 | cs=__sincostab.x[k+2]; /* */ | |
136 | ccs=__sincostab.x[k+3]; /* */ | |
e4d82761 | 137 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 138 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); |
e4d82761 | 139 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 140 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 UD |
141 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
142 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
50944bca UD |
143 | ADD2(ds,dss,d,dd,ds,dss,r,s); |
144 | ||
e4d82761 | 145 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 146 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); |
e4d82761 | 147 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 148 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); |
e4d82761 | 149 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 150 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); |
e4d82761 | 151 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 152 | |
e4d82761 UD |
153 | MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); |
154 | MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
155 | ||
156 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
157 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); | |
158 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
50944bca | 159 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 | 160 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
50944bca UD |
161 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
162 | ADD2(ds,dss,d,dd,ds,dss,r,s); | |
e4d82761 UD |
163 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
164 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); | |
165 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
166 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); | |
167 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
168 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); | |
169 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
170 | MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); | |
171 | MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
50944bca UD |
172 | ADD2(e,ee,dc,dcc,e,ee,r,s); |
173 | SUB2(cs,ccs,e,ee,e,ee,r,s); | |
174 | ||
e4d82761 UD |
175 | v[0]=e; |
176 | v[1]=ee; | |
177 | } | |
178 | /**********************************************************************/ | |
179 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ | |
180 | /* as Double-Length number and store it in array v */ | |
181 | /**********************************************************************/ | |
31d3cc00 UD |
182 | void |
183 | SECTION | |
184 | __docos(double x, double dx, double v[]) { | |
e4d82761 UD |
185 | double y,yy,p,w[2]; |
186 | if (x>0) {y=x; yy=dx;} | |
187 | else {y=-x; yy=-dx;} | |
188 | if (y<0.5*hp0.x) /* y< PI/4 */ | |
a1a87169 | 189 | {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];} |
e4d82761 UD |
190 | else if (y<1.5*hp0.x) { /* y< 3/4 * PI */ |
191 | p=hp0.x-y; /* p = PI/2 - y */ | |
192 | yy=hp1.x-yy; | |
193 | y=p+yy; | |
194 | yy=(p-y)+yy; | |
ca58f1db | 195 | if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} |
a1a87169 UD |
196 | /* cos(x) = sin ( 90 - x ) */ |
197 | else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; | |
e4d82761 UD |
198 | } |
199 | } | |
200 | else { /* y>= 3/4 * PI */ | |
201 | p=2.0*hp0.x-y; /* p = PI- y */ | |
50944bca | 202 | yy=2.0*hp1.x-yy; |
e4d82761 UD |
203 | y=p+yy; |
204 | yy=(p-y)+yy; | |
ca58f1db | 205 | __dubcos(y,yy,w); |
e4d82761 UD |
206 | v[0]=-w[0]; |
207 | v[1]=-w[1]; | |
208 | } | |
209 | } |