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 */
23 static 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 */
34 static 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 
45 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
46 
47 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
48 
49 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
50 
51 
52 /*--- poly_atan() -----------------------------------------------------------+
53  |                                                                           |
54  +---------------------------------------------------------------------------*/
poly_atan(FPU_REG * st0_ptr,u_char st0_tag,FPU_REG * st1_ptr,u_char st1_tag)55 void	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