1 /*
2 * p_atan.c
3 *
4 * Compute the tan of a FPU_REG, using a polynomial approximation.
5 *
6 *
7 * Copyright (C) 1992, 1993 W. Metzenthen, 22 Parker St, Ormond,
8 * Vic 3163, Australia.
9 * E-mail apm233m@vaxc.cc.monash.edu.au
10 * All rights reserved.
11 *
12 * This copyright notice covers the redistribution and use of the
13 * FPU emulator developed by W. Metzenthen. It covers only its use
14 * in the 386BSD operating system. Any other use is not permitted
15 * under this copyright.
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 * 2. Redistributions in binary form must include information specifying
23 * that source code for the emulator is freely available and include
24 * either:
25 * a) an offer to provide the source code for a nominal distribution
26 * fee, or
27 * b) list at least two alternative methods whereby the source
28 * can be obtained, e.g. a publically accessible bulletin board
29 * and an anonymous ftp site from which the software can be
30 * downloaded.
31 * 3. All advertising materials specifically mentioning features or use of
32 * this emulator must acknowledge that it was developed by W. Metzenthen.
33 * 4. The name of W. Metzenthen may not be used to endorse or promote
34 * products derived from this software without specific prior written
35 * permission.
36 *
37 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
38 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
39 * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
40 * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
41 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
42 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
43 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
44 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
45 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
46 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
47 *
48 */
49
50 #include "exception.h"
51 #include "reg_constant.h"
52 #include "fpu_emu.h"
53 #include "control_w.h"
54
55
56 #define HIPOWERon 6 /* odd poly, negative terms */
57 static unsigned oddnegterms[HIPOWERon][2] =
58 {
59 { 0x00000000, 0x00000000 }, /* for + 1.0 */
60 { 0x763b6f3d, 0x1adc4428 },
61 { 0x20f0630b, 0x0502909d },
62 { 0x4e825578, 0x0198ce38 },
63 { 0x22b7cb87, 0x008da6e3 },
64 { 0x9b30ca03, 0x00239c79 }
65 } ;
66
67 #define HIPOWERop 6 /* odd poly, positive terms */
68 static unsigned oddplterms[HIPOWERop][2] =
69 {
70 { 0xa6f67cb8, 0x94d910bd },
71 { 0xa02ffab4, 0x0a43cb45 },
72 { 0x04265e6b, 0x02bf5655 },
73 { 0x0a728914, 0x00f280f7 },
74 { 0x6d640e01, 0x004d6556 },
75 { 0xf1dd2dbf, 0x000a530a }
76 };
77
78
79 static unsigned denomterm[2] =
80 { 0xfc4bd208, 0xea2e6612 };
81
82
83
84 /*--- poly_atan() -----------------------------------------------------------+
85 | |
86 +---------------------------------------------------------------------------*/
poly_atan(FPU_REG * arg)87 void poly_atan(FPU_REG *arg)
88 {
89 char recursions = 0;
90 short exponent;
91 FPU_REG odd_poly, even_poly, pos_poly, neg_poly;
92 FPU_REG argSq;
93 long long arg_signif, argSqSq;
94
95
96 #ifdef PARANOID
97 if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
98 { arith_invalid(arg); return; } /* Need a positive number */
99 #endif PARANOID
100
101 exponent = arg->exp - EXP_BIAS;
102
103 if ( arg->tag == TW_Zero )
104 {
105 /* Return 0.0 */
106 reg_move(&CONST_Z, arg);
107 return;
108 }
109
110 if ( exponent >= -2 )
111 {
112 /* argument is in the range [0.25 .. 1.0] */
113 if ( exponent >= 0 )
114 {
115 #ifdef PARANOID
116 if ( (exponent == 0) &&
117 (arg->sigl == 0) && (arg->sigh == 0x80000000) )
118 #endif PARANOID
119 {
120 reg_move(&CONST_PI4, arg);
121 return;
122 }
123 #ifdef PARANOID
124 EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
125 #endif PARANOID
126 }
127
128 /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
129 /* convert the argument by an identity for atan */
130 if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
131 {
132 FPU_REG numerator, denom;
133
134 recursions++;
135
136 arg_signif = *(long long *)&(arg->sigl);
137 if ( exponent < -1 )
138 {
139 if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
140 arg_signif++; /* round up */
141 }
142 *(long long *)&(numerator.sigl) = -arg_signif;
143 numerator.exp = EXP_BIAS - 1;
144 normalize(&numerator); /* 1 - arg */
145
146 arg_signif = *(long long *)&(arg->sigl);
147 if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
148 arg_signif++; /* round up */
149 *(long long *)&(denom.sigl) = arg_signif;
150 denom.sigh |= 0x80000000; /* 1 + arg */
151
152 arg->exp = numerator.exp;
153 reg_u_div(&numerator, &denom, arg, FULL_PRECISION);
154
155 exponent = arg->exp - EXP_BIAS;
156 }
157 }
158
159 *(long long *)&arg_signif = *(long long *)&(arg->sigl);
160
161 #ifdef PARANOID
162 /* This must always be true */
163 if ( exponent >= -1 )
164 {
165 EXCEPTION(EX_INTERNAL|0x120); /* There must be a logic error */
166 }
167 #endif PARANOID
168
169 /* shift the argument right by the required places */
170 if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
171 arg_signif++; /* round up */
172
173 /* Now have arg_signif with binary point at the left
174 .1xxxxxxxx */
175 mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
176 mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
177
178 /* will be a valid positive nr with expon = 0 */
179 *(short *)&(pos_poly.sign) = 0;
180 pos_poly.exp = EXP_BIAS;
181
182 /* Do the basic fixed point polynomial evaluation */
183 polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
184 (unsigned short (*)[4])oddplterms, HIPOWERop-1);
185 mul64((long long *)(&argSq.sigl), (long long *)(&pos_poly.sigl),
186 (long long *)(&pos_poly.sigl));
187
188 /* will be a valid positive nr with expon = 0 */
189 *(short *)&(neg_poly.sign) = 0;
190 neg_poly.exp = EXP_BIAS;
191
192 /* Do the basic fixed point polynomial evaluation */
193 polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
194 (unsigned short (*)[4])oddnegterms, HIPOWERon-1);
195
196 /* Subtract the mantissas */
197 *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
198
199 reg_move(&pos_poly, &odd_poly);
200 poly_add_1(&odd_poly);
201
202 /* The complete odd polynomial */
203 reg_u_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
204
205 /* will be a valid positive nr with expon = 0 */
206 *(short *)&(even_poly.sign) = 0;
207
208 mul64((long long *)(&argSq.sigl),
209 (long long *)(&denomterm), (long long *)(&even_poly.sigl));
210
211 poly_add_1(&even_poly);
212
213 reg_div(&odd_poly, &even_poly, arg, FULL_PRECISION);
214
215 if ( recursions )
216 reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION);
217 }
218
219
220 /* The argument to this function must be polynomial() compatible,
221 i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
222 be normalized.
223 This function adds 1.0 to the (assumed positive) argument. */
poly_add_1(FPU_REG * src)224 void poly_add_1(FPU_REG *src)
225 {
226 /* Rounding in a consistent direction produces better results
227 for the use of this function in poly_atan. Simple truncation
228 is used here instead of round-to-nearest. */
229
230 #ifdef OBSOLETE
231 char round = (src->sigl & 3) == 3;
232 #endif OBSOLETE
233
234 shrx(&src->sigl, 1);
235
236 #ifdef OBSOLETE
237 if ( round ) (*(long long *)&src->sigl)++; /* Round to even */
238 #endif OBSOLETE
239
240 src->sigh |= 0x80000000;
241
242 src->exp = EXP_BIAS;
243
244 }
245