Deleted Added
sdiff udiff text old ( 152535 ) new ( 152596 )
full compact
e_rem_pio2f.c (152535) e_rem_pio2f.c (152596)
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 * Debugged and optimized by Bruce D. Evans.
4 */
5
6/*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 *
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
15 */
16
17#ifndef lint
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 * Debugged and optimized by Bruce D. Evans.
4 */
5
6/*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 *
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
15 */
16
17#ifndef lint
18static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 152535 2005-11-17 02:20:04Z bde $";
18static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 152596 2005-11-19 02:38:27Z bde $";
19#endif
20
21/* __ieee754_rem_pio2f(x,y)
22 *
23 * return the remainder of x rem pi/2 in y[0]+y[1]
24 * use double precision internally
25 * use __kernel_rem_pio2() for large x
26 */

--- 15 unchanged lines hidden (view full) ---

420xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
430x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
440x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
450x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
46};
47
48/*
49 * invpio2: 53 bits of 2/pi
19#endif
20
21/* __ieee754_rem_pio2f(x,y)
22 *
23 * return the remainder of x rem pi/2 in y[0]+y[1]
24 * use double precision internally
25 * use __kernel_rem_pio2() for large x
26 */

--- 15 unchanged lines hidden (view full) ---

420xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
430x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
440x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
450x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
46};
47
48/*
49 * invpio2: 53 bits of 2/pi
50 * e1pio2: 1*pi/2 rounded to 53 bits
51 * e2pio2: 2*pi/2 rounded to 53 bits
52 * e3pio2: 3*pi/2 rounded to 53 bits
53 * e4pio2: 4*pi/2 rounded to 53 bits
54 * pio2_1: first 33 bit of pi/2
55 * pio2_1t: pi/2 - pio2_1
56 */
57
58static const double
59zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
60half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
61two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
62invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
50 * pio2_1: first 33 bit of pi/2
51 * pio2_1t: pi/2 - pio2_1
52 */
53
54static const double
55zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
56half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
57two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
58invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
63e1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
64e2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
65e3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
66e4pio2 = 4*M_PI_2, /* 0x401921FB, 0x54442D18 */
67pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
68pio2_1t = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
69
70 int32_t __ieee754_rem_pio2f(float x, float *y)
71{
72 double z,w,t,r,fn;
73 double tx[3];
74 int32_t e0,i,nx,n,ix,hx;
75
76 GET_FLOAT_WORD(hx,x);
77 ix = hx&0x7fffffff;
59pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
60pio2_1t = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
61
62 int32_t __ieee754_rem_pio2f(float x, float *y)
63{
64 double z,w,t,r,fn;
65 double tx[3];
66 int32_t e0,i,nx,n,ix,hx;
67
68 GET_FLOAT_WORD(hx,x);
69 ix = hx&0x7fffffff;
78 if(ix<=0x3f490fda) /* |x| ~<= pi/4, reduction is null */
79 {y[0] = x; y[1] = 0; return 0;}
80 /* 53 bit pi is good enough for special cases */
81 if(ix<=0x407b53d1) { /* |x| ~<= 5*pi/4 */
82 if(ix<=0x4016cbe3) { /* |x| ~<= 3*pi/4 */
83 if(hx>0) {
84 z = x - e1pio2;
85 n = 1;
86 } else {
87 z = x + e1pio2;
88 n = 3;
89 }
90 y[0] = z;
91 y[1] = z - y[0];
92 return n;
93 } else {
94 if(hx>0)
95 z = x - e2pio2;
96 else
97 z = x + e2pio2;
98 y[0] = z;
99 y[1] = z - y[0];
100 return 2;
101 }
102 }
103 if(ix<=0x40e231d5) { /* |x| ~<= 9*pi/4*/
104 if(ix<=0x40afeddf) { /* |x| ~<= 7*pi/4 */
105 if(hx>0) {
106 z = x - e3pio2;
107 n = 3;
108 } else {
109 z = x + e3pio2;
110 n = 1;
111 }
112 y[0] = z;
113 y[1] = z - y[0];
114 return n;
115 } else {
116 if(hx>0)
117 z = x - e4pio2;
118 else
119 z = x + e4pio2;
120 y[0] = z;
121 y[1] = z - y[0];
122 return 0;
123 }
124 }
125 /* 33+53 bit pi is good enough for medium size */
126 if(ix<=0x49490f80) { /* |x| ~<= 2^19*(pi/2), medium size */
127 t = fabsf(x);
128 n = (int32_t) (t*invpio2+half);
129 fn = (double)n;
130 r = t-fn*pio2_1;
131 w = fn*pio2_1t;
132 y[0] = r-w;

--- 29 unchanged lines hidden --- 		70 /* 33+53 bit pi is good enough for medium size */
71 if(ix<=0x49490f80) { /* |x| ~<= 2^19*(pi/2), medium size */
72 t = fabsf(x);
73 n = (int32_t) (t*invpio2+half);
74 fn = (double)n;
75 r = t-fn*pio2_1;
76 w = fn*pio2_1t;
77 y[0] = r-w;

--- 29 unchanged lines hidden ---