Commit | Line | Data |
---|---|---|
f24c459c JH |
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 | #define HIPOWERop 6 /* odd poly, positive terms */ | |
67 | static unsigned oddplterms[HIPOWERop][2] = | |
68 | { | |
69 | {0xa6f67cb8, 0x94d910bd}, | |
70 | {0xa02ffab4, 0x0a43cb45}, | |
71 | {0x04265e6b, 0x02bf5655}, | |
72 | {0x0a728914, 0x00f280f7}, | |
73 | {0x6d640e01, 0x004d6556}, | |
74 | {0xf1dd2dbf, 0x000a530a} | |
75 | }; | |
76 | ||
77 | ||
78 | static unsigned denomterm[2] = | |
79 | {0xfc4bd208, 0xea2e6612}; | |
80 | ||
81 | ||
82 | ||
83 | /*--- poly_atan() -----------------------------------------------------------+ | |
84 | | | | |
85 | +---------------------------------------------------------------------------*/ | |
86 | void | |
87 | 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); | |
99 | return; | |
100 | } /* Need a positive number */ | |
101 | #endif /* PARANOID */ | |
102 | ||
103 | exponent = arg->exp - EXP_BIAS; | |
104 | ||
105 | if (arg->tag == TW_Zero) { | |
106 | /* Return 0.0 */ | |
107 | reg_move(&CONST_Z, arg); | |
108 | return; | |
109 | } | |
110 | if (exponent >= -2) { | |
111 | /* argument is in the range [0.25 .. 1.0] */ | |
112 | if (exponent >= 0) { | |
113 | #ifdef PARANOID | |
114 | if ((exponent == 0) && | |
115 | (arg->sigl == 0) && (arg->sigh == 0x80000000)) | |
116 | #endif /* PARANOID */ | |
117 | { | |
118 | reg_move(&CONST_PI4, arg); | |
119 | return; | |
120 | } | |
121 | #ifdef PARANOID | |
122 | EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic | |
123 | * error */ | |
124 | #endif /* PARANOID */ | |
125 | } | |
126 | /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */ | |
127 | /* convert the argument by an identity for atan */ | |
128 | if ((exponent >= -1) || (arg->sigh > 0xd413ccd0)) { | |
129 | FPU_REG numerator, denom; | |
130 | ||
131 | recursions++; | |
132 | ||
133 | arg_signif = *(long long *) &(arg->sigl); | |
134 | if (exponent < -1) { | |
135 | if (shrx(&arg_signif, -1 - exponent) >= (unsigned)0x80000000) | |
136 | arg_signif++; /* round up */ | |
137 | } | |
138 | *(long long *) &(numerator.sigl) = -arg_signif; | |
139 | numerator.exp = EXP_BIAS - 1; | |
140 | normalize(&numerator); /* 1 - arg */ | |
141 | ||
142 | arg_signif = *(long long *) &(arg->sigl); | |
143 | if (shrx(&arg_signif, -exponent) >= (unsigned)0x80000000) | |
144 | arg_signif++; /* round up */ | |
145 | *(long long *) &(denom.sigl) = arg_signif; | |
146 | denom.sigh |= 0x80000000; /* 1 + arg */ | |
147 | ||
148 | arg->exp = numerator.exp; | |
149 | reg_u_div(&numerator, &denom, arg, FULL_PRECISION); | |
150 | ||
151 | exponent = arg->exp - EXP_BIAS; | |
152 | } | |
153 | } | |
154 | *(long long *) &arg_signif = *(long long *) &(arg->sigl); | |
155 | ||
156 | #ifdef PARANOID | |
157 | /* This must always be true */ | |
158 | if (exponent >= -1) { | |
159 | EXCEPTION(EX_INTERNAL | 0x120); /* There must be a logic error */ | |
160 | } | |
161 | #endif /* PARANOID */ | |
162 | ||
163 | /* shift the argument right by the required places */ | |
164 | if (shrx(&arg_signif, -1 - exponent) >= (unsigned)0x80000000) | |
165 | arg_signif++; /* round up */ | |
166 | ||
167 | /* Now have arg_signif with binary point at the left .1xxxxxxxx */ | |
168 | mul64(&arg_signif, &arg_signif, (long long *) (&argSq.sigl)); | |
169 | mul64((long long *) (&argSq.sigl), (long long *) (&argSq.sigl), &argSqSq); | |
170 | ||
171 | /* will be a valid positive nr with expon = 0 */ | |
172 | *(short *) &(pos_poly.sign) = 0; | |
173 | pos_poly.exp = EXP_BIAS; | |
174 | ||
175 | /* Do the basic fixed point polynomial evaluation */ | |
176 | polynomial(&pos_poly.sigl, (unsigned *) &argSqSq, | |
177 | (unsigned short (*)[4]) oddplterms, HIPOWERop - 1); | |
178 | mul64((long long *) (&argSq.sigl), (long long *) (&pos_poly.sigl), | |
179 | (long long *) (&pos_poly.sigl)); | |
180 | ||
181 | /* will be a valid positive nr with expon = 0 */ | |
182 | *(short *) &(neg_poly.sign) = 0; | |
183 | neg_poly.exp = EXP_BIAS; | |
184 | ||
185 | /* Do the basic fixed point polynomial evaluation */ | |
186 | polynomial(&neg_poly.sigl, (unsigned *) &argSqSq, | |
187 | (unsigned short (*)[4]) oddnegterms, HIPOWERon - 1); | |
188 | ||
189 | /* Subtract the mantissas */ | |
190 | *((long long *) (&pos_poly.sigl)) -= *((long long *) (&neg_poly.sigl)); | |
191 | ||
192 | reg_move(&pos_poly, &odd_poly); | |
193 | poly_add_1(&odd_poly); | |
194 | ||
195 | /* The complete odd polynomial */ | |
196 | reg_u_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION); | |
197 | ||
198 | /* will be a valid positive nr with expon = 0 */ | |
199 | *(short *) &(even_poly.sign) = 0; | |
200 | ||
201 | mul64((long long *) (&argSq.sigl), | |
202 | (long long *) (&denomterm), (long long *) (&even_poly.sigl)); | |
203 | ||
204 | poly_add_1(&even_poly); | |
205 | ||
206 | reg_div(&odd_poly, &even_poly, arg, FULL_PRECISION); | |
207 | ||
208 | if (recursions) | |
209 | reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION); | |
210 | } | |
211 | ||
212 | ||
213 | /* The argument to this function must be polynomial() compatible, | |
214 | i.e. have an exponent (not checked) of EXP_BIAS-1 but need not | |
215 | be normalized. | |
216 | This function adds 1.0 to the (assumed positive) argument. */ | |
217 | void | |
218 | poly_add_1(FPU_REG * src) | |
219 | { | |
220 | /* Rounding in a consistent direction produces better results | |
221 | for the use of this function in poly_atan. Simple truncation | |
222 | is used here instead of round-to-nearest. */ | |
223 | ||
224 | #ifdef OBSOLETE | |
225 | char round = (src->sigl & 3) == 3; | |
226 | #endif /* OBSOLETE */ | |
227 | ||
228 | shrx(&src->sigl, 1); | |
229 | ||
230 | #ifdef OBSOLETE | |
231 | if (round) | |
232 | (*(long long *) &src->sigl)++; /* Round to even */ | |
233 | #endif /* OBSOLETE */ | |
234 | ||
235 | src->sigh |= 0x80000000; | |
236 | ||
237 | src->exp = EXP_BIAS; | |
238 | ||
239 | } |