1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 */
20/********************************************************************/
21/*                                                                  */
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"
37#include "mydefs.h"
38#include "sincos.tbl"
39#include "dla.h"
40#include "dosincos.h"
41#include "math_private.h"
42
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                          */
47/*(x+dx) between 0 and PI/4                                            */
48/***********************************************************************/
49
50void __dubsin(double x, double dx, double v[]) {
51  double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
52    sn,ssn,cs,ccs,ds,dss,dc,dcc;
53#if 0
54  double xx,y,yy,z,zz;
55#endif
56  mynumber u;
57  int4 k;
58
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    */
71  ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
72  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* series    */
73  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
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) */
77
78  MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor    */
79  ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
80  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* series    */
81  ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
82  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* for cos   */
83  ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
84  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* dc=cos(t) */
85
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);
88  SUB2(e,ee,dc,dcc,e,ee,r,s);
89  ADD2(e,ee,sn,ssn,e,ee,r,s);                    /* e+ee=sin(x+dx) */
90
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
101void __dubcos(double x, double dx, double v[]) {
102  double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
103    sn,ssn,cs,ccs,ds,dss,dc,dcc;
104#if 0
105  double xx,y,yy,z,zz;
106#endif
107  mynumber u;
108  int4 k;
109  u.x=x+big.x;
110  k = u.i[LOW_HALF]<<2;
111  x=x-(u.x-big.x);
112  d=x+dx;
113  dd=(x-d)+dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
114  MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
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);
120  ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
121  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
122  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
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);
125  ADD2(ds,dss,d,dd,ds,dss,r,s);
126
127  MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
128  ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
129  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
130  ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
131  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
132  ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
133  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
134
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);
141  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
142  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
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);
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);
154  ADD2(e,ee,dc,dcc,e,ee,r,s);
155  SUB2(cs,ccs,e,ee,e,ee,r,s);
156
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/**********************************************************************/
164void __docos(double x, double dx, double v[]) {
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    */
169           {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
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;
175       if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
176                                       /* cos(x) = sin ( 90 -  x ) */
177         else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
178	 }
179     }
180  else { /* y>= 3/4 * PI */
181    p=2.0*hp0.x-y;    /* p = PI- y */
182    yy=2.0*hp1.x-yy;
183    y=p+yy;
184    yy=(p-y)+yy;
185    __dubcos(y,yy,w);
186    v[0]=-w[0];
187    v[1]=-w[1];
188  }
189}
190