]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 AJ |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001 Free Software Foundation | |
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 | |
17 | * along with this program; if not, write to the Free Software | |
50944bca | 18 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
e4d82761 UD |
19 | */ |
20 | /********************************************************************/ | |
50944bca | 21 | /* */ |
e4d82761 UD |
22 | /* MODULE_NAME: dosincos.c */ |
23 | /* */ | |
24 | /* */ | |
25 | /* FUNCTIONS: dubsin */ | |
26 | /* dubcos */ | |
27 | /* docos */ | |
28 | /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */ | |
29 | /* sincos.tbl */ | |
30 | /* */ | |
31 | /* Routines compute sin() and cos() as Double-Length numbers */ | |
32 | /********************************************************************/ | |
33 | ||
34 | ||
35 | ||
36 | #include "endian.h" | |
50944bca | 37 | #include "mydefs.h" |
e4d82761 UD |
38 | #include "sincos.tbl" |
39 | #include "dla.h" | |
40 | #include "dosincos.h" | |
15b3c029 AJ |
41 | #include "math_private.h" |
42 | ||
e4d82761 UD |
43 | /***********************************************************************/ |
44 | /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ | |
45 | /* as Double-Length number and store it at array v .It computes it by */ | |
46 | /* arithmetic action on Double-Length numbers */ | |
50944bca | 47 | /*(x+dx) between 0 and PI/4 */ |
e4d82761 UD |
48 | /***********************************************************************/ |
49 | ||
ca58f1db | 50 | void __dubsin(double x, double dx, double v[]) { |
50944bca | 51 | double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, |
e4d82761 | 52 | sn,ssn,cs,ccs,ds,dss,dc,dcc; |
50944bca UD |
53 | #if 0 |
54 | double xx,y,yy,z,zz; | |
55 | #endif | |
e4d82761 UD |
56 | mynumber u; |
57 | int4 k; | |
50944bca | 58 | |
e4d82761 UD |
59 | u.x=x+big.x; |
60 | k = u.i[LOW_HALF]<<2; | |
61 | x=x-(u.x-big.x); | |
62 | d=x+dx; | |
63 | dd=(x-d)+dx; | |
64 | /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ | |
65 | MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); | |
66 | sn=sincos.x[k]; /* */ | |
67 | ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ | |
68 | cs=sincos.x[k+2]; /* */ | |
69 | ccs=sincos.x[k+3]; /* */ | |
70 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */ | |
50944bca | 71 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); |
e4d82761 | 72 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */ |
50944bca | 73 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 UD |
74 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */ |
75 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
76 | ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */ | |
50944bca | 77 | |
e4d82761 | 78 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */ |
50944bca | 79 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); |
e4d82761 | 80 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */ |
50944bca | 81 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); |
e4d82761 | 82 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */ |
50944bca | 83 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); |
e4d82761 | 84 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */ |
50944bca | 85 | |
e4d82761 UD |
86 | MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); |
87 | MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
50944bca | 88 | SUB2(e,ee,dc,dcc,e,ee,r,s); |
e4d82761 | 89 | ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */ |
50944bca | 90 | |
e4d82761 UD |
91 | v[0]=e; |
92 | v[1]=ee; | |
93 | } | |
94 | /**********************************************************************/ | |
95 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ | |
96 | /* as Double-Length number and store it in array v .It computes it by */ | |
97 | /* arithmetic action on Double-Length numbers */ | |
98 | /*(x+dx) between 0 and PI/4 */ | |
99 | /**********************************************************************/ | |
100 | ||
ca58f1db | 101 | void __dubcos(double x, double dx, double v[]) { |
50944bca | 102 | double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, |
e4d82761 | 103 | sn,ssn,cs,ccs,ds,dss,dc,dcc; |
50944bca UD |
104 | #if 0 |
105 | double xx,y,yy,z,zz; | |
106 | #endif | |
e4d82761 | 107 | mynumber u; |
50944bca UD |
108 | int4 k; |
109 | u.x=x+big.x; | |
e4d82761 UD |
110 | k = u.i[LOW_HALF]<<2; |
111 | x=x-(u.x-big.x); | |
50944bca | 112 | d=x+dx; |
e4d82761 | 113 | dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ |
50944bca | 114 | MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); |
e4d82761 UD |
115 | sn=sincos.x[k]; /* */ |
116 | ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ | |
117 | cs=sincos.x[k+2]; /* */ | |
118 | ccs=sincos.x[k+3]; /* */ | |
119 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
50944bca | 120 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); |
e4d82761 | 121 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 122 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 UD |
123 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
124 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
50944bca UD |
125 | ADD2(ds,dss,d,dd,ds,dss,r,s); |
126 | ||
e4d82761 | 127 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 128 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); |
e4d82761 | 129 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 130 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); |
e4d82761 | 131 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 132 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); |
e4d82761 | 133 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
50944bca | 134 | |
e4d82761 UD |
135 | MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); |
136 | MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
137 | ||
138 | MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
139 | ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); | |
140 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); | |
50944bca | 141 | ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); |
e4d82761 | 142 | MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
50944bca UD |
143 | MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); |
144 | ADD2(ds,dss,d,dd,ds,dss,r,s); | |
e4d82761 UD |
145 | MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); |
146 | ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); | |
147 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
148 | ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); | |
149 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
150 | ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); | |
151 | MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
152 | MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); | |
153 | MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc); | |
50944bca UD |
154 | ADD2(e,ee,dc,dcc,e,ee,r,s); |
155 | SUB2(cs,ccs,e,ee,e,ee,r,s); | |
156 | ||
e4d82761 UD |
157 | v[0]=e; |
158 | v[1]=ee; | |
159 | } | |
160 | /**********************************************************************/ | |
161 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ | |
162 | /* as Double-Length number and store it in array v */ | |
163 | /**********************************************************************/ | |
ca58f1db | 164 | void __docos(double x, double dx, double v[]) { |
e4d82761 UD |
165 | double y,yy,p,w[2]; |
166 | if (x>0) {y=x; yy=dx;} | |
167 | else {y=-x; yy=-dx;} | |
168 | if (y<0.5*hp0.x) /* y< PI/4 */ | |
ca58f1db | 169 | {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];} |
e4d82761 UD |
170 | else if (y<1.5*hp0.x) { /* y< 3/4 * PI */ |
171 | p=hp0.x-y; /* p = PI/2 - y */ | |
172 | yy=hp1.x-yy; | |
173 | y=p+yy; | |
174 | yy=(p-y)+yy; | |
ca58f1db | 175 | if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} |
e4d82761 | 176 | /* cos(x) = sin ( 90 - x ) */ |
ca58f1db | 177 | else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; |
e4d82761 UD |
178 | } |
179 | } | |
180 | else { /* y>= 3/4 * PI */ | |
181 | p=2.0*hp0.x-y; /* p = PI- y */ | |
50944bca | 182 | yy=2.0*hp1.x-yy; |
e4d82761 UD |
183 | y=p+yy; |
184 | yy=(p-y)+yy; | |
ca58f1db | 185 | __dubcos(y,yy,w); |
e4d82761 UD |
186 | v[0]=-w[0]; |
187 | v[1]=-w[1]; | |
188 | } | |
189 | } |