1 /*
2 * poly_tan.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 HIPOWERop 3 /* odd poly, positive terms */
57 static unsigned short oddplterms[HIPOWERop][4] =
58 {
59 { 0x846a, 0x42d1, 0xb544, 0x921f},
60 { 0x6fb2, 0x0215, 0x95c0, 0x099c},
61 { 0xfce6, 0x0cc8, 0x1c9a, 0x0000}
62 };
63
64 #define HIPOWERon 2 /* odd poly, negative terms */
65 static unsigned short oddnegterms[HIPOWERon][4] =
66 {
67 { 0x6906, 0xe205, 0x25c8, 0x8838},
68 { 0x1dd7, 0x3fe3, 0x944e, 0x002c}
69 };
70
71 #define HIPOWERep 2 /* even poly, positive terms */
72 static unsigned short evenplterms[HIPOWERep][4] =
73 {
74 { 0xdb8f, 0x3761, 0x1432, 0x2acf},
75 { 0x16eb, 0x13c1, 0x3099, 0x0003}
76 };
77
78 #define HIPOWERen 2 /* even poly, negative terms */
79 static unsigned short evennegterms[HIPOWERen][4] =
80 {
81 { 0x3a7c, 0xe4c5, 0x7f87, 0x2945},
82 { 0x572b, 0x664c, 0xc543, 0x018c}
83 };
84
85
86 /*--- poly_tan() ------------------------------------------------------------+
87 | |
88 +---------------------------------------------------------------------------*/
poly_tan(FPU_REG * arg,FPU_REG * y_reg)89 void poly_tan(FPU_REG *arg, FPU_REG *y_reg)
90 {
91 char invert = 0;
92 short exponent;
93 FPU_REG odd_poly, even_poly, pos_poly, neg_poly;
94 FPU_REG argSq;
95 long long arg_signif, argSqSq;
96
97
98 exponent = arg->exp - EXP_BIAS;
99
100 if ( arg->tag == TW_Zero )
101 {
102 /* Return 0.0 */
103 reg_move(&CONST_Z, y_reg);
104 return;
105 }
106
107 if ( exponent >= -1 )
108 {
109 /* argument is in the range [0.5 .. 1.0] */
110 if ( exponent >= 0 )
111 {
112 #ifdef PARANOID
113 if ( (exponent == 0) &&
114 (arg->sigl == 0) && (arg->sigh == 0x80000000) )
115 #endif PARANOID
116 {
117 arith_overflow(y_reg);
118 return;
119 }
120 #ifdef PARANOID
121 EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
122 return;
123 #endif PARANOID
124 }
125 /* The argument is in the range [0.5 .. 1.0) */
126 /* Convert the argument to a number in the range (0.0 .. 0.5] */
127 *((long long *)(&arg->sigl)) = - *((long long *)(&arg->sigl));
128 normalize(arg); /* Needed later */
129 exponent = arg->exp - EXP_BIAS;
130 invert = 1;
131 }
132
133 #ifdef PARANOID
134 if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
135 { arith_invalid(y_reg); return; } /* Need a positive number */
136 #endif PARANOID
137
138 *(long long *)&arg_signif = *(long long *)&(arg->sigl);
139 if ( exponent < -1 )
140 {
141 /* shift the argument right by the required places */
142 if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
143 arg_signif++; /* round up */
144 }
145
146 mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
147 mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
148
149 /* will be a valid positive nr with expon = 0 */
150 *(short *)&(pos_poly.sign) = 0;
151 pos_poly.exp = EXP_BIAS;
152
153 /* Do the basic fixed point polynomial evaluation */
154 polynomial(&pos_poly.sigl, (unsigned *)&argSqSq, oddplterms, HIPOWERop-1);
155
156 /* will be a valid positive nr with expon = 0 */
157 *(short *)&(neg_poly.sign) = 0;
158 neg_poly.exp = EXP_BIAS;
159
160 /* Do the basic fixed point polynomial evaluation */
161 polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, oddnegterms, HIPOWERon-1);
162 mul64((long long *)(&argSq.sigl), (long long *)(&neg_poly.sigl),
163 (long long *)(&neg_poly.sigl));
164
165 /* Subtract the mantissas */
166 *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
167
168 /* Convert to 64 bit signed-compatible */
169 pos_poly.exp -= 1;
170
171 reg_move(&pos_poly, &odd_poly);
172 normalize(&odd_poly);
173
174 reg_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
175 reg_u_add(&odd_poly, arg, &odd_poly, FULL_PRECISION); /* This is just the odd polynomial */
176
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, evenplterms, HIPOWERep-1);
184 mul64((long long *)(&argSq.sigl),
185 (long long *)(&pos_poly.sigl), (long long *)(&pos_poly.sigl));
186
187 /* will be a valid positive nr with expon = 0 */
188 *(short *)&(neg_poly.sign) = 0;
189 neg_poly.exp = EXP_BIAS;
190
191 /* Do the basic fixed point polynomial evaluation */
192 polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, evennegterms, HIPOWERen-1);
193
194 /* Subtract the mantissas */
195 *((long long *)(&neg_poly.sigl)) -= *((long long *)(&pos_poly.sigl));
196 /* and multiply by argSq */
197
198 /* Convert argSq to a valid reg number */
199 *(short *)&(argSq.sign) = 0;
200 argSq.exp = EXP_BIAS - 1;
201 normalize(&argSq);
202
203 /* Convert to 64 bit signed-compatible */
204 neg_poly.exp -= 1;
205
206 reg_move(&neg_poly, &even_poly);
207 normalize(&even_poly);
208
209 reg_mul(&even_poly, &argSq, &even_poly, FULL_PRECISION);
210 reg_add(&even_poly, &argSq, &even_poly, FULL_PRECISION);
211 reg_sub(&CONST_1, &even_poly, &even_poly, FULL_PRECISION); /* This is just the even polynomial */
212
213 /* Now ready to copy the results */
214 if ( invert )
215 { reg_div(&even_poly, &odd_poly, y_reg, FULL_PRECISION); }
216 else
217 { reg_div(&odd_poly, &even_poly, y_reg, FULL_PRECISION); }
218
219 }
220