~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

TOMOYO Linux Cross Reference
Linux/arch/x86/math-emu/poly_atan.c

Version: ~ [ linux-5.16-rc3 ] ~ [ linux-5.15.5 ] ~ [ linux-5.14.21 ] ~ [ linux-5.13.19 ] ~ [ linux-5.12.19 ] ~ [ linux-5.11.22 ] ~ [ linux-5.10.82 ] ~ [ linux-5.9.16 ] ~ [ linux-5.8.18 ] ~ [ linux-5.7.19 ] ~ [ linux-5.6.19 ] ~ [ linux-5.5.19 ] ~ [ linux-5.4.162 ] ~ [ linux-5.3.18 ] ~ [ linux-5.2.21 ] ~ [ linux-5.1.21 ] ~ [ linux-5.0.21 ] ~ [ linux-4.20.17 ] ~ [ linux-4.19.218 ] ~ [ linux-4.18.20 ] ~ [ linux-4.17.19 ] ~ [ linux-4.16.18 ] ~ [ linux-4.15.18 ] ~ [ linux-4.14.256 ] ~ [ linux-4.13.16 ] ~ [ linux-4.12.14 ] ~ [ linux-4.11.12 ] ~ [ linux-4.10.17 ] ~ [ linux-4.9.291 ] ~ [ linux-4.8.17 ] ~ [ linux-4.7.10 ] ~ [ linux-4.6.7 ] ~ [ linux-4.5.7 ] ~ [ linux-4.4.293 ] ~ [ linux-4.3.6 ] ~ [ linux-4.2.8 ] ~ [ linux-4.1.52 ] ~ [ linux-4.0.9 ] ~ [ linux-3.18.140 ] ~ [ linux-3.16.85 ] ~ [ linux-3.14.79 ] ~ [ linux-3.12.74 ] ~ [ linux-3.10.108 ] ~ [ linux-2.6.32.71 ] ~ [ linux-2.6.0 ] ~ [ linux-2.4.37.11 ] ~ [ unix-v6-master ] ~ [ ccs-tools-1.8.5 ] ~ [ policy-sample ] ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

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

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

kernel.org | git.kernel.org | LWN.net | Project Home | Wiki (Japanese) | Wiki (English) | SVN repository | Mail admin

Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.

osdn.jp