Branch data Line data Source code
1 : : /* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
2 : : /*
3 : : * ====================================================
4 : : * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 : : * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
6 : : *
7 : : * Developed at SunSoft, 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 : : * Optimized by Bruce D. Evans.
14 : : */
15 : :
16 : : #include "cdefs-compat.h"
17 : : //__FBSDID("$FreeBSD: src/lib/msun/ld80/e_rem_pio2l.h,v 1.3 2011/06/18 13:56:33 benl Exp $");
18 : :
19 : : /* ld80 version of __ieee754_rem_pio2l(x,y)
20 : : *
21 : : * return the remainder of x rem pi/2 in y[0]+y[1]
22 : : * use __kernel_rem_pio2()
23 : : */
24 : :
25 : : #include <float.h>
26 : : #include <openlibm_math.h>
27 : :
28 : : #include "math_private.h"
29 : :
30 : : #define BIAS (LDBL_MAX_EXP - 1)
31 : :
32 : : /*
33 : : * invpio2: 64 bits of 2/pi
34 : : * pio2_1: first 39 bits of pi/2
35 : : * pio2_1t: pi/2 - pio2_1
36 : : * pio2_2: second 39 bits of pi/2
37 : : * pio2_2t: pi/2 - (pio2_1+pio2_2)
38 : : * pio2_3: third 39 bits of pi/2
39 : : * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
40 : : */
41 : :
42 : : static const double
43 : : zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
44 : : two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
45 : : pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
46 : : pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
47 : : pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */
48 : :
49 : : #if defined(__amd64__) || defined(__i386__)
50 : : /* Long double constants are slow on these arches, and broken on i386. */
51 : : static const volatile double
52 : : invpio2hi = 6.3661977236758138e-01, /* 0x145f306dc9c883.0p-53 */
53 : : invpio2lo = -3.9356538861223811e-17, /* -0x16b00000000000.0p-107 */
54 : : pio2_1thi = -1.0746346554971943e-12, /* -0x12e7b9676733af.0p-92 */
55 : : pio2_1tlo = 8.8451028997905949e-29, /* 0x1c080000000000.0p-146 */
56 : : pio2_2thi = 6.3683171635109499e-25, /* 0x18a2e03707344a.0p-133 */
57 : : pio2_2tlo = 2.3183081793789774e-41, /* 0x10280000000000.0p-187 */
58 : : pio2_3thi = -2.7529965190440717e-37, /* -0x176b7ed8fbbacc.0p-174 */
59 : : pio2_3tlo = -4.2006647512740502e-54; /* -0x19c00000000000.0p-230 */
60 : : #define invpio2 ((long double)invpio2hi + invpio2lo)
61 : : #define pio2_1t ((long double)pio2_1thi + pio2_1tlo)
62 : : #define pio2_2t ((long double)pio2_2thi + pio2_2tlo)
63 : : #define pio2_3t ((long double)pio2_3thi + pio2_3tlo)
64 : : #else
65 : : static const long double
66 : : invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */
67 : : pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
68 : : pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */
69 : : pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
70 : : #endif
71 : :
72 : : //VBS
73 : : //static inline __always_inline int
74 : : //__ieee754_rem_pio2l(long double x, long double *y)
75 : :
76 : : static inline int
77 : 0 : __ieee754_rem_pio2l(long double x, long double *y)
78 : : {
79 : : union IEEEl2bits u,u1;
80 : : long double z,w,t,r,fn;
81 : : double tx[3],ty[2];
82 : : int e0,ex,i,j,nx,n;
83 : : int16_t expsign;
84 : :
85 : 0 : u.e = x;
86 : 0 : expsign = u.xbits.expsign;
87 : 0 : ex = expsign & 0x7fff;
88 [ # # # # : 0 : if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
# # ]
89 : : /* |x| ~< 2^25*(pi/2), medium size */
90 : : /* Use a specialized rint() to get fn. Assume round-to-nearest. */
91 : 0 : fn = x*invpio2+0x1.8p63;
92 : 0 : fn = fn-0x1.8p63;
93 : : #ifdef HAVE_EFFICIENT_IRINT
94 : : n = irint(fn);
95 : : #else
96 : 0 : n = fn;
97 : : #endif
98 : 0 : r = x-fn*pio2_1;
99 : 0 : w = fn*pio2_1t; /* 1st round good to 102 bit */
100 : : {
101 : : union IEEEl2bits u2;
102 : : int ex1;
103 : 0 : j = ex;
104 : 0 : y[0] = r-w;
105 : 0 : u2.e = y[0];
106 : 0 : ex1 = u2.xbits.expsign & 0x7fff;
107 : 0 : i = j-ex1;
108 [ # # ]: 0 : if(i>22) { /* 2nd iteration needed, good to 141 */
109 : 0 : t = r;
110 : 0 : w = fn*pio2_2;
111 : 0 : r = t-w;
112 : 0 : w = fn*pio2_2t-((t-r)-w);
113 : 0 : y[0] = r-w;
114 : 0 : u2.e = y[0];
115 : 0 : ex1 = u2.xbits.expsign & 0x7fff;
116 : 0 : i = j-ex1;
117 [ # # ]: 0 : if(i>61) { /* 3rd iteration need, 180 bits acc */
118 : 0 : t = r; /* will cover all possible cases */
119 : 0 : w = fn*pio2_3;
120 : 0 : r = t-w;
121 : 0 : w = fn*pio2_3t-((t-r)-w);
122 : 0 : y[0] = r-w;
123 : : }
124 : : }
125 : : }
126 : 0 : y[1] = (r-y[0])-w;
127 : 0 : return n;
128 : : }
129 : : /*
130 : : * all other (large) arguments
131 : : */
132 [ # # ]: 0 : if(ex==0x7fff) { /* x is inf or NaN */
133 : 0 : y[0]=y[1]=x-x; return 0;
134 : : }
135 : : /* set z = scalbn(|x|,ilogb(x)-23) */
136 : 0 : u1.e = x;
137 : 0 : e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
138 : 0 : u1.xbits.expsign = ex - e0;
139 : 0 : z = u1.e;
140 [ # # ]: 0 : for(i=0;i<2;i++) {
141 : 0 : tx[i] = (double)((int32_t)(z));
142 : 0 : z = (z-tx[i])*two24;
143 : : }
144 : 0 : tx[2] = z;
145 : 0 : nx = 3;
146 [ # # ]: 0 : while(tx[nx-1]==zero) nx--; /* skip zero term */
147 : 0 : n = __kernel_rem_pio2(tx,ty,e0,nx,2);
148 : 0 : r = (long double)ty[0] + ty[1];
149 : 0 : w = ty[1] - (r - ty[0]);
150 [ # # ]: 0 : if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
151 : 0 : y[0] = r; y[1] = w; return n;
152 : : }
|
