105a0b428SJohn Marino /* e_rem_pio2f.c -- float version of e_rem_pio2.c
205a0b428SJohn Marino * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
305a0b428SJohn Marino */
405a0b428SJohn Marino
505a0b428SJohn Marino /*
605a0b428SJohn Marino * ====================================================
705a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
805a0b428SJohn Marino *
905a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
1005a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
1105a0b428SJohn Marino * software is freely granted, provided that this notice
1205a0b428SJohn Marino * is preserved.
1305a0b428SJohn Marino * ====================================================
1405a0b428SJohn Marino */
1505a0b428SJohn Marino
1605a0b428SJohn Marino /* __ieee754_rem_pio2f(x,y)
1705a0b428SJohn Marino *
1805a0b428SJohn Marino * return the remainder of x rem pi/2 in y[0]+y[1]
1905a0b428SJohn Marino * use __kernel_rem_pio2f()
2005a0b428SJohn Marino */
2105a0b428SJohn Marino
2205a0b428SJohn Marino #include "math.h"
2305a0b428SJohn Marino #include "math_private.h"
2405a0b428SJohn Marino
2505a0b428SJohn Marino /*
2605a0b428SJohn Marino * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
2705a0b428SJohn Marino */
2805a0b428SJohn Marino static const int32_t two_over_pi[] = {
2905a0b428SJohn Marino 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
3005a0b428SJohn Marino 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
3105a0b428SJohn Marino 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
3205a0b428SJohn Marino 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
3305a0b428SJohn Marino 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
3405a0b428SJohn Marino 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
3505a0b428SJohn Marino 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
3605a0b428SJohn Marino 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
3705a0b428SJohn Marino 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
3805a0b428SJohn Marino 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
3905a0b428SJohn Marino 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
4005a0b428SJohn Marino 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
4105a0b428SJohn Marino 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
4205a0b428SJohn Marino 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
4305a0b428SJohn Marino 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
4405a0b428SJohn Marino 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
4505a0b428SJohn Marino 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
4605a0b428SJohn Marino 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
4705a0b428SJohn Marino 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
4805a0b428SJohn Marino 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
4905a0b428SJohn Marino 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
5005a0b428SJohn Marino 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
5105a0b428SJohn Marino };
5205a0b428SJohn Marino
5305a0b428SJohn Marino /* This array is like the one in e_rem_pio2.c, but the numbers are
5405a0b428SJohn Marino single precision and the last 8 bits are forced to 0. */
5505a0b428SJohn Marino static const int32_t npio2_hw[] = {
5605a0b428SJohn Marino 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
5705a0b428SJohn Marino 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
5805a0b428SJohn Marino 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
5905a0b428SJohn Marino 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
6005a0b428SJohn Marino 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
6105a0b428SJohn Marino 0x4242c700, 0x42490f00
6205a0b428SJohn Marino };
6305a0b428SJohn Marino
6405a0b428SJohn Marino /*
6505a0b428SJohn Marino * invpio2: 24 bits of 2/pi
6605a0b428SJohn Marino * pio2_1: first 17 bit of pi/2
6705a0b428SJohn Marino * pio2_1t: pi/2 - pio2_1
6805a0b428SJohn Marino * pio2_2: second 17 bit of pi/2
6905a0b428SJohn Marino * pio2_2t: pi/2 - (pio2_1+pio2_2)
7005a0b428SJohn Marino * pio2_3: third 17 bit of pi/2
7105a0b428SJohn Marino * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
7205a0b428SJohn Marino */
7305a0b428SJohn Marino
7405a0b428SJohn Marino static const float
7505a0b428SJohn Marino zero = 0.0000000000e+00, /* 0x00000000 */
7605a0b428SJohn Marino half = 5.0000000000e-01, /* 0x3f000000 */
7705a0b428SJohn Marino two8 = 2.5600000000e+02, /* 0x43800000 */
7805a0b428SJohn Marino invpio2 = 6.3661980629e-01, /* 0x3f22f984 */
7905a0b428SJohn Marino pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */
8005a0b428SJohn Marino pio2_1t = 1.0804334124e-05, /* 0x37354443 */
8105a0b428SJohn Marino pio2_2 = 1.0804273188e-05, /* 0x37354400 */
8205a0b428SJohn Marino pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */
8305a0b428SJohn Marino pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */
8405a0b428SJohn Marino pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
8505a0b428SJohn Marino
8605a0b428SJohn Marino int32_t
__ieee754_rem_pio2f(float x,float * y)8705a0b428SJohn Marino __ieee754_rem_pio2f(float x, float *y)
8805a0b428SJohn Marino {
8905a0b428SJohn Marino float z,w,t,r,fn;
9005a0b428SJohn Marino float tx[3];
9105a0b428SJohn Marino int32_t e0,i,j,nx,n,ix,hx;
9205a0b428SJohn Marino
9305a0b428SJohn Marino GET_FLOAT_WORD(hx,x);
9405a0b428SJohn Marino ix = hx&0x7fffffff;
9505a0b428SJohn Marino if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
9605a0b428SJohn Marino {y[0] = x; y[1] = 0; return 0;}
9705a0b428SJohn Marino if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
9805a0b428SJohn Marino if(hx>0) {
9905a0b428SJohn Marino z = x - pio2_1;
10005a0b428SJohn Marino if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
10105a0b428SJohn Marino y[0] = z - pio2_1t;
10205a0b428SJohn Marino y[1] = (z-y[0])-pio2_1t;
10305a0b428SJohn Marino } else { /* near pi/2, use 24+24+24 bit pi */
10405a0b428SJohn Marino z -= pio2_2;
10505a0b428SJohn Marino y[0] = z - pio2_2t;
10605a0b428SJohn Marino y[1] = (z-y[0])-pio2_2t;
10705a0b428SJohn Marino }
10805a0b428SJohn Marino return 1;
10905a0b428SJohn Marino } else { /* negative x */
11005a0b428SJohn Marino z = x + pio2_1;
11105a0b428SJohn Marino if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
11205a0b428SJohn Marino y[0] = z + pio2_1t;
11305a0b428SJohn Marino y[1] = (z-y[0])+pio2_1t;
11405a0b428SJohn Marino } else { /* near pi/2, use 24+24+24 bit pi */
11505a0b428SJohn Marino z += pio2_2;
11605a0b428SJohn Marino y[0] = z + pio2_2t;
11705a0b428SJohn Marino y[1] = (z-y[0])+pio2_2t;
11805a0b428SJohn Marino }
11905a0b428SJohn Marino return -1;
12005a0b428SJohn Marino }
12105a0b428SJohn Marino }
12205a0b428SJohn Marino if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
12305a0b428SJohn Marino t = fabsf(x);
12405a0b428SJohn Marino n = (int32_t) (t*invpio2+half);
12505a0b428SJohn Marino fn = (float)n;
12605a0b428SJohn Marino r = t-fn*pio2_1;
12705a0b428SJohn Marino w = fn*pio2_1t; /* 1st round good to 40 bit */
128*74b7c7a8SJohn Marino if(n<32&&(ix&0xffffff00)!=(u_int32_t)npio2_hw[n-1]) {
12905a0b428SJohn Marino y[0] = r-w; /* quick check no cancellation */
13005a0b428SJohn Marino } else {
13105a0b428SJohn Marino u_int32_t high;
13205a0b428SJohn Marino j = ix>>23;
13305a0b428SJohn Marino y[0] = r-w;
13405a0b428SJohn Marino GET_FLOAT_WORD(high,y[0]);
13505a0b428SJohn Marino i = j-((high>>23)&0xff);
13605a0b428SJohn Marino if(i>8) { /* 2nd iteration needed, good to 57 */
13705a0b428SJohn Marino t = r;
13805a0b428SJohn Marino w = fn*pio2_2;
13905a0b428SJohn Marino r = t-w;
14005a0b428SJohn Marino w = fn*pio2_2t-((t-r)-w);
14105a0b428SJohn Marino y[0] = r-w;
14205a0b428SJohn Marino GET_FLOAT_WORD(high,y[0]);
14305a0b428SJohn Marino i = j-((high>>23)&0xff);
14405a0b428SJohn Marino if(i>25) { /* 3rd iteration need, 74 bits acc */
14505a0b428SJohn Marino t = r; /* will cover all possible cases */
14605a0b428SJohn Marino w = fn*pio2_3;
14705a0b428SJohn Marino r = t-w;
14805a0b428SJohn Marino w = fn*pio2_3t-((t-r)-w);
14905a0b428SJohn Marino y[0] = r-w;
15005a0b428SJohn Marino }
15105a0b428SJohn Marino }
15205a0b428SJohn Marino }
15305a0b428SJohn Marino y[1] = (r-y[0])-w;
15405a0b428SJohn Marino if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
15505a0b428SJohn Marino else return n;
15605a0b428SJohn Marino }
15705a0b428SJohn Marino /*
15805a0b428SJohn Marino * all other (large) arguments
15905a0b428SJohn Marino */
16005a0b428SJohn Marino if(ix>=0x7f800000) { /* x is inf or NaN */
16105a0b428SJohn Marino y[0]=y[1]=x-x; return 0;
16205a0b428SJohn Marino }
16305a0b428SJohn Marino /* set z = scalbn(|x|,ilogb(x)-7) */
16405a0b428SJohn Marino e0 = (ix>>23)-134; /* e0 = ilogb(z)-7; */
16505a0b428SJohn Marino SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
16605a0b428SJohn Marino for(i=0;i<2;i++) {
16705a0b428SJohn Marino tx[i] = (float)((int32_t)(z));
16805a0b428SJohn Marino z = (z-tx[i])*two8;
16905a0b428SJohn Marino }
17005a0b428SJohn Marino tx[2] = z;
17105a0b428SJohn Marino nx = 3;
17205a0b428SJohn Marino while(tx[nx-1]==zero) nx--; /* skip zero term */
17305a0b428SJohn Marino n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
17405a0b428SJohn Marino if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
17505a0b428SJohn Marino return n;
17605a0b428SJohn Marino }
177