e_rem_pio2l.h revision 2117
150477Speter/* @(#)e_rem_pio2.c 5.1 93/09/24 */ 212029Sache/* 3127474Stjr * ==================================================== 4127474Stjr * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5123682Sache * 6118459Smtm * Developed at SunPro, a Sun Microsystems, Inc. business. 786072Sache * Permission to use, copy, modify, and distribute this 888473Sphantom * software is freely granted, provided that this notice 9117259Sache * is preserved. 1088473Sphantom * ==================================================== 1177977Sache */ 12125208Sache 13118652Sache#ifndef lint 1477977Sachestatic char rcsid[] = "$Id: e_rem_pio2.c,v 1.5 1994/08/18 23:05:56 jtc Exp $"; 1577977Sache#endif 16196790Sache 1788473Sphantom/* __ieee754_rem_pio2(x,y) 1877977Sache * 1977977Sache * return the remainder of x rem pi/2 in y[0]+y[1] 2088473Sphantom * use __kernel_rem_pio2() 2170484Sphantom */ 2277977Sache 2370484Sphantom#include "math.h" 24174904Sache#include "math_private.h" 2552388Sache 2677977Sache/* 27118174Sache * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 28122151Sdavidxu */ 29115628Sache#ifdef __STDC__ 30123657Sachestatic const int32_t two_over_pi[] = { 3147831Sfoxfair#else 3223222Swoschstatic int32_t two_over_pi[] = { 33291794Sbdrewery#endif 3412029Sache0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 35136611Sru0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 36136611Sru0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 37136611Sru0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 38136611Sru0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 39136611Sru0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 40136611Sru0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 41136611Sru0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 42136611Sru0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 43136611Sru0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 44136611Sru0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 45136611Sru}; 46136611Sru 47136611Sru#ifdef __STDC__ 48136611Srustatic const int32_t npio2_hw[] = { 49174904Sache#else 50196790Sachestatic int32_t npio2_hw[] = { 51136611Sru#endif 5288473Sphantom0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 53136611Sru0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 54136611Sru0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 55193908Sedwin0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 56193908Sedwin0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 57193908Sedwin0x404858EB, 0x404921FB, 58164131Sdes}; 5993885Sphantom 60136611Sru/* 61136611Sru * invpio2: 53 bits of 2/pi 6293885Sphantom * pio2_1: first 33 bit of pi/2 63136611Sru * pio2_1t: pi/2 - pio2_1 64136611Sru * pio2_2: second 33 bit of pi/2 65105445Sache * pio2_2t: pi/2 - (pio2_1+pio2_2) 66136611Sru * pio2_3: third 33 bit of pi/2 67136611Sru * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 6888473Sphantom */ 69136611Sru 70136611Sru#ifdef __STDC__ 71123682Sachestatic const double 72136611Sru#else 73136611Srustatic double 74196790Sache#endif 75196790Sachezero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 76196790Sachehalf = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 77143126Srutwo24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 78136611Sruinvpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 79136611Srupio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 80193908Sedwinpio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 81193908Sedwinpio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 82134437Stjrpio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 83134437Stjrpio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 84128526Stjrpio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 85196790Sache 86164131Sdes#ifdef __STDC__ 87128526Stjr int32_t __ieee754_rem_pio2(double x, double *y) 88128526Stjr#else 89136611Sru int32_t __ieee754_rem_pio2(x,y) 9078002Sache double x,y[]; 9112029Sache#endif 92136611Sru{ 93136611Sru double z,w,t,r,fn; 94136611Sru double tx[3]; 95136611Sru int32_t e0,i,j,nx,n,ix,hx; 96136611Sru u_int32_t low; 97136611Sru 9823222Swosch GET_HIGH_WORD(hx,x); /* high word of x */ 9923222Swosch ix = hx&0x7fffffff; 10041760Sdillon if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 10112029Sache {y[0] = x; y[1] = 0; return 0;} 10212029Sache if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 103 t = fabs(x); 104 n = (int32_t) (t*invpio2+half); 105 fn = (double)n; 106 r = t-fn*pio2_1; 107 w = fn*pio2_1t; /* 1st round good to 85 bit */ 108 if(n<32&&ix!=npio2_hw[n-1]) { 109 y[0] = r-w; /* quick check no cancellation */ 110 } else { 111 u_int32_t high; 112 j = ix>>20; 113 y[0] = r-w; 114 GET_HIGH_WORD(high,y[0]); 115 i = j-((high>>20)&0x7ff); 116 if(i>16) { /* 2nd iteration needed, good to 118 */ 117 t = r; 118 w = fn*pio2_2; 119 r = t-w; 120 w = fn*pio2_2t-((t-r)-w); 121 y[0] = r-w; 122 GET_HIGH_WORD(high,y[0]); 123 i = j-((high>>20)&0x7ff); 124 if(i>49) { /* 3rd iteration need, 151 bits acc */ 125 t = r; /* will cover all possible cases */ 126 w = fn*pio2_3; 127 r = t-w; 128 w = fn*pio2_3t-((t-r)-w); 129 y[0] = r-w; 130 } 131 } 132 } 133 y[1] = (r-y[0])-w; 134 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 135 else return n; 136 } 137 /* 138 * all other (large) arguments 139 */ 140 if(ix>=0x7ff00000) { /* x is inf or NaN */ 141 y[0]=y[1]=x-x; return 0; 142 } 143 /* set z = scalbn(|x|,ilogb(x)-23) */ 144 GET_LOW_WORD(low,x); 145 SET_LOW_WORD(z,low); 146 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 147 SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 148 for(i=0;i<2;i++) { 149 tx[i] = (double)((int32_t)(z)); 150 z = (z-tx[i])*two24; 151 } 152 tx[2] = z; 153 nx = 3; 154 while(tx[nx-1]==zero) nx--; /* skip zero term */ 155 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 156 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 157 return n; 158} 159