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