• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*---------------------------------------------------------------------------+
2  |  poly_tan.c                                                               |
3  |                                                                           |
4  | Compute the tan of a FPU_REG, using a polynomial approximation.           |
5  |                                                                           |
6  | Copyright (C) 1992,1993,1994,1997,1999                                    |
7  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
8  |                       Australia.  E-mail   billm@melbpc.org.au            |
9  |                                                                           |
10  |                                                                           |
11  +---------------------------------------------------------------------------*/
12 
13 #include "exception.h"
14 #include "reg_constant.h"
15 #include "fpu_emu.h"
16 #include "fpu_system.h"
17 #include "control_w.h"
18 #include "poly.h"
19 
20 #define	HiPOWERop	3	/* odd poly, positive terms */
21 static const unsigned long long oddplterm[HiPOWERop] = {
22 	0x0000000000000000LL,
23 	0x0051a1cf08fca228LL,
24 	0x0000000071284ff7LL
25 };
26 
27 #define	HiPOWERon	2	/* odd poly, negative terms */
28 static const unsigned long long oddnegterm[HiPOWERon] = {
29 	0x1291a9a184244e80LL,
30 	0x0000583245819c21LL
31 };
32 
33 #define	HiPOWERep	2	/* even poly, positive terms */
34 static const unsigned long long evenplterm[HiPOWERep] = {
35 	0x0e848884b539e888LL,
36 	0x00003c7f18b887daLL
37 };
38 
39 #define	HiPOWERen	2	/* even poly, negative terms */
40 static const unsigned long long evennegterm[HiPOWERen] = {
41 	0xf1f0200fd51569ccLL,
42 	0x003afb46105c4432LL
43 };
44 
45 static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
46 
47 /*--- poly_tan() ------------------------------------------------------------+
48  |                                                                           |
49  +---------------------------------------------------------------------------*/
poly_tan(FPU_REG * st0_ptr)50 void poly_tan(FPU_REG *st0_ptr)
51 {
52 	long int exponent;
53 	int invert;
54 	Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
55 	    argSignif, fix_up;
56 	unsigned long adj;
57 
58 	exponent = exponent(st0_ptr);
59 
60 #ifdef PARANOID
61 	if (signnegative(st0_ptr)) {	/* Can't hack a number < 0.0 */
62 		arith_invalid(0);
63 		return;
64 	}			/* Need a positive number */
65 #endif /* PARANOID */
66 
67 	/* Split the problem into two domains, smaller and larger than pi/4 */
68 	if ((exponent == 0)
69 	    || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
70 		/* The argument is greater than (approx) pi/4 */
71 		invert = 1;
72 		accum.lsw = 0;
73 		XSIG_LL(accum) = significand(st0_ptr);
74 
75 		if (exponent == 0) {
76 			/* The argument is >= 1.0 */
77 			/* Put the binary point at the left. */
78 			XSIG_LL(accum) <<= 1;
79 		}
80 		/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
81 		XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
82 		/* This is a special case which arises due to rounding. */
83 		if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
84 			FPU_settag0(TAG_Valid);
85 			significand(st0_ptr) = 0x8a51e04daabda360LL;
86 			setexponent16(st0_ptr,
87 				      (0x41 + EXTENDED_Ebias) | SIGN_Negative);
88 			return;
89 		}
90 
91 		argSignif.lsw = accum.lsw;
92 		XSIG_LL(argSignif) = XSIG_LL(accum);
93 		exponent = -1 + norm_Xsig(&argSignif);
94 	} else {
95 		invert = 0;
96 		argSignif.lsw = 0;
97 		XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
98 
99 		if (exponent < -1) {
100 			/* shift the argument right by the required places */
101 			if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >=
102 			    0x80000000U)
103 				XSIG_LL(accum)++;	/* round up */
104 		}
105 	}
106 
107 	XSIG_LL(argSq) = XSIG_LL(accum);
108 	argSq.lsw = accum.lsw;
109 	mul_Xsig_Xsig(&argSq, &argSq);
110 	XSIG_LL(argSqSq) = XSIG_LL(argSq);
111 	argSqSq.lsw = argSq.lsw;
112 	mul_Xsig_Xsig(&argSqSq, &argSqSq);
113 
114 	/* Compute the negative terms for the numerator polynomial */
115 	accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
116 	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
117 			HiPOWERon - 1);
118 	mul_Xsig_Xsig(&accumulatoro, &argSq);
119 	negate_Xsig(&accumulatoro);
120 	/* Add the positive terms */
121 	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
122 			HiPOWERop - 1);
123 
124 	/* Compute the positive terms for the denominator polynomial */
125 	accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
126 	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
127 			HiPOWERep - 1);
128 	mul_Xsig_Xsig(&accumulatore, &argSq);
129 	negate_Xsig(&accumulatore);
130 	/* Add the negative terms */
131 	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
132 			HiPOWERen - 1);
133 	/* Multiply by arg^2 */
134 	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
135 	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
136 	/* de-normalize and divide by 2 */
137 	shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
138 	negate_Xsig(&accumulatore);	/* This does 1 - accumulator */
139 
140 	/* Now find the ratio. */
141 	if (accumulatore.msw == 0) {
142 		/* accumulatoro must contain 1.0 here, (actually, 0) but it
143 		   really doesn't matter what value we use because it will
144 		   have negligible effect in later calculations
145 		 */
146 		XSIG_LL(accum) = 0x8000000000000000LL;
147 		accum.lsw = 0;
148 	} else {
149 		div_Xsig(&accumulatoro, &accumulatore, &accum);
150 	}
151 
152 	/* Multiply by 1/3 * arg^3 */
153 	mul64_Xsig(&accum, &XSIG_LL(argSignif));
154 	mul64_Xsig(&accum, &XSIG_LL(argSignif));
155 	mul64_Xsig(&accum, &XSIG_LL(argSignif));
156 	mul64_Xsig(&accum, &twothirds);
157 	shr_Xsig(&accum, -2 * (exponent + 1));
158 
159 	/* tan(arg) = arg + accum */
160 	add_two_Xsig(&accum, &argSignif, &exponent);
161 
162 	if (invert) {
163 		/* We now have the value of tan(pi_2 - arg) where pi_2 is an
164 		   approximation for pi/2
165 		 */
166 		/* The next step is to fix the answer to compensate for the
167 		   error due to the approximation used for pi/2
168 		 */
169 
170 		/* This is (approx) delta, the error in our approx for pi/2
171 		   (see above). It has an exponent of -65
172 		 */
173 		XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
174 		fix_up.lsw = 0;
175 
176 		if (exponent == 0)
177 			adj = 0xffffffff;	/* We want approx 1.0 here, but
178 						   this is close enough. */
179 		else if (exponent > -30) {
180 			adj = accum.msw >> -(exponent + 1);	/* tan */
181 			adj = mul_32_32(adj, adj);	/* tan^2 */
182 		} else
183 			adj = 0;
184 		adj = mul_32_32(0x898cc517, adj);	/* delta * tan^2 */
185 
186 		fix_up.msw += adj;
187 		if (!(fix_up.msw & 0x80000000)) {	/* did fix_up overflow ? */
188 			/* Yes, we need to add an msb */
189 			shr_Xsig(&fix_up, 1);
190 			fix_up.msw |= 0x80000000;
191 			shr_Xsig(&fix_up, 64 + exponent);
192 		} else
193 			shr_Xsig(&fix_up, 65 + exponent);
194 
195 		add_two_Xsig(&accum, &fix_up, &exponent);
196 
197 		/* accum now contains tan(pi/2 - arg).
198 		   Use tan(arg) = 1.0 / tan(pi/2 - arg)
199 		 */
200 		accumulatoro.lsw = accumulatoro.midw = 0;
201 		accumulatoro.msw = 0x80000000;
202 		div_Xsig(&accumulatoro, &accum, &accum);
203 		exponent = -exponent - 1;
204 	}
205 
206 	/* Transfer the result */
207 	round_Xsig(&accum);
208 	FPU_settag0(TAG_Valid);
209 	significand(st0_ptr) = XSIG_LL(accum);
210 	setexponent16(st0_ptr, exponent + EXTENDED_Ebias);	/* Result is positive. */
211 
212 }
213