1b2441318SGreg Kroah-Hartman // SPDX-License-Identifier: GPL-2.0
2da957e11SThomas Gleixner /*---------------------------------------------------------------------------+
3da957e11SThomas Gleixner | poly_atan.c |
4da957e11SThomas Gleixner | |
5da957e11SThomas Gleixner | Compute the arctan of a FPU_REG, using a polynomial approximation. |
6da957e11SThomas Gleixner | |
7da957e11SThomas Gleixner | Copyright (C) 1992,1993,1994,1997 |
8da957e11SThomas Gleixner | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
9da957e11SThomas Gleixner | E-mail billm@suburbia.net |
10da957e11SThomas Gleixner | |
11da957e11SThomas Gleixner | |
12da957e11SThomas Gleixner +---------------------------------------------------------------------------*/
13da957e11SThomas Gleixner
14da957e11SThomas Gleixner #include "exception.h"
15da957e11SThomas Gleixner #include "reg_constant.h"
16da957e11SThomas Gleixner #include "fpu_emu.h"
17da957e11SThomas Gleixner #include "fpu_system.h"
18da957e11SThomas Gleixner #include "status_w.h"
19da957e11SThomas Gleixner #include "control_w.h"
20da957e11SThomas Gleixner #include "poly.h"
21da957e11SThomas Gleixner
22da957e11SThomas Gleixner #define HIPOWERon 6 /* odd poly, negative terms */
233d0d14f9SIngo Molnar static const unsigned long long oddnegterms[HIPOWERon] = {
24da957e11SThomas Gleixner 0x0000000000000000LL, /* Dummy (not for - 1.0) */
25da957e11SThomas Gleixner 0x015328437f756467LL,
26da957e11SThomas Gleixner 0x0005dda27b73dec6LL,
27da957e11SThomas Gleixner 0x0000226bf2bfb91aLL,
28da957e11SThomas Gleixner 0x000000ccc439c5f7LL,
29da957e11SThomas Gleixner 0x0000000355438407LL
30da957e11SThomas Gleixner };
31da957e11SThomas Gleixner
32da957e11SThomas Gleixner #define HIPOWERop 6 /* odd poly, positive terms */
333d0d14f9SIngo Molnar static const unsigned long long oddplterms[HIPOWERop] = {
34da957e11SThomas Gleixner /* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
35da957e11SThomas Gleixner 0x0db55a71875c9ac2LL,
36da957e11SThomas Gleixner 0x0029fce2d67880b0LL,
37da957e11SThomas Gleixner 0x0000dfd3908b4596LL,
38da957e11SThomas Gleixner 0x00000550fd61dab4LL,
39da957e11SThomas Gleixner 0x0000001c9422b3f9LL,
40da957e11SThomas Gleixner 0x000000003e3301e1LL
41da957e11SThomas Gleixner };
42da957e11SThomas Gleixner
43da957e11SThomas Gleixner static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
44da957e11SThomas Gleixner
45da957e11SThomas Gleixner static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
46da957e11SThomas Gleixner
47da957e11SThomas Gleixner static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
48da957e11SThomas Gleixner
49da957e11SThomas Gleixner /*--- poly_atan() -----------------------------------------------------------+
50da957e11SThomas Gleixner | |
51da957e11SThomas Gleixner +---------------------------------------------------------------------------*/
poly_atan(FPU_REG * st0_ptr,u_char st0_tag,FPU_REG * st1_ptr,u_char st1_tag)52da957e11SThomas Gleixner void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
53da957e11SThomas Gleixner FPU_REG *st1_ptr, u_char st1_tag)
54da957e11SThomas Gleixner {
553d0d14f9SIngo Molnar u_char transformed, inverted, sign1, sign2;
56da957e11SThomas Gleixner int exponent;
57da957e11SThomas Gleixner long int dummy_exp;
583d0d14f9SIngo Molnar Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
59da957e11SThomas Gleixner u_char tag;
60da957e11SThomas Gleixner
61da957e11SThomas Gleixner sign1 = getsign(st0_ptr);
62da957e11SThomas Gleixner sign2 = getsign(st1_ptr);
633d0d14f9SIngo Molnar if (st0_tag == TAG_Valid) {
64da957e11SThomas Gleixner exponent = exponent(st0_ptr);
653d0d14f9SIngo Molnar } else {
66da957e11SThomas Gleixner /* This gives non-compatible stack contents... */
67da957e11SThomas Gleixner FPU_to_exp16(st0_ptr, st0_ptr);
68da957e11SThomas Gleixner exponent = exponent16(st0_ptr);
69da957e11SThomas Gleixner }
703d0d14f9SIngo Molnar if (st1_tag == TAG_Valid) {
71da957e11SThomas Gleixner exponent -= exponent(st1_ptr);
723d0d14f9SIngo Molnar } else {
73da957e11SThomas Gleixner /* This gives non-compatible stack contents... */
74da957e11SThomas Gleixner FPU_to_exp16(st1_ptr, st1_ptr);
75da957e11SThomas Gleixner exponent -= exponent16(st1_ptr);
76da957e11SThomas Gleixner }
77da957e11SThomas Gleixner
78da957e11SThomas Gleixner if ((exponent < 0) || ((exponent == 0) &&
79da957e11SThomas Gleixner ((st0_ptr->sigh < st1_ptr->sigh) ||
80da957e11SThomas Gleixner ((st0_ptr->sigh == st1_ptr->sigh) &&
813d0d14f9SIngo Molnar (st0_ptr->sigl < st1_ptr->sigl))))) {
82da957e11SThomas Gleixner inverted = 1;
83da957e11SThomas Gleixner Numer.lsw = Denom.lsw = 0;
84da957e11SThomas Gleixner XSIG_LL(Numer) = significand(st0_ptr);
85da957e11SThomas Gleixner XSIG_LL(Denom) = significand(st1_ptr);
863d0d14f9SIngo Molnar } else {
87da957e11SThomas Gleixner inverted = 0;
88da957e11SThomas Gleixner exponent = -exponent;
89da957e11SThomas Gleixner Numer.lsw = Denom.lsw = 0;
90da957e11SThomas Gleixner XSIG_LL(Numer) = significand(st1_ptr);
91da957e11SThomas Gleixner XSIG_LL(Denom) = significand(st0_ptr);
92da957e11SThomas Gleixner }
93da957e11SThomas Gleixner div_Xsig(&Numer, &Denom, &argSignif);
94da957e11SThomas Gleixner exponent += norm_Xsig(&argSignif);
95da957e11SThomas Gleixner
96da957e11SThomas Gleixner if ((exponent >= -1)
973d0d14f9SIngo Molnar || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
98da957e11SThomas Gleixner /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
99da957e11SThomas Gleixner /* Convert the argument by an identity for atan */
100da957e11SThomas Gleixner transformed = 1;
101da957e11SThomas Gleixner
1023d0d14f9SIngo Molnar if (exponent >= 0) {
103da957e11SThomas Gleixner #ifdef PARANOID
104da957e11SThomas Gleixner if (!((exponent == 0) &&
105da957e11SThomas Gleixner (argSignif.lsw == 0) && (argSignif.midw == 0) &&
1063d0d14f9SIngo Molnar (argSignif.msw == 0x80000000))) {
107da957e11SThomas Gleixner EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
108da957e11SThomas Gleixner return;
109da957e11SThomas Gleixner }
110da957e11SThomas Gleixner #endif /* PARANOID */
111da957e11SThomas Gleixner argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
1123d0d14f9SIngo Molnar } else {
113da957e11SThomas Gleixner Numer.lsw = Denom.lsw = argSignif.lsw;
114da957e11SThomas Gleixner XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
115da957e11SThomas Gleixner
116da957e11SThomas Gleixner if (exponent < -1)
117da957e11SThomas Gleixner shr_Xsig(&Numer, -1 - exponent);
118da957e11SThomas Gleixner negate_Xsig(&Numer);
119da957e11SThomas Gleixner
120da957e11SThomas Gleixner shr_Xsig(&Denom, -exponent);
121da957e11SThomas Gleixner Denom.msw |= 0x80000000;
122da957e11SThomas Gleixner
123da957e11SThomas Gleixner div_Xsig(&Numer, &Denom, &argSignif);
124da957e11SThomas Gleixner
125da957e11SThomas Gleixner exponent = -1 + norm_Xsig(&argSignif);
126da957e11SThomas Gleixner }
1273d0d14f9SIngo Molnar } else {
128da957e11SThomas Gleixner transformed = 0;
129da957e11SThomas Gleixner }
130da957e11SThomas Gleixner
1313d0d14f9SIngo Molnar argSq.lsw = argSignif.lsw;
1323d0d14f9SIngo Molnar argSq.midw = argSignif.midw;
133da957e11SThomas Gleixner argSq.msw = argSignif.msw;
134da957e11SThomas Gleixner mul_Xsig_Xsig(&argSq, &argSq);
135da957e11SThomas Gleixner
1363d0d14f9SIngo Molnar argSqSq.lsw = argSq.lsw;
1373d0d14f9SIngo Molnar argSqSq.midw = argSq.midw;
1383d0d14f9SIngo Molnar argSqSq.msw = argSq.msw;
139da957e11SThomas Gleixner mul_Xsig_Xsig(&argSqSq, &argSqSq);
140da957e11SThomas Gleixner
141da957e11SThomas Gleixner accumulatore.lsw = argSq.lsw;
142da957e11SThomas Gleixner XSIG_LL(accumulatore) = XSIG_LL(argSq);
143da957e11SThomas Gleixner
144da957e11SThomas Gleixner shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
145da957e11SThomas Gleixner shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
146da957e11SThomas Gleixner
147da957e11SThomas Gleixner /* Now have argSq etc with binary point at the left
148da957e11SThomas Gleixner .1xxxxxxxx */
149da957e11SThomas Gleixner
150da957e11SThomas Gleixner /* Do the basic fixed point polynomial evaluation */
151da957e11SThomas Gleixner accumulator.msw = accumulator.midw = accumulator.lsw = 0;
152da957e11SThomas Gleixner polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
153da957e11SThomas Gleixner oddplterms, HIPOWERop - 1);
154da957e11SThomas Gleixner mul64_Xsig(&accumulator, &XSIG_LL(argSq));
155da957e11SThomas Gleixner negate_Xsig(&accumulator);
1563d0d14f9SIngo Molnar polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
1573d0d14f9SIngo Molnar HIPOWERon - 1);
158da957e11SThomas Gleixner negate_Xsig(&accumulator);
159da957e11SThomas Gleixner add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
160da957e11SThomas Gleixner
161da957e11SThomas Gleixner mul64_Xsig(&accumulatore, &denomterm);
162da957e11SThomas Gleixner shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
163da957e11SThomas Gleixner accumulatore.msw |= 0x80000000;
164da957e11SThomas Gleixner
165da957e11SThomas Gleixner div_Xsig(&accumulator, &accumulatore, &accumulator);
166da957e11SThomas Gleixner
167da957e11SThomas Gleixner mul_Xsig_Xsig(&accumulator, &argSignif);
168da957e11SThomas Gleixner mul_Xsig_Xsig(&accumulator, &argSq);
169da957e11SThomas Gleixner
170da957e11SThomas Gleixner shr_Xsig(&accumulator, 3);
171da957e11SThomas Gleixner negate_Xsig(&accumulator);
172da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &argSignif);
173da957e11SThomas Gleixner
1743d0d14f9SIngo Molnar if (transformed) {
175da957e11SThomas Gleixner /* compute pi/4 - accumulator */
176da957e11SThomas Gleixner shr_Xsig(&accumulator, -1 - exponent);
177da957e11SThomas Gleixner negate_Xsig(&accumulator);
178da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif);
179da957e11SThomas Gleixner exponent = -1;
180da957e11SThomas Gleixner }
181da957e11SThomas Gleixner
1823d0d14f9SIngo Molnar if (inverted) {
183da957e11SThomas Gleixner /* compute pi/2 - accumulator */
184da957e11SThomas Gleixner shr_Xsig(&accumulator, -exponent);
185da957e11SThomas Gleixner negate_Xsig(&accumulator);
186da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif);
187da957e11SThomas Gleixner exponent = 0;
188da957e11SThomas Gleixner }
189da957e11SThomas Gleixner
1903d0d14f9SIngo Molnar if (sign1) {
191da957e11SThomas Gleixner /* compute pi - accumulator */
192da957e11SThomas Gleixner shr_Xsig(&accumulator, 1 - exponent);
193da957e11SThomas Gleixner negate_Xsig(&accumulator);
194da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif);
195da957e11SThomas Gleixner exponent = 1;
196da957e11SThomas Gleixner }
197da957e11SThomas Gleixner
198da957e11SThomas Gleixner exponent += round_Xsig(&accumulator);
199da957e11SThomas Gleixner
200da957e11SThomas Gleixner significand(st1_ptr) = XSIG_LL(accumulator);
201da957e11SThomas Gleixner setexponent16(st1_ptr, exponent);
202da957e11SThomas Gleixner
203da957e11SThomas Gleixner tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
204da957e11SThomas Gleixner FPU_settagi(1, tag);
205da957e11SThomas Gleixner
206da957e11SThomas Gleixner set_precision_flag_up(); /* We do not really know if up or down,
207da957e11SThomas Gleixner use this as the default. */
208da957e11SThomas Gleixner
209da957e11SThomas Gleixner }
210