1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2019 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
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.
15 *
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/>.
18 */
19/********************************************************************/
20/* */
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"
36#include "mydefs.h"
37#include <dla.h>
38#include "dosincos.h"
39#include <math_private.h>
40
41#ifndef SECTION
42# define SECTION
43#endif
44
45extern const union
46{
47 int4 i[880];
48 double x[440];
49} __sincostab attribute_hidden;
50
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/***********************************************************************/
57
58void
59SECTION
60__dubsin (double x, double dx, double v[])
61{
62 double r, s, c, cc, d, dd, d2, dd2, e, ee,
63 sn, ssn, cs, ccs, ds, dss, dc, dcc;
64#ifndef DLA_FMS
65 double p, hx, tx, hy, ty, q;
66#endif
67 mynumber u;
68 int4 k;
69
70 u.x = x + big.x;
71 k = u.i[LOW_HALF] << 2;
72 x = x - (u.x - big.x);
73 d = x + dx;
74 dd = (x - d) + dx;
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);
89
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);
98
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) */
103
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
114void
115SECTION
116__dubcos (double x, double dx, double v[])
117{
118 double r, s, c, cc, d, dd, d2, dd2, e, ee,
119 sn, ssn, cs, ccs, ds, dss, dc, dcc;
120#ifndef DLA_FMS
121 double p, hx, tx, hy, ty, q;
122#endif
123 mynumber u;
124 int4 k;
125 u.x = x + big.x;
126 k = u.i[LOW_HALF] << 2;
127 x = x - (u.x - big.x);
128 d = x + dx;
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);
142
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);
150
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);
153
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);
172
173 v[0] = e;
174 v[1] = ee;
175}
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/**********************************************************************/
180void
181SECTION
182__docos (double x, double dx, double v[])
183{
184 double y, yy, p, w[2];
185 if (x > 0)
186 {
187 y = x; yy = dx;
188 }
189 else
190 {
191 y = -x; yy = -dx;
192 }
193 if (y < 0.5 * hp0.x) /* y< PI/4 */
194 {
195 __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
196 }
197 else if (y < 1.5 * hp0.x) /* y< 3/4 * PI */
198 {
199 p = hp0.x - y; /* p = PI/2 - y */
200 yy = hp1.x - yy;
201 y = p + yy;
202 yy = (p - y) + yy;
203 if (y > 0)
204 {
205 __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
206 }
207 /* cos(x) = sin ( 90 - x ) */
208 else
209 {
210 __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
211 }
212 }
213 else /* y>= 3/4 * PI */
214 {
215 p = 2.0 * hp0.x - y; /* p = PI- y */
216 yy = 2.0 * hp1.x - yy;
217 y = p + yy;
218 yy = (p - y) + yy;
219 __dubcos (y, yy, w);
220 v[0] = -w[0];
221 v[1] = -w[1];
222 }
223}
224