1/*---------------------------------------------------------------------------+
2 |  poly_atan.c                                                              |
3 |                                                                           |
4 | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
5 |                                                                           |
6 | Copyright (C) 1992,1993,1994,1997                                         |
7 |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
8 |                  E-mail   billm@suburbia.net                              |
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 "status_w.h"
18#include "control_w.h"
19#include "poly.h"
20
21
22#define	HIPOWERon	6	/* odd poly, negative terms */
23static const unsigned long long oddnegterms[HIPOWERon] =
24{
25  0x0000000000000000LL, /* Dummy (not for - 1.0) */
26  0x015328437f756467LL,
27  0x0005dda27b73dec6LL,
28  0x0000226bf2bfb91aLL,
29  0x000000ccc439c5f7LL,
30  0x0000000355438407LL
31} ;
32
33#define	HIPOWERop	6	/* odd poly, positive terms */
34static const unsigned long long oddplterms[HIPOWERop] =
35{
36/*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
37  0x0db55a71875c9ac2LL,
38  0x0029fce2d67880b0LL,
39  0x0000dfd3908b4596LL,
40  0x00000550fd61dab4LL,
41  0x0000001c9422b3f9LL,
42  0x000000003e3301e1LL
43};
44
45static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
46
47static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
48
49static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
50
51
52/*--- poly_atan() -----------------------------------------------------------+
53 |                                                                           |
54 +---------------------------------------------------------------------------*/
55void	poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
56		  FPU_REG *st1_ptr, u_char st1_tag)
57{
58  u_char	transformed, inverted,
59                sign1, sign2;
60  int           exponent;
61  long int   	dummy_exp;
62  Xsig          accumulator, Numer, Denom, accumulatore, argSignif,
63                argSq, argSqSq;
64  u_char        tag;
65
66  sign1 = getsign(st0_ptr);
67  sign2 = getsign(st1_ptr);
68  if ( st0_tag == TAG_Valid )
69    {
70      exponent = exponent(st0_ptr);
71    }
72  else
73    {
74      /* This gives non-compatible stack contents... */
75      FPU_to_exp16(st0_ptr, st0_ptr);
76      exponent = exponent16(st0_ptr);
77    }
78  if ( st1_tag == TAG_Valid )
79    {
80      exponent -= exponent(st1_ptr);
81    }
82  else
83    {
84      /* This gives non-compatible stack contents... */
85      FPU_to_exp16(st1_ptr, st1_ptr);
86      exponent -= exponent16(st1_ptr);
87    }
88
89  if ( (exponent < 0) || ((exponent == 0) &&
90			  ((st0_ptr->sigh < st1_ptr->sigh) ||
91			   ((st0_ptr->sigh == st1_ptr->sigh) &&
92			    (st0_ptr->sigl < st1_ptr->sigl))) ) )
93    {
94      inverted = 1;
95      Numer.lsw = Denom.lsw = 0;
96      XSIG_LL(Numer) = significand(st0_ptr);
97      XSIG_LL(Denom) = significand(st1_ptr);
98    }
99  else
100    {
101      inverted = 0;
102      exponent = -exponent;
103      Numer.lsw = Denom.lsw = 0;
104      XSIG_LL(Numer) = significand(st1_ptr);
105      XSIG_LL(Denom) = significand(st0_ptr);
106     }
107  div_Xsig(&Numer, &Denom, &argSignif);
108  exponent += norm_Xsig(&argSignif);
109
110  if ( (exponent >= -1)
111      || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
112    {
113      /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
114      /* Convert the argument by an identity for atan */
115      transformed = 1;
116
117      if ( exponent >= 0 )
118	{
119#ifdef PARANOID
120	  if ( !( (exponent == 0) &&
121		 (argSignif.lsw == 0) && (argSignif.midw == 0) &&
122		 (argSignif.msw == 0x80000000) ) )
123	    {
124	      EXCEPTION(EX_INTERNAL|0x104);  /* There must be a logic error */
125	      return;
126	    }
127#endif /* PARANOID */
128	  argSignif.msw = 0;   /* Make the transformed arg -> 0.0 */
129	}
130      else
131	{
132	  Numer.lsw = Denom.lsw = argSignif.lsw;
133	  XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
134
135	  if ( exponent < -1 )
136	    shr_Xsig(&Numer, -1-exponent);
137	  negate_Xsig(&Numer);
138
139	  shr_Xsig(&Denom, -exponent);
140	  Denom.msw |= 0x80000000;
141
142	  div_Xsig(&Numer, &Denom, &argSignif);
143
144	  exponent = -1 + norm_Xsig(&argSignif);
145	}
146    }
147  else
148    {
149      transformed = 0;
150    }
151
152  argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw;
153  argSq.msw = argSignif.msw;
154  mul_Xsig_Xsig(&argSq, &argSq);
155
156  argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw;
157  mul_Xsig_Xsig(&argSqSq, &argSqSq);
158
159  accumulatore.lsw = argSq.lsw;
160  XSIG_LL(accumulatore) = XSIG_LL(argSq);
161
162  shr_Xsig(&argSq, 2*(-1-exponent-1));
163  shr_Xsig(&argSqSq, 4*(-1-exponent-1));
164
165  /* Now have argSq etc with binary point at the left
166     .1xxxxxxxx */
167
168  /* Do the basic fixed point polynomial evaluation */
169  accumulator.msw = accumulator.midw = accumulator.lsw = 0;
170  polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
171		   oddplterms, HIPOWERop-1);
172  mul64_Xsig(&accumulator, &XSIG_LL(argSq));
173  negate_Xsig(&accumulator);
174  polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
175  negate_Xsig(&accumulator);
176  add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
177
178  mul64_Xsig(&accumulatore, &denomterm);
179  shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
180  accumulatore.msw |= 0x80000000;
181
182  div_Xsig(&accumulator, &accumulatore, &accumulator);
183
184  mul_Xsig_Xsig(&accumulator, &argSignif);
185  mul_Xsig_Xsig(&accumulator, &argSq);
186
187  shr_Xsig(&accumulator, 3);
188  negate_Xsig(&accumulator);
189  add_Xsig_Xsig(&accumulator, &argSignif);
190
191  if ( transformed )
192    {
193      /* compute pi/4 - accumulator */
194      shr_Xsig(&accumulator, -1-exponent);
195      negate_Xsig(&accumulator);
196      add_Xsig_Xsig(&accumulator, &pi_signif);
197      exponent = -1;
198    }
199
200  if ( inverted )
201    {
202      /* compute pi/2 - accumulator */
203      shr_Xsig(&accumulator, -exponent);
204      negate_Xsig(&accumulator);
205      add_Xsig_Xsig(&accumulator, &pi_signif);
206      exponent = 0;
207    }
208
209  if ( sign1 )
210    {
211      /* compute pi - accumulator */
212      shr_Xsig(&accumulator, 1 - exponent);
213      negate_Xsig(&accumulator);
214      add_Xsig_Xsig(&accumulator, &pi_signif);
215      exponent = 1;
216    }
217
218  exponent += round_Xsig(&accumulator);
219
220  significand(st1_ptr) = XSIG_LL(accumulator);
221  setexponent16(st1_ptr, exponent);
222
223  tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
224  FPU_settagi(1, tag);
225
226  set_precision_flag_up();  /* We do not really know if up or down,
227			       use this as the default. */
228
229}
230