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