00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
00015 #endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 #include "math.h"
00134 #include "mathP.h"
00135
00136 #ifdef __STDC__
00137 static const int init_jk[] = {2,3,4,6};
00138 #else
00139 static int init_jk[] = {2,3,4,6};
00140 #endif
00141
00142 #ifdef __STDC__
00143 static const double PIo2[] = {
00144 #else
00145 static double PIo2[] = {
00146 #endif
00147 1.57079625129699707031e+00,
00148 7.54978941586159635335e-08,
00149 5.39030252995776476554e-15,
00150 3.28200341580791294123e-22,
00151 1.27065575308067607349e-29,
00152 1.22933308981111328932e-36,
00153 2.73370053816464559624e-44,
00154 2.16741683877804819444e-51,
00155 };
00156
00157 #ifdef __STDC__
00158 static const double
00159 #else
00160 static double
00161 #endif
00162 zero = 0.0,
00163 one = 1.0,
00164 two24 = 1.67772160000000000000e+07,
00165 twon24 = 5.96046447753906250000e-08;
00166
00167 #ifdef __STDC__
00168 int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
00169 #else
00170 int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
00171 double x[], y[]; int e0,nx,prec; int32_t ipio2[];
00172 #endif
00173 {
00174 int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
00175 double z,fw,f[20],fq[20],q[20];
00176
00177
00178 jk = init_jk[prec];
00179 jp = jk;
00180
00181
00182 jx = nx-1;
00183 jv = (e0-3)/24; if(jv<0) jv=0;
00184 q0 = e0-24*(jv+1);
00185
00186
00187 j = jv-jx; m = jx+jk;
00188 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
00189
00190
00191 for (i=0;i<=jk;i++) {
00192 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
00193 }
00194
00195 jz = jk;
00196 recompute:
00197
00198 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
00199 fw = (double)((int32_t)(twon24* z));
00200 iq[i] = (int32_t)(z-two24*fw);
00201 z = q[j-1]+fw;
00202 }
00203
00204
00205 z = scalbn(z,q0);
00206 z -= 8.0*floor(z*0.125);
00207 n = (int32_t) z;
00208 z -= (double)n;
00209 ih = 0;
00210 if(q0>0) {
00211 i = (iq[jz-1]>>(24-q0)); n += i;
00212 iq[jz-1] -= i<<(24-q0);
00213 ih = iq[jz-1]>>(23-q0);
00214 }
00215 else if(q0==0) ih = iq[jz-1]>>23;
00216 else if(z>=0.5) ih=2;
00217
00218 if(ih>0) {
00219 n += 1; carry = 0;
00220 for(i=0;i<jz ;i++) {
00221 j = iq[i];
00222 if(carry==0) {
00223 if(j!=0) {
00224 carry = 1; iq[i] = 0x1000000- j;
00225 }
00226 } else iq[i] = 0xffffff - j;
00227 }
00228 if(q0>0) {
00229 switch(q0) {
00230 case 1:
00231 iq[jz-1] &= 0x7fffff; break;
00232 case 2:
00233 iq[jz-1] &= 0x3fffff; break;
00234 }
00235 }
00236 if(ih==2) {
00237 z = one - z;
00238 if(carry!=0) z -= scalbn(one,q0);
00239 }
00240 }
00241
00242
00243 if(z==zero) {
00244 j = 0;
00245 for (i=jz-1;i>=jk;i--) j |= iq[i];
00246 if(j==0) {
00247 for(k=1;iq[jk-k]==0;k++);
00248
00249 for(i=jz+1;i<=jz+k;i++) {
00250 f[jx+i] = (double) ipio2[jv+i];
00251 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
00252 q[i] = fw;
00253 }
00254 jz += k;
00255 goto recompute;
00256 }
00257 }
00258
00259
00260 if(z==0.0) {
00261 jz -= 1; q0 -= 24;
00262 while(iq[jz]==0) { jz--; q0-=24;}
00263 } else {
00264 z = scalbn(z,-q0);
00265 if(z>=two24) {
00266 fw = (double)((int32_t)(twon24*z));
00267 iq[jz] = (int32_t)(z-two24*fw);
00268 jz += 1; q0 += 24;
00269 iq[jz] = (int32_t) fw;
00270 } else iq[jz] = (int32_t) z ;
00271 }
00272
00273
00274 fw = scalbn(one,q0);
00275 for(i=jz;i>=0;i--) {
00276 q[i] = fw*(double)iq[i]; fw*=twon24;
00277 }
00278
00279
00280 for(i=jz;i>=0;i--) {
00281 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
00282 fq[jz-i] = fw;
00283 }
00284
00285
00286 switch(prec) {
00287 case 0:
00288 fw = 0.0;
00289 for (i=jz;i>=0;i--) fw += fq[i];
00290 y[0] = (ih==0)? fw: -fw;
00291 break;
00292 case 1:
00293 case 2:
00294 fw = 0.0;
00295 for (i=jz;i>=0;i--) fw += fq[i];
00296 y[0] = (ih==0)? fw: -fw;
00297 fw = fq[0]-fw;
00298 for (i=1;i<=jz;i++) fw += fq[i];
00299 y[1] = (ih==0)? fw: -fw;
00300 break;
00301 case 3:
00302 for (i=jz;i>0;i--) {
00303 fw = fq[i-1]+fq[i];
00304 fq[i] += fq[i-1]-fw;
00305 fq[i-1] = fw;
00306 }
00307 for (i=jz;i>1;i--) {
00308 fw = fq[i-1]+fq[i];
00309 fq[i] += fq[i-1]-fw;
00310 fq[i-1] = fw;
00311 }
00312 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
00313 if(ih==0) {
00314 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
00315 } else {
00316 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
00317 }
00318 }
00319 return n&7;
00320 }