1/*
2 * IBM Accurate Mathematical Library
3 * Written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
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
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19/*******************************************************************/
20/* */
21/* MODULE_NAME: branred.c */
22/* */
23/* FUNCTIONS: branred */
24/* */
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"
38#include <math.h>
39#include <math_private.h>
40
41#ifndef SECTION
42# define SECTION
43#endif
44
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/*******************************************************************/
52int
53SECTION
54__branred(double x, double *a, double *aa)
55{
56 int i,k;
57 mynumber u,gor;
58 double r[6],s,t,sum,b,bb,sum1,sum2,b1,bb1,b2,bb2,x1,x2,t1,t2;
59
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;
66 k = (u.i[HIGH_HALF]>>20)&2047;
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);
72 for (i=0;i<6;i++)
73 { r[i] = x1*toverp[k+i]*gor.x; gor.x *= tm24.x; }
74 for (i=0;i<3;i++) {
75 s=(r[i]+big.x)-big.x;
76 sum+=s;
77 r[i]-=s;
78 }
79 t=0;
80 for (i=0;i<6;i++)
81 t+=r[5-i];
82 bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5];
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;
94
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++) {
105 s=(r[i]+big.x)-big.x;
106 sum+=s;
107 r[i]-=s;
108 }
109 t=0;
110 for (i=0;i<6;i++)
111 t+=r[5-i];
112 bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5];
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;
120
121 b2=b;
122 bb2=bb;
123 sum2=sum;
124
125 sum=sum1+sum2;
126 b=b1+b2;
127 bb = (fabs(b1)>fabs(b2))? (b1-b)+b2 : (b2-b)+b1;
128 if (b > 0.5)
129 {b-=1.0; sum+=1.0;}
130 else if (b < -0.5)
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}
145