]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * Written by International Business Machines Corp. |
688903eb | 4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
e4d82761 UD |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify | |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
50944bca | 10 | * |
e4d82761 UD |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
e4d82761 UD |
18 | */ |
19 | /*******************************************************************/ | |
20 | /* */ | |
50944bca | 21 | /* MODULE_NAME: branred.c */ |
e4d82761 UD |
22 | /* */ |
23 | /* FUNCTIONS: branred */ | |
50944bca | 24 | /* */ |
e4d82761 UD |
25 | /* FILES NEEDED: branred.h mydefs.h endian.h mpa.h */ |
26 | /* mha.c */ | |
27 | /* */ | |
28 | /* Routine branred() performs range reduction of a double number */ | |
29 | /* x into Double length number a+aa,such that */ | |
30 | /* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ | |
31 | /* Routine returns the integer (n mod 4) of the above description */ | |
32 | /* of x. */ | |
33 | /*******************************************************************/ | |
34 | ||
35 | #include "endian.h" | |
36 | #include "mydefs.h" | |
37 | #include "branred.h" | |
0e9be4db | 38 | #include <math.h> |
1ed0291c | 39 | #include <math_private.h> |
e4d82761 | 40 | |
31d3cc00 UD |
41 | #ifndef SECTION |
42 | # define SECTION | |
43 | #endif | |
44 | ||
e4d82761 UD |
45 | |
46 | /*******************************************************************/ | |
47 | /* Routine branred() performs range reduction of a double number */ | |
48 | /* x into Double length number a+aa,such that */ | |
49 | /* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ | |
50 | /* Routine return integer (n mod 4) */ | |
51 | /*******************************************************************/ | |
31d3cc00 UD |
52 | int |
53 | SECTION | |
54 | __branred(double x, double *a, double *aa) | |
e4d82761 | 55 | { |
50944bca | 56 | int i,k; |
50944bca | 57 | mynumber u,gor; |
50944bca UD |
58 | double r[6],s,t,sum,b,bb,sum1,sum2,b1,bb1,b2,bb2,x1,x2,t1,t2; |
59 | ||
e4d82761 UD |
60 | x*=tm600.x; |
61 | t=x*split; /* split x to two numbers */ | |
62 | x1=t-(t-x); | |
63 | x2=x-x1; | |
64 | sum=0; | |
65 | u.x = x1; | |
50944bca | 66 | k = (u.i[HIGH_HALF]>>20)&2047; |
e4d82761 UD |
67 | k = (k-450)/24; |
68 | if (k<0) | |
69 | k=0; | |
70 | gor.x = t576.x; | |
71 | gor.i[HIGH_HALF] -= ((k*24)<<20); | |
50944bca | 72 | for (i=0;i<6;i++) |
e4d82761 UD |
73 | { r[i] = x1*toverp[k+i]*gor.x; gor.x *= tm24.x; } |
74 | for (i=0;i<3;i++) { | |
50944bca | 75 | s=(r[i]+big.x)-big.x; |
e4d82761 UD |
76 | sum+=s; |
77 | r[i]-=s; | |
78 | } | |
79 | t=0; | |
50944bca | 80 | for (i=0;i<6;i++) |
e4d82761 | 81 | t+=r[5-i]; |
50944bca | 82 | bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5]; |
e4d82761 UD |
83 | s=(t+big.x)-big.x; |
84 | sum+=s; | |
85 | t-=s; | |
86 | b=t+bb; | |
87 | bb=(t-b)+bb; | |
88 | s=(sum+big1.x)-big1.x; | |
89 | sum-=s; | |
90 | b1=b; | |
91 | bb1=bb; | |
92 | sum1=sum; | |
93 | sum=0; | |
50944bca | 94 | |
e4d82761 UD |
95 | u.x = x2; |
96 | k = (u.i[HIGH_HALF]>>20)&2047; | |
97 | k = (k-450)/24; | |
98 | if (k<0) | |
99 | k=0; | |
100 | gor.x = t576.x; | |
101 | gor.i[HIGH_HALF] -= ((k*24)<<20); | |
102 | for (i=0;i<6;i++) | |
103 | { r[i] = x2*toverp[k+i]*gor.x; gor.x *= tm24.x; } | |
104 | for (i=0;i<3;i++) { | |
50944bca | 105 | s=(r[i]+big.x)-big.x; |
e4d82761 UD |
106 | sum+=s; |
107 | r[i]-=s; | |
108 | } | |
109 | t=0; | |
110 | for (i=0;i<6;i++) | |
111 | t+=r[5-i]; | |
50944bca | 112 | bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5]; |
e4d82761 UD |
113 | s=(t+big.x)-big.x; |
114 | sum+=s; | |
115 | t-=s; | |
116 | b=t+bb; | |
117 | bb=(t-b)+bb; | |
118 | s=(sum+big1.x)-big1.x; | |
119 | sum-=s; | |
50944bca | 120 | |
e4d82761 UD |
121 | b2=b; |
122 | bb2=bb; | |
123 | sum2=sum; | |
50944bca | 124 | |
e4d82761 UD |
125 | sum=sum1+sum2; |
126 | b=b1+b2; | |
0e9be4db | 127 | bb = (fabs(b1)>fabs(b2))? (b1-b)+b2 : (b2-b)+b1; |
50944bca | 128 | if (b > 0.5) |
e4d82761 | 129 | {b-=1.0; sum+=1.0;} |
50944bca | 130 | else if (b < -0.5) |
e4d82761 UD |
131 | {b+=1.0; sum-=1.0;} |
132 | s=b+(bb+bb1+bb2); | |
133 | t=((b-s)+bb)+(bb1+bb2); | |
134 | b=s*split; | |
135 | t1=b-(b-s); | |
136 | t2=s-t1; | |
137 | b=s*hp0.x; | |
138 | bb=(((t1*mp1.x-b)+t1*mp2.x)+t2*mp1.x)+(t2*mp2.x+s*hp1.x+t*hp0.x); | |
139 | s=b+bb; | |
140 | t=(b-s)+bb; | |
141 | *a=s; | |
142 | *aa=t; | |
143 | return ((int) sum)&3; /* return quater of unit circle */ | |
144 | } |