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