1 // SPDX-License-Identifier: GPL-2.0 1
2 /*--------------------------------------------
3 | poly_atan.c
4 |
5 | Compute the arctan of a FPU_REG, using a po
6 |
7 | Copyright (C) 1992,1993,1994,1997
8 | W. Metzenthen, 22 Parker S
9 | E-mail billm@suburbia.ne
10 |
11 |
12 +--------------------------------------------
13
14 #include "exception.h"
15 #include "reg_constant.h"
16 #include "fpu_emu.h"
17 #include "fpu_system.h"
18 #include "status_w.h"
19 #include "control_w.h"
20 #include "poly.h"
21
22 #define HIPOWERon 6 /* odd poly, n
23 static const unsigned long long oddnegterms[HI
24 0x0000000000000000LL, /* Dummy (not
25 0x015328437f756467LL,
26 0x0005dda27b73dec6LL,
27 0x0000226bf2bfb91aLL,
28 0x000000ccc439c5f7LL,
29 0x0000000355438407LL
30 };
31
32 #define HIPOWERop 6 /* odd poly, p
33 static const unsigned long long oddplterms[HIP
34 /* 0xaaaaaaaaaaaaaaabLL, transferred to fixe
35 0x0db55a71875c9ac2LL,
36 0x0029fce2d67880b0LL,
37 0x0000dfd3908b4596LL,
38 0x00000550fd61dab4LL,
39 0x0000001c9422b3f9LL,
40 0x000000003e3301e1LL
41 };
42
43 static const unsigned long long denomterm = 0x
44
45 static const Xsig fixedpterm = MK_XSIG(0xaaaaa
46
47 static const Xsig pi_signif = MK_XSIG(0xc90fda
48
49 /*--- poly_atan() ----------------------------
50 |
51 +--------------------------------------------
52 void poly_atan(FPU_REG *st0_ptr, u_char st0_ta
53 FPU_REG *st1_ptr, u_char st1_ta
54 {
55 u_char transformed, inverted, sign1, s
56 int exponent;
57 long int dummy_exp;
58 Xsig accumulator, Numer, Denom, accumu
59 u_char tag;
60
61 sign1 = getsign(st0_ptr);
62 sign2 = getsign(st1_ptr);
63 if (st0_tag == TAG_Valid) {
64 exponent = exponent(st0_ptr);
65 } else {
66 /* This gives non-compatible s
67 FPU_to_exp16(st0_ptr, st0_ptr)
68 exponent = exponent16(st0_ptr)
69 }
70 if (st1_tag == TAG_Valid) {
71 exponent -= exponent(st1_ptr);
72 } else {
73 /* This gives non-compatible s
74 FPU_to_exp16(st1_ptr, st1_ptr)
75 exponent -= exponent16(st1_ptr
76 }
77
78 if ((exponent < 0) || ((exponent == 0)
79 ((st0_ptr->sigh
80 ((st0_ptr->sig
81 (st0_ptr->sig
82 inverted = 1;
83 Numer.lsw = Denom.lsw = 0;
84 XSIG_LL(Numer) = significand(s
85 XSIG_LL(Denom) = significand(s
86 } else {
87 inverted = 0;
88 exponent = -exponent;
89 Numer.lsw = Denom.lsw = 0;
90 XSIG_LL(Numer) = significand(s
91 XSIG_LL(Denom) = significand(s
92 }
93 div_Xsig(&Numer, &Denom, &argSignif);
94 exponent += norm_Xsig(&argSignif);
95
96 if ((exponent >= -1)
97 || ((exponent == -2) && (argSignif
98 /* The argument is greater tha
99 /* Convert the argument by an
100 transformed = 1;
101
102 if (exponent >= 0) {
103 #ifdef PARANOID
104 if (!((exponent == 0)
105 (argSignif.lsw =
106 (argSignif.msw =
107 EXCEPTION(EX_I
108 return;
109 }
110 #endif /* PARANOID */
111 argSignif.msw = 0;
112 } else {
113 Numer.lsw = Denom.lsw
114 XSIG_LL(Numer) = XSIG_
115
116 if (exponent < -1)
117 shr_Xsig(&Nume
118 negate_Xsig(&Numer);
119
120 shr_Xsig(&Denom, -expo
121 Denom.msw |= 0x8000000
122
123 div_Xsig(&Numer, &Deno
124
125 exponent = -1 + norm_X
126 }
127 } else {
128 transformed = 0;
129 }
130
131 argSq.lsw = argSignif.lsw;
132 argSq.midw = argSignif.midw;
133 argSq.msw = argSignif.msw;
134 mul_Xsig_Xsig(&argSq, &argSq);
135
136 argSqSq.lsw = argSq.lsw;
137 argSqSq.midw = argSq.midw;
138 argSqSq.msw = argSq.msw;
139 mul_Xsig_Xsig(&argSqSq, &argSqSq);
140
141 accumulatore.lsw = argSq.lsw;
142 XSIG_LL(accumulatore) = XSIG_LL(argSq)
143
144 shr_Xsig(&argSq, 2 * (-1 - exponent -
145 shr_Xsig(&argSqSq, 4 * (-1 - exponent
146
147 /* Now have argSq etc with binary poin
148 .1xxxxxxxx */
149
150 /* Do the basic fixed point polynomial
151 accumulator.msw = accumulator.midw = a
152 polynomial_Xsig(&accumulator, &XSIG_LL
153 oddplterms, HIPOWERop
154 mul64_Xsig(&accumulator, &XSIG_LL(argS
155 negate_Xsig(&accumulator);
156 polynomial_Xsig(&accumulator, &XSIG_LL
157 HIPOWERon - 1);
158 negate_Xsig(&accumulator);
159 add_two_Xsig(&accumulator, &fixedpterm
160
161 mul64_Xsig(&accumulatore, &denomterm);
162 shr_Xsig(&accumulatore, 1 + 2 * (-1 -
163 accumulatore.msw |= 0x80000000;
164
165 div_Xsig(&accumulator, &accumulatore,
166
167 mul_Xsig_Xsig(&accumulator, &argSignif
168 mul_Xsig_Xsig(&accumulator, &argSq);
169
170 shr_Xsig(&accumulator, 3);
171 negate_Xsig(&accumulator);
172 add_Xsig_Xsig(&accumulator, &argSignif
173
174 if (transformed) {
175 /* compute pi/4 - accumulator
176 shr_Xsig(&accumulator, -1 - ex
177 negate_Xsig(&accumulator);
178 add_Xsig_Xsig(&accumulator, &p
179 exponent = -1;
180 }
181
182 if (inverted) {
183 /* compute pi/2 - accumulator
184 shr_Xsig(&accumulator, -expone
185 negate_Xsig(&accumulator);
186 add_Xsig_Xsig(&accumulator, &p
187 exponent = 0;
188 }
189
190 if (sign1) {
191 /* compute pi - accumulator */
192 shr_Xsig(&accumulator, 1 - exp
193 negate_Xsig(&accumulator);
194 add_Xsig_Xsig(&accumulator, &p
195 exponent = 1;
196 }
197
198 exponent += round_Xsig(&accumulator);
199
200 significand(st1_ptr) = XSIG_LL(accumul
201 setexponent16(st1_ptr, exponent);
202
203 tag = FPU_round(st1_ptr, 1, 0, FULL_PR
204 FPU_settagi(1, tag);
205
206 set_precision_flag_up(); /* We
207 use
208
209 }
210
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.