1
/*---------------------------------------------------------------------------+
4
| Compute the arctan of a FPU_REG, using a polynomial approximation. |
6
| Copyright (C) 1992,1993,1994,1997 |
7
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
8
| E-mail billm@suburbia.net |
11
+---------------------------------------------------------------------------*/
13
#include "exception.h"
14
#include "reg_constant.h"
16
#include "fpu_system.h"
18
#include "control_w.h"
21
#define HIPOWERon 6 /* odd poly, negative terms */
22
static const unsigned long long oddnegterms[HIPOWERon] = {
23
0x0000000000000000LL, /* Dummy (not for - 1.0) */
31
#define HIPOWERop 6 /* odd poly, positive terms */
32
static const unsigned long long oddplterms[HIPOWERop] = {
33
/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
42
static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
44
static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
46
static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
48
/*--- poly_atan() -----------------------------------------------------------+
50
+---------------------------------------------------------------------------*/
51
void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
52
FPU_REG *st1_ptr, u_char st1_tag)
54
u_char transformed, inverted, sign1, sign2;
57
Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
60
sign1 = getsign(st0_ptr);
61
sign2 = getsign(st1_ptr);
62
if (st0_tag == TAG_Valid) {
63
exponent = exponent(st0_ptr);
65
/* This gives non-compatible stack contents... */
66
FPU_to_exp16(st0_ptr, st0_ptr);
67
exponent = exponent16(st0_ptr);
69
if (st1_tag == TAG_Valid) {
70
exponent -= exponent(st1_ptr);
72
/* This gives non-compatible stack contents... */
73
FPU_to_exp16(st1_ptr, st1_ptr);
74
exponent -= exponent16(st1_ptr);
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))))) {
82
Numer.lsw = Denom.lsw = 0;
83
XSIG_LL(Numer) = significand(st0_ptr);
84
XSIG_LL(Denom) = significand(st1_ptr);
88
Numer.lsw = Denom.lsw = 0;
89
XSIG_LL(Numer) = significand(st1_ptr);
90
XSIG_LL(Denom) = significand(st0_ptr);
92
div_Xsig(&Numer, &Denom, &argSignif);
93
exponent += norm_Xsig(&argSignif);
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 */
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 */
109
#endif /* PARANOID */
110
argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
112
Numer.lsw = Denom.lsw = argSignif.lsw;
113
XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
116
shr_Xsig(&Numer, -1 - exponent);
119
shr_Xsig(&Denom, -exponent);
120
Denom.msw |= 0x80000000;
122
div_Xsig(&Numer, &Denom, &argSignif);
124
exponent = -1 + norm_Xsig(&argSignif);
130
argSq.lsw = argSignif.lsw;
131
argSq.midw = argSignif.midw;
132
argSq.msw = argSignif.msw;
133
mul_Xsig_Xsig(&argSq, &argSq);
135
argSqSq.lsw = argSq.lsw;
136
argSqSq.midw = argSq.midw;
137
argSqSq.msw = argSq.msw;
138
mul_Xsig_Xsig(&argSqSq, &argSqSq);
140
accumulatore.lsw = argSq.lsw;
141
XSIG_LL(accumulatore) = XSIG_LL(argSq);
143
shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
144
shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
146
/* Now have argSq etc with binary point at the left
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,
157
negate_Xsig(&accumulator);
158
add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
160
mul64_Xsig(&accumulatore, &denomterm);
161
shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
162
accumulatore.msw |= 0x80000000;
164
div_Xsig(&accumulator, &accumulatore, &accumulator);
166
mul_Xsig_Xsig(&accumulator, &argSignif);
167
mul_Xsig_Xsig(&accumulator, &argSq);
169
shr_Xsig(&accumulator, 3);
170
negate_Xsig(&accumulator);
171
add_Xsig_Xsig(&accumulator, &argSignif);
174
/* compute pi/4 - accumulator */
175
shr_Xsig(&accumulator, -1 - exponent);
176
negate_Xsig(&accumulator);
177
add_Xsig_Xsig(&accumulator, &pi_signif);
182
/* compute pi/2 - accumulator */
183
shr_Xsig(&accumulator, -exponent);
184
negate_Xsig(&accumulator);
185
add_Xsig_Xsig(&accumulator, &pi_signif);
190
/* compute pi - accumulator */
191
shr_Xsig(&accumulator, 1 - exponent);
192
negate_Xsig(&accumulator);
193
add_Xsig_Xsig(&accumulator, &pi_signif);
197
exponent += round_Xsig(&accumulator);
199
significand(st1_ptr) = XSIG_LL(accumulator);
200
setexponent16(st1_ptr, exponent);
202
tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
205
set_precision_flag_up(); /* We do not really know if up or down,
206
use this as the default. */