2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 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 <https://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
;
68 k
= u
.i
[LOW_HALF
] << 2;
69 x
= x
- (u
.x
- big
.x
);
72 /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
73 MUL2 (d
, dd
, d
, dd
, d2
, dd2
, c
, cc
);
74 sn
= __sincostab
.x
[k
]; /* */
75 ssn
= __sincostab
.x
[k
+ 1]; /* sin(Xi) and cos(Xi) */
76 cs
= __sincostab
.x
[k
+ 2]; /* */
77 ccs
= __sincostab
.x
[k
+ 3]; /* */
78 /* Taylor series for sin ds=sin(t) */
79 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, c
, cc
);
80 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
81 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
82 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
83 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
84 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, c
, cc
);
85 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
87 /* Taylor series for cos dc=cos(t) */
88 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, c
, cc
);
89 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
90 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
91 ADD2 (dc
, dcc
, c4
.x
, cc4
.x
, dc
, dcc
, r
, s
);
92 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
93 ADD2 (dc
, dcc
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
94 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
96 MUL2 (cs
, ccs
, ds
, dss
, e
, ee
, c
, cc
);
97 MUL2 (dc
, dcc
, sn
, ssn
, dc
, dcc
, c
, cc
);
98 SUB2 (e
, ee
, dc
, dcc
, e
, ee
, r
, s
);
99 ADD2 (e
, ee
, sn
, ssn
, e
, ee
, r
, s
); /* e+ee=sin(x+dx) */
104 /**********************************************************************/
105 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
106 /* as Double-Length number and store it in array v .It computes it by */
107 /* arithmetic action on Double-Length numbers */
108 /*(x+dx) between 0 and PI/4 */
109 /**********************************************************************/
113 __dubcos (double x
, double dx
, double v
[])
115 double r
, s
, c
, cc
, d
, dd
, d2
, dd2
, e
, ee
,
116 sn
, ssn
, cs
, ccs
, ds
, dss
, dc
, dcc
;
120 k
= u
.i
[LOW_HALF
] << 2;
121 x
= x
- (u
.x
- big
.x
);
123 dd
= (x
- d
) + dx
; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
124 MUL2 (d
, dd
, d
, dd
, d2
, dd2
, c
, cc
);
125 sn
= __sincostab
.x
[k
]; /* */
126 ssn
= __sincostab
.x
[k
+ 1]; /* sin(Xi) and cos(Xi) */
127 cs
= __sincostab
.x
[k
+ 2]; /* */
128 ccs
= __sincostab
.x
[k
+ 3]; /* */
129 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, c
, cc
);
130 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
131 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
132 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
133 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
134 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, c
, cc
);
135 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
137 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, c
, cc
);
138 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
139 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
140 ADD2 (dc
, dcc
, c4
.x
, cc4
.x
, dc
, dcc
, r
, s
);
141 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
142 ADD2 (dc
, dcc
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
143 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
145 MUL2 (cs
, ccs
, ds
, dss
, e
, ee
, c
, cc
);
146 MUL2 (dc
, dcc
, sn
, ssn
, dc
, dcc
, c
, cc
);
148 MUL2 (d2
, dd2
, s7
.x
, ss7
.x
, ds
, dss
, c
, cc
);
149 ADD2 (ds
, dss
, s5
.x
, ss5
.x
, ds
, dss
, r
, s
);
150 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
151 ADD2 (ds
, dss
, s3
.x
, ss3
.x
, ds
, dss
, r
, s
);
152 MUL2 (d2
, dd2
, ds
, dss
, ds
, dss
, c
, cc
);
153 MUL2 (d
, dd
, ds
, dss
, ds
, dss
, c
, cc
);
154 ADD2 (ds
, dss
, d
, dd
, ds
, dss
, r
, s
);
155 MUL2 (d2
, dd2
, c8
.x
, cc8
.x
, dc
, dcc
, c
, cc
);
156 ADD2 (dc
, dcc
, c6
.x
, cc6
.x
, dc
, dcc
, r
, s
);
157 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
158 ADD2 (dc
, dcc
, c4
.x
, cc4
.x
, dc
, dcc
, r
, s
);
159 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
160 ADD2 (dc
, dcc
, c2
.x
, cc2
.x
, dc
, dcc
, r
, s
);
161 MUL2 (d2
, dd2
, dc
, dcc
, dc
, dcc
, c
, cc
);
162 MUL2 (sn
, ssn
, ds
, dss
, e
, ee
, c
, cc
);
163 MUL2 (dc
, dcc
, cs
, ccs
, dc
, dcc
, c
, cc
);
164 ADD2 (e
, ee
, dc
, dcc
, e
, ee
, r
, s
);
165 SUB2 (cs
, ccs
, e
, ee
, e
, ee
, r
, s
);
170 /**********************************************************************/
171 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
172 /* as Double-Length number and store it in array v */
173 /**********************************************************************/
176 __docos (double x
, double dx
, double v
[])
178 double y
, yy
, p
, w
[2];
187 if (y
< 0.5 * hp0
.x
) /* y< PI/4 */
189 __dubcos (y
, yy
, w
); v
[0] = w
[0]; v
[1] = w
[1];
191 else if (y
< 1.5 * hp0
.x
) /* y< 3/4 * PI */
193 p
= hp0
.x
- y
; /* p = PI/2 - y */
199 __dubsin (y
, yy
, w
); v
[0] = w
[0]; v
[1] = w
[1];
201 /* cos(x) = sin ( 90 - x ) */
204 __dubsin (-y
, -yy
, w
); v
[0] = -w
[0]; v
[1] = -w
[1];
207 else /* y>= 3/4 * PI */
209 p
= 2.0 * hp0
.x
- y
; /* p = PI- y */
210 yy
= 2.0 * hp1
.x
- yy
;