2 * IBM Accurate Mathematical Library
3 * Copyright (C) 2001-2023 Free Software Foundation, Inc.
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.
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.
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/>.
18 /*******************************************************************/
20 /* MODULE_NAME: branred.c */
22 /* FUNCTIONS: branred */
24 /* FILES NEEDED: branred.h mydefs.h endian.h mpa.h */
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 */
32 /*******************************************************************/
38 #include <math_private.h>
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 /*******************************************************************/
53 __branred(double x
, double *a
, double *aa
)
57 double r
[6],s
,t
,sum
,b
,bb
,sum1
,sum2
,b1
,bb1
,b2
,bb2
,x1
,x2
,t1
,t2
;
60 t
=x
*split
; /* split x to two numbers */
65 k
= (u
.i
[HIGH_HALF
]>>20)&2047;
70 gor
.i
[HIGH_HALF
] -= ((k
*24)<<20);
72 { r
[i
] = x1
*toverp
[k
+i
]*gor
.x
; gor
.x
*= tm24
.x
; }
81 bb
=(((((r
[0]-t
)+r
[1])+r
[2])+r
[3])+r
[4])+r
[5];
87 s
=(sum
+big1
.x
)-big1
.x
;
95 k
= (u
.i
[HIGH_HALF
]>>20)&2047;
100 gor
.i
[HIGH_HALF
] -= ((k
*24)<<20);
102 { r
[i
] = x2
*toverp
[k
+i
]*gor
.x
; gor
.x
*= tm24
.x
; }
104 s
=(r
[i
]+big
.x
)-big
.x
;
111 bb
=(((((r
[0]-t
)+r
[1])+r
[2])+r
[3])+r
[4])+r
[5];
117 s
=(sum
+big1
.x
)-big1
.x
;
126 bb
= (fabs(b1
)>fabs(b2
))? (b1
-b
)+b2
: (b2
-b
)+b1
;
132 t
=((b
-s
)+bb
)+(bb1
+bb2
);
137 bb
=(((t1
*mp1
.x
-b
)+t1
*mp2
.x
)+t2
*mp1
.x
)+(t2
*mp2
.x
+s
*hp1
.x
+t
*hp0
.x
);
142 return ((int) sum
)&3; /* return quarter of unit circle */