Loading...
1/*---------------------------------------------------------------------------+
2 | poly_atan.c |
3 | |
4 | Compute the arctan of a FPU_REG, using a polynomial approximation. |
5 | |
6 | Copyright (C) 1992,1993,1994,1997 |
7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
8 | E-mail billm@suburbia.net |
9 | |
10 | |
11 +---------------------------------------------------------------------------*/
12
13#include "exception.h"
14#include "reg_constant.h"
15#include "fpu_emu.h"
16#include "fpu_system.h"
17#include "status_w.h"
18#include "control_w.h"
19#include "poly.h"
20
21#define HIPOWERon 6 /* odd poly, negative terms */
22static const unsigned long long oddnegterms[HIPOWERon] = {
23 0x0000000000000000LL, /* Dummy (not for - 1.0) */
24 0x015328437f756467LL,
25 0x0005dda27b73dec6LL,
26 0x0000226bf2bfb91aLL,
27 0x000000ccc439c5f7LL,
28 0x0000000355438407LL
29};
30
31#define HIPOWERop 6 /* odd poly, positive terms */
32static const unsigned long long oddplterms[HIPOWERop] = {
33/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
34 0x0db55a71875c9ac2LL,
35 0x0029fce2d67880b0LL,
36 0x0000dfd3908b4596LL,
37 0x00000550fd61dab4LL,
38 0x0000001c9422b3f9LL,
39 0x000000003e3301e1LL
40};
41
42static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
43
44static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
45
46static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
47
48/*--- poly_atan() -----------------------------------------------------------+
49 | |
50 +---------------------------------------------------------------------------*/
51void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
52 FPU_REG *st1_ptr, u_char st1_tag)
53{
54 u_char transformed, inverted, sign1, sign2;
55 int exponent;
56 long int dummy_exp;
57 Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
58 u_char tag;
59
60 sign1 = getsign(st0_ptr);
61 sign2 = getsign(st1_ptr);
62 if (st0_tag == TAG_Valid) {
63 exponent = exponent(st0_ptr);
64 } else {
65 /* This gives non-compatible stack contents... */
66 FPU_to_exp16(st0_ptr, st0_ptr);
67 exponent = exponent16(st0_ptr);
68 }
69 if (st1_tag == TAG_Valid) {
70 exponent -= exponent(st1_ptr);
71 } else {
72 /* This gives non-compatible stack contents... */
73 FPU_to_exp16(st1_ptr, st1_ptr);
74 exponent -= exponent16(st1_ptr);
75 }
76
77 if ((exponent < 0) || ((exponent == 0) &&
78 ((st0_ptr->sigh < st1_ptr->sigh) ||
79 ((st0_ptr->sigh == st1_ptr->sigh) &&
80 (st0_ptr->sigl < st1_ptr->sigl))))) {
81 inverted = 1;
82 Numer.lsw = Denom.lsw = 0;
83 XSIG_LL(Numer) = significand(st0_ptr);
84 XSIG_LL(Denom) = significand(st1_ptr);
85 } else {
86 inverted = 0;
87 exponent = -exponent;
88 Numer.lsw = Denom.lsw = 0;
89 XSIG_LL(Numer) = significand(st1_ptr);
90 XSIG_LL(Denom) = significand(st0_ptr);
91 }
92 div_Xsig(&Numer, &Denom, &argSignif);
93 exponent += norm_Xsig(&argSignif);
94
95 if ((exponent >= -1)
96 || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
97 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
98 /* Convert the argument by an identity for atan */
99 transformed = 1;
100
101 if (exponent >= 0) {
102#ifdef PARANOID
103 if (!((exponent == 0) &&
104 (argSignif.lsw == 0) && (argSignif.midw == 0) &&
105 (argSignif.msw == 0x80000000))) {
106 EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
107 return;
108 }
109#endif /* PARANOID */
110 argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
111 } else {
112 Numer.lsw = Denom.lsw = argSignif.lsw;
113 XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
114
115 if (exponent < -1)
116 shr_Xsig(&Numer, -1 - exponent);
117 negate_Xsig(&Numer);
118
119 shr_Xsig(&Denom, -exponent);
120 Denom.msw |= 0x80000000;
121
122 div_Xsig(&Numer, &Denom, &argSignif);
123
124 exponent = -1 + norm_Xsig(&argSignif);
125 }
126 } else {
127 transformed = 0;
128 }
129
130 argSq.lsw = argSignif.lsw;
131 argSq.midw = argSignif.midw;
132 argSq.msw = argSignif.msw;
133 mul_Xsig_Xsig(&argSq, &argSq);
134
135 argSqSq.lsw = argSq.lsw;
136 argSqSq.midw = argSq.midw;
137 argSqSq.msw = argSq.msw;
138 mul_Xsig_Xsig(&argSqSq, &argSqSq);
139
140 accumulatore.lsw = argSq.lsw;
141 XSIG_LL(accumulatore) = XSIG_LL(argSq);
142
143 shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
144 shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
145
146 /* Now have argSq etc with binary point at the left
147 .1xxxxxxxx */
148
149 /* Do the basic fixed point polynomial evaluation */
150 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
151 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
152 oddplterms, HIPOWERop - 1);
153 mul64_Xsig(&accumulator, &XSIG_LL(argSq));
154 negate_Xsig(&accumulator);
155 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
156 HIPOWERon - 1);
157 negate_Xsig(&accumulator);
158 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
159
160 mul64_Xsig(&accumulatore, &denomterm);
161 shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
162 accumulatore.msw |= 0x80000000;
163
164 div_Xsig(&accumulator, &accumulatore, &accumulator);
165
166 mul_Xsig_Xsig(&accumulator, &argSignif);
167 mul_Xsig_Xsig(&accumulator, &argSq);
168
169 shr_Xsig(&accumulator, 3);
170 negate_Xsig(&accumulator);
171 add_Xsig_Xsig(&accumulator, &argSignif);
172
173 if (transformed) {
174 /* compute pi/4 - accumulator */
175 shr_Xsig(&accumulator, -1 - exponent);
176 negate_Xsig(&accumulator);
177 add_Xsig_Xsig(&accumulator, &pi_signif);
178 exponent = -1;
179 }
180
181 if (inverted) {
182 /* compute pi/2 - accumulator */
183 shr_Xsig(&accumulator, -exponent);
184 negate_Xsig(&accumulator);
185 add_Xsig_Xsig(&accumulator, &pi_signif);
186 exponent = 0;
187 }
188
189 if (sign1) {
190 /* compute pi - accumulator */
191 shr_Xsig(&accumulator, 1 - exponent);
192 negate_Xsig(&accumulator);
193 add_Xsig_Xsig(&accumulator, &pi_signif);
194 exponent = 1;
195 }
196
197 exponent += round_Xsig(&accumulator);
198
199 significand(st1_ptr) = XSIG_LL(accumulator);
200 setexponent16(st1_ptr, exponent);
201
202 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
203 FPU_settagi(1, tag);
204
205 set_precision_flag_up(); /* We do not really know if up or down,
206 use this as the default. */
207
208}
1// SPDX-License-Identifier: GPL-2.0
2/*---------------------------------------------------------------------------+
3 | poly_atan.c |
4 | |
5 | Compute the arctan of a FPU_REG, using a polynomial approximation. |
6 | |
7 | Copyright (C) 1992,1993,1994,1997 |
8 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
9 | E-mail billm@suburbia.net |
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, negative terms */
23static const unsigned long long oddnegterms[HIPOWERon] = {
24 0x0000000000000000LL, /* Dummy (not for - 1.0) */
25 0x015328437f756467LL,
26 0x0005dda27b73dec6LL,
27 0x0000226bf2bfb91aLL,
28 0x000000ccc439c5f7LL,
29 0x0000000355438407LL
30};
31
32#define HIPOWERop 6 /* odd poly, positive terms */
33static const unsigned long long oddplterms[HIPOWERop] = {
34/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
35 0x0db55a71875c9ac2LL,
36 0x0029fce2d67880b0LL,
37 0x0000dfd3908b4596LL,
38 0x00000550fd61dab4LL,
39 0x0000001c9422b3f9LL,
40 0x000000003e3301e1LL
41};
42
43static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
44
45static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
46
47static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
48
49/*--- poly_atan() -----------------------------------------------------------+
50 | |
51 +---------------------------------------------------------------------------*/
52void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
53 FPU_REG *st1_ptr, u_char st1_tag)
54{
55 u_char transformed, inverted, sign1, sign2;
56 int exponent;
57 long int dummy_exp;
58 Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
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 stack contents... */
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 stack contents... */
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 < st1_ptr->sigh) ||
80 ((st0_ptr->sigh == st1_ptr->sigh) &&
81 (st0_ptr->sigl < st1_ptr->sigl))))) {
82 inverted = 1;
83 Numer.lsw = Denom.lsw = 0;
84 XSIG_LL(Numer) = significand(st0_ptr);
85 XSIG_LL(Denom) = significand(st1_ptr);
86 } else {
87 inverted = 0;
88 exponent = -exponent;
89 Numer.lsw = Denom.lsw = 0;
90 XSIG_LL(Numer) = significand(st1_ptr);
91 XSIG_LL(Denom) = significand(st0_ptr);
92 }
93 div_Xsig(&Numer, &Denom, &argSignif);
94 exponent += norm_Xsig(&argSignif);
95
96 if ((exponent >= -1)
97 || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
98 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
99 /* Convert the argument by an identity for atan */
100 transformed = 1;
101
102 if (exponent >= 0) {
103#ifdef PARANOID
104 if (!((exponent == 0) &&
105 (argSignif.lsw == 0) && (argSignif.midw == 0) &&
106 (argSignif.msw == 0x80000000))) {
107 EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
108 return;
109 }
110#endif /* PARANOID */
111 argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
112 } else {
113 Numer.lsw = Denom.lsw = argSignif.lsw;
114 XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
115
116 if (exponent < -1)
117 shr_Xsig(&Numer, -1 - exponent);
118 negate_Xsig(&Numer);
119
120 shr_Xsig(&Denom, -exponent);
121 Denom.msw |= 0x80000000;
122
123 div_Xsig(&Numer, &Denom, &argSignif);
124
125 exponent = -1 + norm_Xsig(&argSignif);
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 - 1));
145 shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
146
147 /* Now have argSq etc with binary point at the left
148 .1xxxxxxxx */
149
150 /* Do the basic fixed point polynomial evaluation */
151 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
152 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
153 oddplterms, HIPOWERop - 1);
154 mul64_Xsig(&accumulator, &XSIG_LL(argSq));
155 negate_Xsig(&accumulator);
156 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
157 HIPOWERon - 1);
158 negate_Xsig(&accumulator);
159 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
160
161 mul64_Xsig(&accumulatore, &denomterm);
162 shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
163 accumulatore.msw |= 0x80000000;
164
165 div_Xsig(&accumulator, &accumulatore, &accumulator);
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 - exponent);
177 negate_Xsig(&accumulator);
178 add_Xsig_Xsig(&accumulator, &pi_signif);
179 exponent = -1;
180 }
181
182 if (inverted) {
183 /* compute pi/2 - accumulator */
184 shr_Xsig(&accumulator, -exponent);
185 negate_Xsig(&accumulator);
186 add_Xsig_Xsig(&accumulator, &pi_signif);
187 exponent = 0;
188 }
189
190 if (sign1) {
191 /* compute pi - accumulator */
192 shr_Xsig(&accumulator, 1 - exponent);
193 negate_Xsig(&accumulator);
194 add_Xsig_Xsig(&accumulator, &pi_signif);
195 exponent = 1;
196 }
197
198 exponent += round_Xsig(&accumulator);
199
200 significand(st1_ptr) = XSIG_LL(accumulator);
201 setexponent16(st1_ptr, exponent);
202
203 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
204 FPU_settagi(1, tag);
205
206 set_precision_flag_up(); /* We do not really know if up or down,
207 use this as the default. */
208
209}