PageRenderTime 39ms CodeModel.GetById 13ms app.highlight 20ms RepoModel.GetById 2ms app.codeStats 0ms

/arch/i386/math-emu/poly_tan.c

https://bitbucket.org/evzijst/gittest
C | 222 lines | 143 code | 30 blank | 49 comment | 20 complexity | 1c3d0f20cb3715e3fb21c48d9a9e5e06 MD5 | raw file
  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}