source: trunk/sys/libm/e_rem_pio2.c @ 9

Last change on this file since 9 was 1, checked in by alain, 7 years ago

First import

File size: 4.7 KB
Line 
1
2/* @(#)e_rem_pio2.c 5.1 93/09/24 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 *
13 */
14
15/* __ieee754_rem_pio2(x,y)
16 *
17 * return the remainder of x rem pi/2 in y[0]+y[1]
18 * use __kernel_rem_pio2()
19 */
20
21#include <libm/fdlibm.h>
22
23/*
24 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
25 */
26#ifdef __STDC__
27static const int two_over_pi[] = {
28#else
29static int two_over_pi[] = {
30#endif
310xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
320x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
330x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
340xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
350x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
360x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
370x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
380xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
390x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
400x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
410x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
42};
43
44#ifdef __STDC__
45static const int npio2_hw[] = {
46#else
47static int npio2_hw[] = {
48#endif
490x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
500x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
510x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
520x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
530x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
540x404858EB, 0x404921FB,
55};
56
57/*
58 * invpio2:  53 bits of 2/pi
59 * pio2_1:   first  33 bit of pi/2
60 * pio2_1t:  pi/2 - pio2_1
61 * pio2_2:   second 33 bit of pi/2
62 * pio2_2t:  pi/2 - (pio2_1+pio2_2)
63 * pio2_3:   third  33 bit of pi/2
64 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
65 */
66
67#ifdef __STDC__
68static const double 
69#else
70static double 
71#endif
72zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
73half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
74two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
75invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
76pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
77pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
78pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
79pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
80pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
81pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
82
83#ifdef __STDC__
84        int __ieee754_rem_pio2(double x, double *y)
85#else
86        int __ieee754_rem_pio2(x,y)
87        double x,y[];
88#endif
89{
90        double z,w,t,r,fn;
91        double tx[3];
92        int e0,i,j,nx,n,ix,hx,i0;
93
94        i0 = ((*(int*)&two24)>>30)^1;   /* high word index */
95        hx = *(i0+(int*)&x);            /* high word of x */
96        ix = hx&0x7fffffff;
97        if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
98            {y[0] = x; y[1] = 0; return 0;}
99        if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
100            t  = fabs(x);
101            n  = (int) (t*invpio2+half);
102            fn = (double)n;
103            r  = t-fn*pio2_1;
104            w  = fn*pio2_1t;    /* 1st round good to 85 bit */
105            if(n<32&&ix!=npio2_hw[n-1]) {       
106                y[0] = r-w;     /* quick check no cancellation */
107            } else {
108                j  = ix>>20;
109                y[0] = r-w; 
110                i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);
111                if(i>16) {  /* 2nd iteration needed, good to 118 */
112                    t  = r;
113                    w  = fn*pio2_2;     
114                    r  = t-w;
115                    w  = fn*pio2_2t-((t-r)-w); 
116                    y[0] = r-w;
117                    i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);
118                    if(i>49)  { /* 3rd iteration need, 151 bits acc */
119                        t  = r; /* will cover all possible cases */
120                        w  = fn*pio2_3; 
121                        r  = t-w;
122                        w  = fn*pio2_3t-((t-r)-w);     
123                        y[0] = r-w;
124                    }
125                }
126            }
127            y[1] = (r-y[0])-w;
128            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
129            else         return n;
130        }
131    /*
132     * all other (large) arguments
133     */
134        if(ix>=0x7ff00000) {            /* x is inf or NaN */
135            y[0]=y[1]=x-x; return 0;
136        }
137    /* set z = scalbn(|x|,ilogb(x)-23) */
138        *(1-i0+(int*)&z) = *(1-i0+(int*)&x);
139        e0      = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */
140        *(i0+(int*)&z) = ix - (e0<<20);
141        for(i=0;i<2;i++) {
142                tx[i] = (double)((int)(z));
143                z     = (z-tx[i])*two24;
144        }
145        tx[2] = z;
146        nx = 3;
147        while(tx[nx-1]==zero) nx--;     /* skip zero term */
148        n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
149        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
150        return n;
151}
Note: See TracBrowser for help on using the repository browser.