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