1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
/*---------------------------------------------------------------------------+
| p_atan.c |
| |
| Compute the tan of a FPU_REG, using a polynomial approximation. |
| |
| Copyright (C) 1992,1993 |
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
| Australia. E-mail billm@vaxc.cc.monash.edu.au |
| |
| |
+---------------------------------------------------------------------------*/
#include "exception.h"
#include "reg_constant.h"
#include "fpu_emu.h"
#include "control_w.h"
#define HIPOWERon 6 /* odd poly, negative terms */
static unsigned const oddnegterms[HIPOWERon][2] =
{
{ 0x00000000, 0x00000000 }, /* for + 1.0 */
{ 0x763b6f3d, 0x1adc4428 },
{ 0x20f0630b, 0x0502909d },
{ 0x4e825578, 0x0198ce38 },
{ 0x22b7cb87, 0x008da6e3 },
{ 0x9b30ca03, 0x00239c79 }
} ;
#define HIPOWERop 6 /* odd poly, positive terms */
static unsigned const oddplterms[HIPOWERop][2] =
{
{ 0xa6f67cb8, 0x94d910bd },
{ 0xa02ffab4, 0x0a43cb45 },
{ 0x04265e6b, 0x02bf5655 },
{ 0x0a728914, 0x00f280f7 },
{ 0x6d640e01, 0x004d6556 },
{ 0xf1dd2dbf, 0x000a530a }
};
static unsigned long long const denomterm = 0xea2e6612fc4bd208LL;
/*--- poly_atan() -----------------------------------------------------------+
| |
+---------------------------------------------------------------------------*/
void poly_atan(FPU_REG *arg)
{
char recursions = 0;
short exponent;
FPU_REG odd_poly, even_poly, pos_poly, neg_poly, ratio;
FPU_REG argSq;
unsigned long long arg_signif, argSqSq;
#ifdef PARANOID
if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
{ arith_invalid(arg); return; } /* Need a positive number */
#endif PARANOID
exponent = arg->exp - EXP_BIAS;
if ( arg->tag == TW_Zero )
{
/* Return 0.0 */
reg_move(&CONST_Z, arg);
return;
}
if ( exponent >= -2 )
{
/* argument is in the range [0.25 .. 1.0] */
if ( exponent >= 0 )
{
#ifdef PARANOID
if ( (exponent == 0) &&
(arg->sigl == 0) && (arg->sigh == 0x80000000) )
#endif PARANOID
{
reg_move(&CONST_PI4, arg);
return;
}
#ifdef PARANOID
EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
return;
#endif PARANOID
}
/* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
/* convert the argument by an identity for atan */
if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
{
FPU_REG numerator, denom;
recursions++;
arg_signif = significand(arg);
if ( exponent < -1 )
{
if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
arg_signif++; /* round up */
}
significand(&numerator) = -arg_signif;
numerator.exp = EXP_BIAS - 1;
normalize(&numerator); /* 1 - arg */
arg_signif = significand(arg);
if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
arg_signif++; /* round up */
significand(&denom) = arg_signif;
denom.sigh |= 0x80000000; /* 1 + arg */
arg->exp = numerator.exp;
reg_u_div(&numerator, &denom, arg, FULL_PRECISION);
exponent = arg->exp - EXP_BIAS;
}
}
arg_signif = significand(arg);
#ifdef PARANOID
/* This must always be true */
if ( exponent >= -1 )
{
EXCEPTION(EX_INTERNAL|0x120); /* There must be a logic error */
}
#endif PARANOID
/* shift the argument right by the required places */
if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
arg_signif++; /* round up */
/* Now have arg_signif with binary point at the left
.1xxxxxxxx */
mul64(&arg_signif, &arg_signif, &significand(&argSq));
mul64(&significand(&argSq), &significand(&argSq), &argSqSq);
/* will be a valid positive nr with expon = 0 */
*(short *)&(pos_poly.sign) = 0;
pos_poly.exp = EXP_BIAS;
/* Do the basic fixed point polynomial evaluation */
polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
(unsigned short (*)[4])oddplterms, HIPOWERop-1);
mul64(&significand(&argSq), &significand(&pos_poly),
&significand(&pos_poly));
/* will be a valid positive nr with expon = 0 */
*(short *)&(neg_poly.sign) = 0;
neg_poly.exp = EXP_BIAS;
/* Do the basic fixed point polynomial evaluation */
polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
(unsigned short (*)[4])oddnegterms, HIPOWERon-1);
/* Subtract the mantissas */
significand(&pos_poly) -= significand(&neg_poly);
reg_move(&pos_poly, &odd_poly);
poly_add_1(&odd_poly);
/* will be a valid positive nr with expon = 0 */
*(short *)&(even_poly.sign) = 0;
mul64(&significand(&argSq), &denomterm, &significand(&even_poly));
poly_add_1(&even_poly);
reg_div(&odd_poly, &even_poly, &ratio, FULL_PRECISION);
reg_u_mul(&ratio, arg, arg, FULL_PRECISION);
if ( recursions )
reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION);
}
/* The argument to this function must be polynomial() compatible,
i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
be normalized.
This function adds 1.0 to the (assumed positive) argument. */
void poly_add_1(FPU_REG *src)
{
/* Rounding in a consistent direction produces better results
for the use of this function in poly_atan. Simple truncation
is used here instead of round-to-nearest. */
#ifdef OBSOLETE
char round = (src->sigl & 3) == 3;
#endif OBSOLETE
shrx(&src->sigl, 1);
#ifdef OBSOLETE
if ( round ) significand(src)++; /* Round to even */
#endif OBSOLETE
src->sigh |= 0x80000000;
src->exp = EXP_BIAS;
}