1 // SPDX-License-Identifier: GPL-2.0 1
2 /*--------------------------------------------
3 | poly_tan.c
4 |
5 | Compute the tan of a FPU_REG, using a polyn
6 |
7 | Copyright (C) 1992,1993,1994,1997,1999
8 | W. Metzenthen, 22 Par
9 | Australia. E-mail
10 |
11 |
12 +--------------------------------------------
13
14 #include "exception.h"
15 #include "reg_constant.h"
16 #include "fpu_emu.h"
17 #include "fpu_system.h"
18 #include "control_w.h"
19 #include "poly.h"
20
21 #define HiPOWERop 3 /* odd poly, p
22 static const unsigned long long oddplterm[HiPO
23 0x0000000000000000LL,
24 0x0051a1cf08fca228LL,
25 0x0000000071284ff7LL
26 };
27
28 #define HiPOWERon 2 /* odd poly, n
29 static const unsigned long long oddnegterm[HiP
30 0x1291a9a184244e80LL,
31 0x0000583245819c21LL
32 };
33
34 #define HiPOWERep 2 /* even poly,
35 static const unsigned long long evenplterm[HiP
36 0x0e848884b539e888LL,
37 0x00003c7f18b887daLL
38 };
39
40 #define HiPOWERen 2 /* even poly,
41 static const unsigned long long evennegterm[Hi
42 0xf1f0200fd51569ccLL,
43 0x003afb46105c4432LL
44 };
45
46 static const unsigned long long twothirds = 0x
47
48 /*--- poly_tan() -----------------------------
49 |
50 +--------------------------------------------
51 void poly_tan(FPU_REG *st0_ptr)
52 {
53 long int exponent;
54 int invert;
55 Xsig argSq, argSqSq, accumulatoro, acc
56 argSignif, fix_up;
57 unsigned long adj;
58
59 exponent = exponent(st0_ptr);
60
61 #ifdef PARANOID
62 if (signnegative(st0_ptr)) { /* Can
63 arith_invalid(0);
64 return;
65 } /* Need a posi
66 #endif /* PARANOID */
67
68 /* Split the problem into two domains,
69 if ((exponent == 0)
70 || ((exponent == -1) && (st0_ptr->
71 /* The argument is greater tha
72 invert = 1;
73 accum.lsw = 0;
74 XSIG_LL(accum) = significand(s
75
76 if (exponent == 0) {
77 /* The argument is >=
78 /* Put the binary poin
79 XSIG_LL(accum) <<= 1;
80 }
81 /* pi/2 in hex is: 1.921fb5444
82 XSIG_LL(accum) = 0x921fb54442d
83 /* This is a special case whic
84 if (XSIG_LL(accum) == 0xffffff
85 FPU_settag0(TAG_Valid)
86 significand(st0_ptr) =
87 setexponent16(st0_ptr,
88 (0x41 +
89 return;
90 }
91
92 argSignif.lsw = accum.lsw;
93 XSIG_LL(argSignif) = XSIG_LL(a
94 exponent = -1 + norm_Xsig(&arg
95 } else {
96 invert = 0;
97 argSignif.lsw = 0;
98 XSIG_LL(accum) = XSIG_LL(argSi
99
100 if (exponent < -1) {
101 /* shift the argument
102 if (FPU_shrx(&XSIG_LL(
103 0x80000000U)
104 XSIG_LL(accum)
105 }
106 }
107
108 XSIG_LL(argSq) = XSIG_LL(accum);
109 argSq.lsw = accum.lsw;
110 mul_Xsig_Xsig(&argSq, &argSq);
111 XSIG_LL(argSqSq) = XSIG_LL(argSq);
112 argSqSq.lsw = argSq.lsw;
113 mul_Xsig_Xsig(&argSqSq, &argSqSq);
114
115 /* Compute the negative terms for the
116 accumulatoro.msw = accumulatoro.midw =
117 polynomial_Xsig(&accumulatoro, &XSIG_L
118 HiPOWERon - 1);
119 mul_Xsig_Xsig(&accumulatoro, &argSq);
120 negate_Xsig(&accumulatoro);
121 /* Add the positive terms */
122 polynomial_Xsig(&accumulatoro, &XSIG_L
123 HiPOWERop - 1);
124
125 /* Compute the positive terms for the
126 accumulatore.msw = accumulatore.midw =
127 polynomial_Xsig(&accumulatore, &XSIG_L
128 HiPOWERep - 1);
129 mul_Xsig_Xsig(&accumulatore, &argSq);
130 negate_Xsig(&accumulatore);
131 /* Add the negative terms */
132 polynomial_Xsig(&accumulatore, &XSIG_L
133 HiPOWERen - 1);
134 /* Multiply by arg^2 */
135 mul64_Xsig(&accumulatore, &XSIG_LL(arg
136 mul64_Xsig(&accumulatore, &XSIG_LL(arg
137 /* de-normalize and divide by 2 */
138 shr_Xsig(&accumulatore, -2 * (1 + expo
139 negate_Xsig(&accumulatore); /* Thi
140
141 /* Now find the ratio. */
142 if (accumulatore.msw == 0) {
143 /* accumulatoro must contain 1
144 really doesn't matter what
145 have negligible effect in l
146 */
147 XSIG_LL(accum) = 0x80000000000
148 accum.lsw = 0;
149 } else {
150 div_Xsig(&accumulatoro, &accum
151 }
152
153 /* Multiply by 1/3 * arg^3 */
154 mul64_Xsig(&accum, &XSIG_LL(argSignif)
155 mul64_Xsig(&accum, &XSIG_LL(argSignif)
156 mul64_Xsig(&accum, &XSIG_LL(argSignif)
157 mul64_Xsig(&accum, &twothirds);
158 shr_Xsig(&accum, -2 * (exponent + 1));
159
160 /* tan(arg) = arg + accum */
161 add_two_Xsig(&accum, &argSignif, &expo
162
163 if (invert) {
164 /* We now have the value of ta
165 approximation for pi/2
166 */
167 /* The next step is to fix the
168 error due to the approximat
169 */
170
171 /* This is (approx) delta, the
172 (see above). It has an expo
173 */
174 XSIG_LL(fix_up) = 0x898cc51701
175 fix_up.lsw = 0;
176
177 if (exponent == 0)
178 adj = 0xffffffff;
179
180 else if (exponent > -30) {
181 adj = accum.msw >> -(e
182 adj = mul_32_32(adj, a
183 } else
184 adj = 0;
185 adj = mul_32_32(0x898cc517, ad
186
187 fix_up.msw += adj;
188 if (!(fix_up.msw & 0x80000000)
189 /* Yes, we need to add
190 shr_Xsig(&fix_up, 1);
191 fix_up.msw |= 0x800000
192 shr_Xsig(&fix_up, 64 +
193 } else
194 shr_Xsig(&fix_up, 65 +
195
196 add_two_Xsig(&accum, &fix_up,
197
198 /* accum now contains tan(pi/2
199 Use tan(arg) = 1.0 / tan(pi
200 */
201 accumulatoro.lsw = accumulator
202 accumulatoro.msw = 0x80000000;
203 div_Xsig(&accumulatoro, &accum
204 exponent = -exponent - 1;
205 }
206
207 /* Transfer the result */
208 round_Xsig(&accum);
209 FPU_settag0(TAG_Valid);
210 significand(st0_ptr) = XSIG_LL(accum);
211 setexponent16(st0_ptr, exponent + EXTE
212
213 }
214
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.