2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
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
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /********************************************************************/
21 /* MODULE_NAME: dosincos.c */
24 /* FUNCTIONS: dubsin */
27 /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */
30 /* Routines compute sin() and cos() as Double-Length numbers */
31 /********************************************************************/
39 #include <math_private.h>
49 } __sincostab attribute_hidden
;
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 */
55 /*(x+dx) between 0 and PI/4 */
56 /***********************************************************************/
60 __dubsin (double x
, double dx
, double v
[])
62 double r
, s
, c
, cc
, d
, dd
, d2
, dd2
, e
, ee
,
63 sn
, ssn
, cs
, ccs
, ds
, dss
, dc
, dcc
;
65 double p
, hx
, tx
, hy
, ty
, q
;
71 k
= u
.i
[LOW_HALF
] << 2;
72 x
= x
- (u
.x
- big
.x
);
75 /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
76 MUL2 (d
, dd
, d
, dd
, d2
, dd2
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
77 sn
= __sincostab
.x
[k
]; /* */
78 ssn
= __sincostab
.x
[k
+ 1]; /* sin(Xi) and cos(Xi) */
79 cs
= __sincostab
.x
[k
+ 2]; /* */
80 ccs
= __sincostab
.x
[k
+ 3]; /* */
81 /* Taylor series for sin ds=sin(t) */
82 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
83 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
84 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
85 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
86 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
87 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
88 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
90 /* Taylor series for cos dc=cos(t) */
91 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
92 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
93 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
94 ADD2 (dc
, dcc
, c4
.x
, cc4
.x
, dc
, dcc
, r
, s
);
95 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
96 ADD2 (dc
, dcc
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
97 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
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
);
101 SUB2 (e
, ee
, dc
, dcc
, e
, ee
, r
, s
);
102 ADD2 (e
, ee
, sn
, ssn
, e
, ee
, r
, s
); /* e+ee=sin(x+dx) */
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 /**********************************************************************/
116 __dubcos (double x
, double dx
, double v
[])
118 double r
, s
, c
, cc
, d
, dd
, d2
, dd2
, e
, ee
,
119 sn
, ssn
, cs
, ccs
, ds
, dss
, dc
, dcc
;
121 double p
, hx
, tx
, hy
, ty
, q
;
126 k
= u
.i
[LOW_HALF
] << 2;
127 x
= x
- (u
.x
- big
.x
);
129 dd
= (x
- d
) + dx
; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
130 MUL2 (d
, dd
, d
, dd
, d2
, dd2
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
131 sn
= __sincostab
.x
[k
]; /* */
132 ssn
= __sincostab
.x
[k
+ 1]; /* sin(Xi) and cos(Xi) */
133 cs
= __sincostab
.x
[k
+ 2]; /* */
134 ccs
= __sincostab
.x
[k
+ 3]; /* */
135 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
136 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
137 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
138 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
139 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
140 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
141 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
143 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
144 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
145 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
146 ADD2 (dc
, dcc
, c4
.x
, cc4
.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
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
149 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
151 MUL2 (cs
, ccs
, ds
, dss
, e
, ee
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
152 MUL2 (dc
, dcc
, sn
, ssn
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
154 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
155 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
156 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
157 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
158 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
159 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
160 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
161 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
162 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
163 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
164 ADD2 (dc
, dcc
, c4
.x
, cc4
.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
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
167 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
168 MUL2 (sn
, ssn
, ds
, dss
, e
, ee
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
169 MUL2 (dc
, dcc
, cs
, ccs
, dc
, dcc
, p
, hx
, tx
, hy
, ty
, q
, c
, cc
);
170 ADD2 (e
, ee
, dc
, dcc
, e
, ee
, r
, s
);
171 SUB2 (cs
, ccs
, e
, ee
, e
, ee
, r
, s
);
176 /**********************************************************************/
177 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
178 /* as Double-Length number and store it in array v */
179 /**********************************************************************/
182 __docos (double x
, double dx
, double v
[])
184 double y
, yy
, p
, w
[2];
193 if (y
< 0.5 * hp0
.x
) /* y< PI/4 */
195 __dubcos (y
, yy
, w
); v
[0] = w
[0]; v
[1] = w
[1];
197 else if (y
< 1.5 * hp0
.x
) /* y< 3/4 * PI */
199 p
= hp0
.x
- y
; /* p = PI/2 - y */
205 __dubsin (y
, yy
, w
); v
[0] = w
[0]; v
[1] = w
[1];
207 /* cos(x) = sin ( 90 - x ) */
210 __dubsin (-y
, -yy
, w
); v
[0] = -w
[0]; v
[1] = -w
[1];
213 else /* y>= 3/4 * PI */
215 p
= 2.0 * hp0
.x
- y
; /* p = PI- y */
216 yy
= 2.0 * hp1
.x
- yy
;