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 */
22static const unsigned long long oddplterm[HiPOWERop] =
23{
24  0x0000000000000000LL,
25  0x0051a1cf08fca228LL,
26  0x0000000071284ff7LL
27};
28
29#define	HiPOWERon	2	/* odd poly, negative terms */
30static const unsigned long long oddnegterm[HiPOWERon] =
31{
32   0x1291a9a184244e80LL,
33   0x0000583245819c21LL
34};
35
36#define	HiPOWERep	2	/* even poly, positive terms */
37static const unsigned long long evenplterm[HiPOWERep] =
38{
39  0x0e848884b539e888LL,
40  0x00003c7f18b887daLL
41};
42
43#define	HiPOWERen	2	/* even poly, negative terms */
44static const unsigned long long evennegterm[HiPOWERen] =
45{
46  0xf1f0200fd51569ccLL,
47  0x003afb46105c4432LL
48};
49
50static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
51
52
53/*--- poly_tan() ------------------------------------------------------------+
54 |                                                                           |
55 +---------------------------------------------------------------------------*/
56void	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