e_rem_pio2l.h revision 176385
139220Sgibbs 239220Sgibbs/* @(#)e_rem_pio2.c 1.4 95/01/18 */ 339220Sgibbs/* 439220Sgibbs * ==================================================== 539220Sgibbs * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 639220Sgibbs * 739220Sgibbs * Developed at SunSoft, a Sun Microsystems, Inc. business. 839220Sgibbs * Permission to use, copy, modify, and distribute this 939220Sgibbs * software is freely granted, provided that this notice 1039220Sgibbs * is preserved. 1139220Sgibbs * ==================================================== 1239220Sgibbs * 1339220Sgibbs */ 1439220Sgibbs 1539220Sgibbs#include <sys/cdefs.h> 1639220Sgibbs__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2.c 176385 2008-02-18 14:02:12Z bde $"); 1739220Sgibbs 1839220Sgibbs/* __ieee754_rem_pio2(x,y) 1939220Sgibbs * 2039220Sgibbs * return the remainder of x rem pi/2 in y[0]+y[1] 2150477Speter * use __kernel_rem_pio2() 2239220Sgibbs */ 2339220Sgibbs 2439220Sgibbs#include "math.h" 2539220Sgibbs#include "math_private.h" 2639220Sgibbs 2739220Sgibbsstatic const int32_t npio2_hw[] = { 2839220Sgibbs0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 2939220Sgibbs0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 3039220Sgibbs0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 3139220Sgibbs0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 3239220Sgibbs0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 3339220Sgibbs0x404858EB, 0x404921FB, 3439220Sgibbs}; 3539220Sgibbs 3639220Sgibbs/* 3739220Sgibbs * invpio2: 53 bits of 2/pi 3839220Sgibbs * pio2_1: first 33 bit of pi/2 3939220Sgibbs * pio2_1t: pi/2 - pio2_1 4039220Sgibbs * pio2_2: second 33 bit of pi/2 4139220Sgibbs * pio2_2t: pi/2 - (pio2_1+pio2_2) 4239220Sgibbs * pio2_3: third 33 bit of pi/2 4339220Sgibbs * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 4439220Sgibbs */ 4539220Sgibbs 4639220Sgibbsstatic const double 4739220Sgibbszero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 4839220Sgibbshalf = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 4939220Sgibbstwo24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 5039220Sgibbsinvpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 5139220Sgibbspio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 5239220Sgibbspio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 5339220Sgibbspio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 5439220Sgibbspio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 5539220Sgibbspio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 5639220Sgibbspio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 5739220Sgibbs 5847507Sgibbs#ifdef INLINE_REM_PIO2 5947507Sgibbsextern inline 6039220Sgibbs#endif 6139220Sgibbsint 6239220Sgibbs__ieee754_rem_pio2(double x, double *y) 6339220Sgibbs{ 6439220Sgibbs double z,w,t,r,fn; 6539220Sgibbs double tx[3]; 6639220Sgibbs int32_t e0,i,j,nx,n,ix,hx; 6739220Sgibbs u_int32_t low; 6839220Sgibbs 6939220Sgibbs GET_HIGH_WORD(hx,x); /* high word of x */ 7039220Sgibbs ix = hx&0x7fffffff; 7139220Sgibbs if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 7239220Sgibbs {y[0] = x; y[1] = 0; return 0;} 7339220Sgibbs if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 7439220Sgibbs if(hx>0) { 7539220Sgibbs z = x - pio2_1; 7639220Sgibbs if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 7739220Sgibbs y[0] = z - pio2_1t; 7839220Sgibbs y[1] = (z-y[0])-pio2_1t; 7939220Sgibbs } else { /* near pi/2, use 33+33+53 bit pi */ 8039220Sgibbs z -= pio2_2; 8139220Sgibbs y[0] = z - pio2_2t; 8239220Sgibbs y[1] = (z-y[0])-pio2_2t; 8339220Sgibbs } 8439220Sgibbs return 1; 8539220Sgibbs } else { /* negative x */ 8639220Sgibbs z = x + pio2_1; 8739220Sgibbs if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 8839220Sgibbs y[0] = z + pio2_1t; 8939220Sgibbs y[1] = (z-y[0])+pio2_1t; 9039220Sgibbs } else { /* near pi/2, use 33+33+53 bit pi */ 9139220Sgibbs z += pio2_2; 9263457Sgibbs y[0] = z + pio2_2t; 9339220Sgibbs y[1] = (z-y[0])+pio2_2t; 9439220Sgibbs } 9539220Sgibbs return -1; 9639220Sgibbs } 9739220Sgibbs } 9863457Sgibbs if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 9963457Sgibbs t = fabs(x); 10039220Sgibbs n = (int32_t) (t*invpio2+half); 10139220Sgibbs fn = (double)n; 10239220Sgibbs r = t-fn*pio2_1; 10339220Sgibbs w = fn*pio2_1t; /* 1st round good to 85 bit */ 10439220Sgibbs if(n<32&&ix!=npio2_hw[n-1]) { 10539220Sgibbs y[0] = r-w; /* quick check no cancellation */ 10639220Sgibbs } else { 10739220Sgibbs u_int32_t high; 10839220Sgibbs j = ix>>20; 10939220Sgibbs y[0] = r-w; 11039220Sgibbs GET_HIGH_WORD(high,y[0]); 11139220Sgibbs i = j-((high>>20)&0x7ff); 11239220Sgibbs if(i>16) { /* 2nd iteration needed, good to 118 */ 11339220Sgibbs t = r; 11439220Sgibbs w = fn*pio2_2; 11539220Sgibbs r = t-w; 11639220Sgibbs w = fn*pio2_2t-((t-r)-w); 11739220Sgibbs y[0] = r-w; 11839220Sgibbs GET_HIGH_WORD(high,y[0]); 11939220Sgibbs i = j-((high>>20)&0x7ff); 12039220Sgibbs if(i>49) { /* 3rd iteration need, 151 bits acc */ 12139220Sgibbs t = r; /* will cover all possible cases */ 12239220Sgibbs w = fn*pio2_3; 12339220Sgibbs r = t-w; 12439220Sgibbs w = fn*pio2_3t-((t-r)-w); 12539220Sgibbs y[0] = r-w; 12639220Sgibbs } 12739220Sgibbs } 12839220Sgibbs } 12939220Sgibbs y[1] = (r-y[0])-w; 13039220Sgibbs if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 13139220Sgibbs else return n; 13239220Sgibbs } 13339220Sgibbs /* 13439220Sgibbs * all other (large) arguments 13539220Sgibbs */ 13639220Sgibbs if(ix>=0x7ff00000) { /* x is inf or NaN */ 13739220Sgibbs y[0]=y[1]=x-x; return 0; 13839220Sgibbs } 13939220Sgibbs /* set z = scalbn(|x|,ilogb(x)-23) */ 14039220Sgibbs GET_LOW_WORD(low,x); 14139220Sgibbs SET_LOW_WORD(z,low); 14239220Sgibbs e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 14339220Sgibbs SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 14439220Sgibbs for(i=0;i<2;i++) { 14539220Sgibbs tx[i] = (double)((int32_t)(z)); 14639220Sgibbs z = (z-tx[i])*two24; 14739220Sgibbs } 14839220Sgibbs tx[2] = z; 14939220Sgibbs nx = 3; 15039220Sgibbs while(tx[nx-1]==zero) nx--; /* skip zero term */ 15139220Sgibbs n = __kernel_rem_pio2(tx,y,e0,nx,2); 15239220Sgibbs if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 15339220Sgibbs return n; 15439220Sgibbs} 15539220Sgibbs