e_rem_pio2l.h revision 8870
1227825Stheraven/* @(#)e_rem_pio2.c 5.1 93/09/24 */ 2227825Stheraven/* 3227825Stheraven * ==================================================== 4227825Stheraven * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5227825Stheraven * 6227825Stheraven * Developed at SunPro, a Sun Microsystems, Inc. business. 7227825Stheraven * Permission to use, copy, modify, and distribute this 8227825Stheraven * software is freely granted, provided that this notice 9227825Stheraven * is preserved. 10227825Stheraven * ==================================================== 11227825Stheraven */ 12227825Stheraven 13227825Stheraven#ifndef lint 14227825Stheravenstatic char rcsid[] = "$Id: e_rem_pio2.c,v 1.2 1995/04/07 23:23:24 bde Exp $"; 15227825Stheraven#endif 16300770Sdim 17227825Stheraven/* __ieee754_rem_pio2(x,y) 18300770Sdim * 19300770Sdim * return the remainder of x rem pi/2 in y[0]+y[1] 20227825Stheraven * use __kernel_rem_pio2() 21300770Sdim */ 22300770Sdim 23300770Sdim#include "math.h" 24300770Sdim#include "math_private.h" 25227825Stheraven 26227825Stheraven/* 27300770Sdim * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 28300770Sdim */ 29300770Sdim#ifdef __STDC__ 30300770Sdimstatic const int32_t two_over_pi[] = { 31227825Stheraven#else 32227825Stheravenstatic int32_t two_over_pi[] = { 33300770Sdim#endif 34300770Sdim0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 35300770Sdim0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 36300770Sdim0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 37300770Sdim0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 38300770Sdim0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 39227825Stheraven0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 40227825Stheraven0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 41300770Sdim0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 42300770Sdim0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 43300770Sdim0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 44300770Sdim0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 45300770Sdim}; 46300770Sdim 47227825Stheraven#ifdef __STDC__ 48227825Stheravenstatic const int32_t npio2_hw[] = { 49300770Sdim#else 50300770Sdimstatic int32_t npio2_hw[] = { 51300770Sdim#endif 52300770Sdim0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 53300770Sdim0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 54300770Sdim0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 55227825Stheraven0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 56227825Stheraven0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 57300770Sdim0x404858EB, 0x404921FB, 58300770Sdim}; 59300770Sdim 60300770Sdim/* 61300770Sdim * invpio2: 53 bits of 2/pi 62300770Sdim * pio2_1: first 33 bit of pi/2 63300770Sdim * pio2_1t: pi/2 - pio2_1 64300770Sdim * pio2_2: second 33 bit of pi/2 65227825Stheraven * pio2_2t: pi/2 - (pio2_1+pio2_2) 66227825Stheraven * pio2_3: third 33 bit of pi/2 67227825Stheraven * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 68300770Sdim */ 69227825Stheraven 70227825Stheraven#ifdef __STDC__ 71227825Stheravenstatic const double 72300770Sdim#else 73227825Stheravenstatic double 74300770Sdim#endif 75300770Sdimzero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 76227825Stheravenhalf = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 77227825Stheraventwo24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 78227825Stheraveninvpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 79300770Sdimpio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 80227825Stheravenpio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 81300770Sdimpio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 82300770Sdimpio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 83227825Stheravenpio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 84227825Stheravenpio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 85227825Stheraven 86300770Sdim#ifdef __STDC__ 87227825Stheraven int32_t __ieee754_rem_pio2(double x, double *y) 88300770Sdim#else 89300770Sdim int32_t __ieee754_rem_pio2(x,y) 90227825Stheraven double x,y[]; 91227825Stheraven#endif 92227825Stheraven{ 93300770Sdim double z,w,t,r,fn; 94227825Stheraven double tx[3]; 95300770Sdim int32_t e0,i,j,nx,n,ix,hx; 96300770Sdim u_int32_t low; 97227825Stheraven 98227825Stheraven GET_HIGH_WORD(hx,x); /* high word of x */ 99227825Stheraven ix = hx&0x7fffffff; 100300770Sdim if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 101227825Stheraven {y[0] = x; y[1] = 0; return 0;} 102300770Sdim if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 103300770Sdim if(hx>0) { 104227825Stheraven z = x - pio2_1; 105227825Stheraven if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 106227825Stheraven y[0] = z - pio2_1t; 107300770Sdim y[1] = (z-y[0])-pio2_1t; 108227825Stheraven } else { /* near pi/2, use 33+33+53 bit pi */ 109300770Sdim z -= pio2_2; 110300770Sdim y[0] = z - pio2_2t; 111227825Stheraven y[1] = (z-y[0])-pio2_2t; 112227825Stheraven } 113227825Stheraven return 1; 114300770Sdim } else { /* negative x */ 115227825Stheraven z = x + pio2_1; 116300770Sdim if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 117300770Sdim y[0] = z + pio2_1t; 118227825Stheraven y[1] = (z-y[0])+pio2_1t; 119227825Stheraven } else { /* near pi/2, use 33+33+53 bit pi */ 120227825Stheraven z += pio2_2; 121300770Sdim y[0] = z + pio2_2t; 122227825Stheraven y[1] = (z-y[0])+pio2_2t; 123300770Sdim } 124300770Sdim return -1; 125227825Stheraven } 126227825Stheraven } 127227825Stheraven if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 128300770Sdim t = fabs(x); 129227825Stheraven n = (int32_t) (t*invpio2+half); 130300770Sdim fn = (double)n; 131300770Sdim r = t-fn*pio2_1; 132227825Stheraven w = fn*pio2_1t; /* 1st round good to 85 bit */ 133227825Stheraven if(n<32&&ix!=npio2_hw[n-1]) { 134227825Stheraven y[0] = r-w; /* quick check no cancellation */ 135300770Sdim } else { 136227825Stheraven u_int32_t high; 137300770Sdim j = ix>>20; 138300770Sdim y[0] = r-w; 139227825Stheraven GET_HIGH_WORD(high,y[0]); 140227825Stheraven i = j-((high>>20)&0x7ff); 141227825Stheraven if(i>16) { /* 2nd iteration needed, good to 118 */ 142227825Stheraven t = r; 143227825Stheraven w = fn*pio2_2; 144232924Stheraven r = t-w; 145227825Stheraven w = fn*pio2_2t-((t-r)-w); 146300770Sdim y[0] = r-w; 147300770Sdim GET_HIGH_WORD(high,y[0]); 148227825Stheraven i = j-((high>>20)&0x7ff); 149227825Stheraven if(i>49) { /* 3rd iteration need, 151 bits acc */ 150227825Stheraven t = r; /* will cover all possible cases */ 151227825Stheraven w = fn*pio2_3; 152232924Stheraven r = t-w; 153227825Stheraven w = fn*pio2_3t-((t-r)-w); 154300770Sdim y[0] = r-w; 155300770Sdim } 156227825Stheraven } 157227825Stheraven } 158227825Stheraven y[1] = (r-y[0])-w; 159227825Stheraven if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 160232924Stheraven else return n; 161227825Stheraven } 162300770Sdim /* 163300770Sdim * all other (large) arguments 164227825Stheraven */ 165227825Stheraven if(ix>=0x7ff00000) { /* x is inf or NaN */ 166227825Stheraven y[0]=y[1]=x-x; return 0; 167227825Stheraven } 168232924Stheraven /* set z = scalbn(|x|,ilogb(x)-23) */ 169227825Stheraven GET_LOW_WORD(low,x); 170300770Sdim SET_LOW_WORD(z,low); 171300770Sdim e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 172227825Stheraven SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 173227825Stheraven for(i=0;i<2;i++) { 174227825Stheraven tx[i] = (double)((int32_t)(z)); 175227825Stheraven z = (z-tx[i])*two24; 176232924Stheraven } 177227825Stheraven tx[2] = z; 178227825Stheraven nx = 3; 179232924Stheraven while(tx[nx-1]==zero) nx--; /* skip zero term */ 180227825Stheraven n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 181227825Stheraven if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 182232924Stheraven return n; 183232924Stheraven} 184227825Stheraven