Subversion Repositories shark

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/* e_rem_pio2f.c -- float version of e_rem_pio2.c
2
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3
 */
4
 
5
/*
6
 * ====================================================
7
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8
 *
9
 * Developed at SunPro, a Sun Microsystems, Inc. business.
10
 * Permission to use, copy, modify, and distribute this
11
 * software is freely granted, provided that this notice
12
 * is preserved.
13
 * ====================================================
14
 */
15
 
16
#ifndef lint
17
static char rcsid[] = "$\Id: e_rem_pio2f.c,v 1.3 1995/05/30 05:48:38 rgrimes Exp $";
18
#endif
19
 
20
/* __ieee754_rem_pio2f(x,y)
21
 *
22
 * return the remainder of x rem pi/2 in y[0]+y[1]
23
 * use __kernel_rem_pio2f()
24
 */
25
 
26
#include "math.h"
27
#include "math_private.h"
28
 
29
/*
30
 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
31
 */
32
#ifdef __STDC__
33
static const int32_t two_over_pi[] = {
34
#else
35
static int32_t two_over_pi[] = {
36
#endif
37
0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
38
0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
39
0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
40
0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
41
0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
42
0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
43
0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
44
0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
45
0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
46
0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
47
0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
48
0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
49
0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
50
0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
51
0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
52
0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
53
0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
54
0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
55
0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
56
0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
57
0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
58
0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
59
};
60
 
61
/* This array is like the one in e_rem_pio2.c, but the numbers are
62
   single precision and the last 8 bits are forced to 0.  */
63
#ifdef __STDC__
64
static const int32_t npio2_hw[] = {
65
#else
66
static int32_t npio2_hw[] = {
67
#endif
68
0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
69
0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
70
0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
71
0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
72
0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
73
0x4242c700, 0x42490f00
74
};
75
 
76
/*
77
 * invpio2:  24 bits of 2/pi
78
 * pio2_1:   first  17 bit of pi/2
79
 * pio2_1t:  pi/2 - pio2_1
80
 * pio2_2:   second 17 bit of pi/2
81
 * pio2_2t:  pi/2 - (pio2_1+pio2_2)
82
 * pio2_3:   third  17 bit of pi/2
83
 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
84
 */
85
 
86
#ifdef __STDC__
87
static const float
88
#else
89
static float
90
#endif
91
zero =  0.0000000000e+00, /* 0x00000000 */
92
half =  5.0000000000e-01, /* 0x3f000000 */
93
two8 =  2.5600000000e+02, /* 0x43800000 */
94
invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
95
pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
96
pio2_1t =  1.0804334124e-05, /* 0x37354443 */
97
pio2_2  =  1.0804273188e-05, /* 0x37354400 */
98
pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
99
pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
100
pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
101
 
102
#ifdef __STDC__
103
        int32_t __ieee754_rem_pio2f(float x, float *y)
104
#else
105
        int32_t __ieee754_rem_pio2f(x,y)
106
        float x,y[];
107
#endif
108
{
109
        float z,w,t,r,fn;
110
        float tx[3];
111
        int32_t e0,i,j,nx,n,ix,hx;
112
 
113
        GET_FLOAT_WORD(hx,x);
114
        ix = hx&0x7fffffff;
115
        if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
116
            {y[0] = x; y[1] = 0; return 0;}
117
        if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
118
            if(hx>0) {
119
                z = x - pio2_1;
120
                if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
121
                    y[0] = z - pio2_1t;
122
                    y[1] = (z-y[0])-pio2_1t;
123
                } else {                /* near pi/2, use 24+24+24 bit pi */
124
                    z -= pio2_2;
125
                    y[0] = z - pio2_2t;
126
                    y[1] = (z-y[0])-pio2_2t;
127
                }
128
                return 1;
129
            } else {    /* negative x */
130
                z = x + pio2_1;
131
                if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
132
                    y[0] = z + pio2_1t;
133
                    y[1] = (z-y[0])+pio2_1t;
134
                } else {                /* near pi/2, use 24+24+24 bit pi */
135
                    z += pio2_2;
136
                    y[0] = z + pio2_2t;
137
                    y[1] = (z-y[0])+pio2_2t;
138
                }
139
                return -1;
140
            }
141
        }
142
        if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
143
            t  = fabsf(x);
144
            n  = (int32_t) (t*invpio2+half);
145
            fn = (float)n;
146
            r  = t-fn*pio2_1;
147
            w  = fn*pio2_1t;    /* 1st round good to 40 bit */
148
            if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) {
149
                y[0] = r-w;     /* quick check no cancellation */
150
            } else {
151
                u_int32_t high;
152
                j  = ix>>23;
153
                y[0] = r-w;
154
                GET_FLOAT_WORD(high,y[0]);
155
                i = j-((high>>23)&0xff);
156
                if(i>8) {  /* 2nd iteration needed, good to 57 */
157
                    t  = r;
158
                    w  = fn*pio2_2;
159
                    r  = t-w;
160
                    w  = fn*pio2_2t-((t-r)-w);
161
                    y[0] = r-w;
162
                    GET_FLOAT_WORD(high,y[0]);
163
                    i = j-((high>>23)&0xff);
164
                    if(i>25)  { /* 3rd iteration need, 74 bits acc */
165
                        t  = r; /* will cover all possible cases */
166
                        w  = fn*pio2_3;
167
                        r  = t-w;
168
                        w  = fn*pio2_3t-((t-r)-w);
169
                        y[0] = r-w;
170
                    }
171
                }
172
            }
173
            y[1] = (r-y[0])-w;
174
            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
175
            else         return n;
176
        }
177
    /*
178
     * all other (large) arguments
179
     */
180
        if(ix>=0x7f800000) {            /* x is inf or NaN */
181
            y[0]=y[1]=x-x; return 0;
182
        }
183
    /* set z = scalbn(|x|,ilogb(x)-7) */
184
        e0      = (ix>>23)-134;         /* e0 = ilogb(z)-7; */
185
        SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
186
        for(i=0;i<2;i++) {
187
                tx[i] = (float)((int32_t)(z));
188
                z     = (z-tx[i])*two8;
189
        }
190
        tx[2] = z;
191
        nx = 3;
192
        while(tx[nx-1]==zero) nx--;     /* skip zero term */
193
        n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
194
        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
195
        return n;
196
}