]> git.ipfire.org Git - thirdparty/glibc.git/blame_incremental - sysdeps/ieee754/dbl-64/dosincos.c
projects / thirdparty / glibc.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Add script to update copyright notices and reformat some to facilitate its use.
[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / dosincos.c
... / ...
CommitLineData
1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2012 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 double r,s,c,cc,d,dd,d2,dd2,e,ee,
62 sn,ssn,cs,ccs,ds,dss,dc,dcc;
63#ifndef DLA_FMS
64 double p,hx,tx,hy,ty,q;
65#endif
66#if 0
67 double xx,y,yy,z,zz;
68#endif
69 mynumber u;
70 int4 k;
71
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;
77 /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
78 MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
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]; /* */
83 MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */
84 ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
85 MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */
86 ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
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) */
90
91 MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */
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); /* series */
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); /* for cos */
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); /* dc=cos(t) */
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 double r,s,c,cc,d,dd,d2,dd2,e,ee,
118 sn,ssn,cs,ccs,ds,dss,dc,dcc;
119#ifndef DLA_FMS
120 double p,hx,tx,hy,ty,q;
121#endif
122#if 0
123 double xx,y,yy,z,zz;
124#endif
125 mynumber u;
126 int4 k;
127 u.x=x+big.x;
128 k = u.i[LOW_HALF]<<2;
129 x=x-(u.x-big.x);
130 d=x+dx;
131 dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
132 MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
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]; /* */
137 MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
138 ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
139 MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
140 ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
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);
143 ADD2(ds,dss,d,dd,ds,dss,r,s);
144
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
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);
159 ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
160 MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
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);
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);
172 ADD2(e,ee,dc,dcc,e,ee,r,s);
173 SUB2(cs,ccs,e,ee,e,ee,r,s);
174
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/**********************************************************************/
182void
183SECTION
184__docos(double x, double dx, double v[]) {
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 */
189 {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
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;
195 if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
196 /* cos(x) = sin ( 90 - x ) */
197 else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
198 }
199 }
200 else { /* y>= 3/4 * PI */
201 p=2.0*hp0.x-y; /* p = PI- y */
202 yy=2.0*hp1.x-yy;
203 y=p+yy;
204 yy=(p-y)+yy;
205 __dubcos(y,yy,w);
206 v[0]=-w[0];
207 v[1]=-w[1];
208 }
209}
A mirror of the official glibc repository