← Back to branch summary

~ubuntu-branches/ubuntu/precise/eglibc/precise

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/dbl-64/e_rem_pio2.c

  • Committer: Package Import Robot
  • Author(s): Matthias Klose
  • Date: 2011-10-04 17:48:26 UTC
  • mfrom: (216.1.23 oneiric)
  • Revision ID: package-import@ubuntu.com-20111004174826-2cyb9ewn3ucymlsx
Tags: 2.13-20ubuntu5
https://launchpad.net/bugs/853688
libc6-dev: Don't break the current {gnat,gcj}-4.4-base versons. LP: #853688.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/dbl-64/e_rem_pio2.c
1
 
#ifdef NOT_NEEDED_ANYMORE
2
 
 
 
1
/* @(#)e_rem_pio2.c 5.1 93/09/24 */
3
2
/*
4
3
 * ====================================================
5
4
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6
5
 *
7
6
 * Developed at SunPro, a Sun Microsystems, Inc. business.
8
7
 * Permission to use, copy, modify, and distribute this
9
 
 * software is freely granted, provided that this notice
 
8
 * software is freely granted, provided that this notice 
10
9
 * is preserved.
11
10
 * ====================================================
12
11
 */
13
12
 
 
13
#if defined(LIBM_SCCS) && !defined(lint)
 
14
static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
 
15
#endif
 
16
 
14
17
/* __ieee754_rem_pio2(x,y)
15
 
 *
16
 
 * return the remainder of x rem pi/2 in y[0]+y[1]
 
18
 * 
 
19
 * return the remainder of x rem pi/2 in y[0]+y[1] 
17
20
 * use __kernel_rem_pio2()
18
21
 */
19
22
 
21
24
#include "math_private.h"
22
25
 
23
26
/*
24
 
 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
 
27
 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
25
28
 */
 
29
#ifdef __STDC__
26
30
static const int32_t two_over_pi[] = {
27
 
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
28
 
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
29
 
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
30
 
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
31
 
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
32
 
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
33
 
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
34
 
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
35
 
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
36
 
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
37
 
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
 
31
#else
 
32
static int32_t two_over_pi[] = {
 
33
#endif
 
34
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
 
35
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
 
36
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
 
37
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
 
38
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
 
39
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
 
40
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
 
41
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
 
42
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
 
43
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
 
44
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
38
45
};
39
46
 
 
47
#ifdef __STDC__
40
48
static const int32_t npio2_hw[] = {
 
49
#else
 
50
static int32_t npio2_hw[] = {
 
51
#endif
41
52
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
42
53
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
43
54
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
56
67
 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
57
68
 */
58
69
 
59
 
static const double
 
70
#ifdef __STDC__
 
71
static const double 
 
72
#else
 
73
static double 
 
74
#endif
60
75
zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
61
76
half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
62
77
two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
68
83
pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
69
84
pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
70
85
 
71
 
int32_t
72
 
__ieee754_rem_pio2(double x, double *y)
 
86
#ifdef __STDC__
 
87
        int32_t __ieee754_rem_pio2(double x, double *y)
 
88
#else
 
89
        int32_t __ieee754_rem_pio2(x,y)
 
90
        double x,y[];
 
91
#endif
73
92
{
74
93
        double z,w,t,r,fn;
75
94
        double tx[3];
81
100
        if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
82
101
            {y[0] = x; y[1] = 0; return 0;}
83
102
        if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
84
 
            if(hx>0) {
 
103
            if(hx>0) { 
85
104
                z = x - pio2_1;
86
 
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
 
105
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
87
106
                    y[0] = z - pio2_1t;
88
107
                    y[1] = (z-y[0])-pio2_1t;
89
108
                } else {                /* near pi/2, use 33+33+53 bit pi */
94
113
                return 1;
95
114
            } else {    /* negative x */
96
115
                z = x + pio2_1;
97
 
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
 
116
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
98
117
                    y[0] = z + pio2_1t;
99
118
                    y[1] = (z-y[0])+pio2_1t;
100
119
                } else {                /* near pi/2, use 33+33+53 bit pi */
111
130
            fn = (double)n;
112
131
            r  = t-fn*pio2_1;
113
132
            w  = fn*pio2_1t;    /* 1st round good to 85 bit */
114
 
            if(n<32&&ix!=npio2_hw[n-1]) {
 
133
            if(n<32&&ix!=npio2_hw[n-1]) {       
115
134
                y[0] = r-w;     /* quick check no cancellation */
116
135
            } else {
117
 
                u_int32_t high;
118
 
                j  = ix>>20;
119
 
                y[0] = r-w;
 
136
                u_int32_t high;
 
137
                j  = ix>>20;
 
138
                y[0] = r-w; 
120
139
                GET_HIGH_WORD(high,y[0]);
121
 
                i = j-((high>>20)&0x7ff);
122
 
                if(i>16) {  /* 2nd iteration needed, good to 118 */
 
140
                i = j-((high>>20)&0x7ff);
 
141
                if(i>16) {  /* 2nd iteration needed, good to 118 */
123
142
                    t  = r;
124
 
                    w  = fn*pio2_2;
 
143
                    w  = fn*pio2_2;     
125
144
                    r  = t-w;
126
 
                    w  = fn*pio2_2t-((t-r)-w);
 
145
                    w  = fn*pio2_2t-((t-r)-w);  
127
146
                    y[0] = r-w;
128
147
                    GET_HIGH_WORD(high,y[0]);
129
148
                    i = j-((high>>20)&0x7ff);
130
149
                    if(i>49)  { /* 3rd iteration need, 151 bits acc */
131
 
                        t  = r; /* will cover all possible cases */
132
 
                        w  = fn*pio2_3;
133
 
                        r  = t-w;
134
 
                        w  = fn*pio2_3t-((t-r)-w);
135
 
                        y[0] = r-w;
 
150
                        t  = r; /* will cover all possible cases */
 
151
                        w  = fn*pio2_3; 
 
152
                        r  = t-w;
 
153
                        w  = fn*pio2_3t-((t-r)-w);      
 
154
                        y[0] = r-w;
136
155
                    }
137
156
                }
138
157
            }
139
158
            y[1] = (r-y[0])-w;
140
 
            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 
159
            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
141
160
            else         return n;
142
161
        }
143
 
    /*
 
162
    /* 
144
163
     * all other (large) arguments
145
164
     */
146
165
        if(ix>=0x7ff00000) {            /* x is inf or NaN */
149
168
    /* set z = scalbn(|x|,ilogb(x)-23) */
150
169
        GET_LOW_WORD(low,x);
151
170
        SET_LOW_WORD(z,low);
152
 
        e0      = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */
 
171
        e0      = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */
153
172
        SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
154
173
        for(i=0;i<2;i++) {
155
174
                tx[i] = (double)((int32_t)(z));
162
181
        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
163
182
        return n;
164
183
}
165
 
 
166
 
#endif