]>
Commit | Line | Data |
---|---|---|
d5efd131 MF |
1 | .file "tancotl.s" |
2 | ||
3 | ||
4 | // Copyright (c) 2000 - 2004, Intel Corporation | |
5 | // All rights reserved. | |
6 | // | |
7 | // Contributed 2000 by the Intel Numerics Group, Intel Corporation | |
8 | // | |
9 | // Redistribution and use in source and binary forms, with or without | |
10 | // modification, are permitted provided that the following conditions are | |
11 | // met: | |
12 | // | |
13 | // * Redistributions of source code must retain the above copyright | |
14 | // notice, this list of conditions and the following disclaimer. | |
15 | // | |
16 | // * Redistributions in binary form must reproduce the above copyright | |
17 | // notice, this list of conditions and the following disclaimer in the | |
18 | // documentation and/or other materials provided with the distribution. | |
19 | // | |
20 | // * The name of Intel Corporation may not be used to endorse or promote | |
21 | // products derived from this software without specific prior written | |
22 | // permission. | |
23 | ||
0347518d MF |
24 | // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
25 | // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | |
d5efd131 | 26 | // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
0347518d | 27 | // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS |
d5efd131 | 28 | // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
0347518d MF |
29 | // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
30 | // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | |
31 | // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY | |
d5efd131 | 32 | // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING |
0347518d MF |
33 | // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
34 | // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
35 | // | |
d5efd131 | 36 | // Intel Corporation is the author of this code, and requests that all |
0347518d | 37 | // problem reports or change requests be submitted to it directly at |
d5efd131 MF |
38 | // http://www.intel.com/software/products/opensource/libraries/num.htm. |
39 | // | |
40 | //********************************************************************* | |
41 | // | |
0347518d | 42 | // History: |
d5efd131 MF |
43 | // |
44 | // 02/02/00 (hand-optimized) | |
45 | // 04/04/00 Unwind support added | |
46 | // 12/28/00 Fixed false invalid flags | |
47 | // 02/06/02 Improved speed | |
48 | // 05/07/02 Changed interface to __libm_pi_by_2_reduce | |
49 | // 05/30/02 Added cotl | |
50 | // 02/10/03 Reordered header: .section, .global, .proc, .align; | |
51 | // used data8 for long double table values | |
52 | // 05/15/03 Reformatted data tables | |
53 | // 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader | |
54 | // | |
55 | //********************************************************************* | |
56 | // | |
57 | // Functions: tanl(x) = tangent(x), for double-extended precision x values | |
58 | // cotl(x) = cotangent(x), for double-extended precision x values | |
59 | // | |
60 | //********************************************************************* | |
61 | // | |
62 | // Resources Used: | |
63 | // | |
64 | // Floating-Point Registers: f8 (Input and Return Value) | |
65 | // f9-f15 | |
66 | // f32-f121 | |
67 | // | |
68 | // General Purpose Registers: | |
69 | // r32-r70 | |
70 | // | |
71 | // Predicate Registers: p6-p15 | |
72 | // | |
73 | //********************************************************************* | |
74 | // | |
75 | // IEEE Special Conditions for tanl: | |
76 | // | |
77 | // Denormal fault raised on denormal inputs | |
78 | // Overflow exceptions do not occur | |
79 | // Underflow exceptions raised when appropriate for tan | |
80 | // (No specialized error handling for this routine) | |
81 | // Inexact raised when appropriate by algorithm | |
82 | // | |
83 | // tanl(SNaN) = QNaN | |
84 | // tanl(QNaN) = QNaN | |
85 | // tanl(inf) = QNaN | |
86 | // tanl(+/-0) = +/-0 | |
87 | // | |
88 | //********************************************************************* | |
89 | // | |
90 | // IEEE Special Conditions for cotl: | |
91 | // | |
92 | // Denormal fault raised on denormal inputs | |
93 | // Overflow exceptions occur at zero and near zero | |
94 | // Underflow exceptions do not occur | |
95 | // Inexact raised when appropriate by algorithm | |
96 | // | |
97 | // cotl(SNaN) = QNaN | |
98 | // cotl(QNaN) = QNaN | |
99 | // cotl(inf) = QNaN | |
100 | // cotl(+/-0) = +/-Inf and error handling is called | |
101 | // | |
102 | //********************************************************************* | |
103 | // | |
104 | // Below are mathematical and algorithmic descriptions for tanl. | |
105 | // For cotl we use next identity cot(x) = -tan(x + Pi/2). | |
106 | // So, to compute cot(x) we just need to increment N (N = N + 1) | |
107 | // and invert sign of the computed result. | |
108 | // | |
109 | //********************************************************************* | |
110 | // | |
111 | // Mathematical Description | |
112 | // | |
113 | // We consider the computation of FPTANL of Arg. Now, given | |
114 | // | |
115 | // Arg = N pi/2 + alpha, |alpha| <= pi/4, | |
116 | // | |
117 | // basic mathematical relationship shows that | |
118 | // | |
119 | // tan( Arg ) = tan( alpha ) if N is even; | |
120 | // = -cot( alpha ) otherwise. | |
121 | // | |
122 | // The value of alpha is obtained by argument reduction and | |
123 | // represented by two working precision numbers r and c where | |
124 | // | |
125 | // alpha = r + c accurately. | |
126 | // | |
127 | // The reduction method is described in a previous write up. | |
128 | // The argument reduction scheme identifies 4 cases. For Cases 2 | |
129 | // and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be | |
130 | // computed very easily by 2 or 3 terms of the Taylor series | |
131 | // expansion as follows: | |
132 | // | |
133 | // Case 2: | |
134 | // ------- | |
135 | // | |
136 | // tan(r + c) = r + c + r^3/3 ...accurately | |
137 | // -cot(r + c) = -1/(r+c) + r/3 ...accurately | |
138 | // | |
139 | // Case 4: | |
140 | // ------- | |
141 | // | |
142 | // tan(r + c) = r + c + r^3/3 + 2r^5/15 ...accurately | |
143 | // -cot(r + c) = -1/(r+c) + r/3 + r^3/45 ...accurately | |
144 | // | |
145 | // | |
146 | // The only cases left are Cases 1 and 3 of the argument reduction | |
147 | // procedure. These two cases will be merged since after the | |
148 | // argument is reduced in either cases, we have the reduced argument | |
149 | // represented as r + c and that the magnitude |r + c| is not small | |
150 | // enough to allow the usage of a very short approximation. | |
151 | // | |
152 | // The greatest challenge of this task is that the second terms of | |
153 | // the Taylor series for tan(r) and -cot(r) | |
154 | // | |
155 | // r + r^3/3 + 2 r^5/15 + ... | |
156 | // | |
157 | // and | |
158 | // | |
159 | // -1/r + r/3 + r^3/45 + ... | |
160 | // | |
161 | // are not very small when |r| is close to pi/4 and the rounding | |
162 | // errors will be a concern if simple polynomial accumulation is | |
163 | // used. When |r| < 2^(-2), however, the second terms will be small | |
164 | // enough (5 bits or so of right shift) that a normal Horner | |
165 | // recurrence suffices. Hence there are two cases that we consider | |
166 | // in the accurate computation of tan(r) and cot(r), |r| <= pi/4. | |
167 | // | |
168 | // Case small_r: |r| < 2^(-2) | |
169 | // -------------------------- | |
170 | // | |
171 | // Since Arg = N pi/4 + r + c accurately, we have | |
172 | // | |
173 | // tan(Arg) = tan(r+c) for N even, | |
174 | // = -cot(r+c) otherwise. | |
175 | // | |
176 | // Here for this case, both tan(r) and -cot(r) can be approximated | |
177 | // by simple polynomials: | |
178 | // | |
179 | // tan(r) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19 | |
180 | // -cot(r) = -1/r + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13 | |
181 | // | |
182 | // accurately. Since |r| is relatively small, tan(r+c) and | |
183 | // -cot(r+c) can be accurately approximated by replacing r with | |
184 | // r+c only in the first two terms of the corresponding polynomials. | |
185 | // | |
186 | // Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to | |
187 | // almost 64 sig. bits, thus | |
188 | // | |
189 | // P1_1 (r+c)^3 = P1_1 r^3 + c * r^2 accurately. | |
190 | // | |
191 | // Hence, | |
192 | // | |
193 | // tan(r+c) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19 | |
194 | // + c*(1 + r^2) | |
195 | // | |
196 | // -cot(r+c) = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13 | |
197 | // + Q1_1*c | |
198 | // | |
199 | // | |
200 | // Case normal_r: 2^(-2) <= |r| <= pi/4 | |
201 | // ------------------------------------ | |
202 | // | |
203 | // This case is more likely than the previous one if one considers | |
204 | // r to be uniformly distributed in [-pi/4 pi/4]. | |
205 | // | |
206 | // The required calculation is either | |
207 | // | |
208 | // tan(r + c) = tan(r) + correction, or | |
209 | // -cot(r + c) = -cot(r) + correction. | |
210 | // | |
211 | // Specifically, | |
212 | // | |
213 | // tan(r + c) = tan(r) + c tan'(r) + O(c^2) | |
214 | // = tan(r) + c sec^2(r) + O(c^2) | |
215 | // = tan(r) + c SEC_sq ...accurately | |
216 | // as long as SEC_sq approximates sec^2(r) | |
217 | // to, say, 5 bits or so. | |
218 | // | |
219 | // Similarly, | |
220 | // | |
221 | // -cot(r + c) = -cot(r) - c cot'(r) + O(c^2) | |
222 | // = -cot(r) + c csc^2(r) + O(c^2) | |
223 | // = -cot(r) + c CSC_sq ...accurately | |
224 | // as long as CSC_sq approximates csc^2(r) | |
225 | // to, say, 5 bits or so. | |
226 | // | |
227 | // We therefore concentrate on accurately calculating tan(r) and | |
228 | // cot(r) for a working-precision number r, |r| <= pi/4 to within | |
229 | // 0.1% or so. | |
230 | // | |
231 | // We will employ a table-driven approach. Let | |
232 | // | |
233 | // r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63 | |
234 | // = sgn_r * ( B + x ) | |
235 | // | |
236 | // where | |
237 | // | |
238 | // B = 2^k * 1.b_1 b_2 ... b_5 1 | |
239 | // x = |r| - B | |
240 | // | |
241 | // Now, | |
242 | // tan(B) + tan(x) | |
243 | // tan( B + x ) = ------------------------ | |
244 | // 1 - tan(B)*tan(x) | |
245 | // | |
246 | // / \ | |
247 | // | tan(B) + tan(x) | | |
248 | ||
249 | // = tan(B) + | ------------------------ - tan(B) | | |
250 | // | 1 - tan(B)*tan(x) | | |
251 | // \ / | |
252 | // | |
253 | // sec^2(B) * tan(x) | |
254 | // = tan(B) + ------------------------ | |
255 | // 1 - tan(B)*tan(x) | |
256 | // | |
257 | // (1/[sin(B)*cos(B)]) * tan(x) | |
258 | // = tan(B) + -------------------------------- | |
259 | // cot(B) - tan(x) | |
260 | // | |
261 | // | |
262 | // Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are | |
263 | // calculated beforehand and stored in a table. Since | |
264 | // | |
265 | // |x| <= 2^k * 2^(-6) <= 2^(-7) (because k = -1, -2) | |
266 | // | |
267 | // a very short polynomial will be sufficient to approximate tan(x) | |
268 | // accurately. The details involved in computing the last expression | |
269 | // will be given in the next section on algorithm description. | |
270 | // | |
271 | // | |
272 | // Now, we turn to the case where cot( B + x ) is needed. | |
273 | // | |
274 | // | |
275 | // 1 - tan(B)*tan(x) | |
276 | // cot( B + x ) = ------------------------ | |
277 | // tan(B) + tan(x) | |
278 | // | |
279 | // / \ | |
280 | // | 1 - tan(B)*tan(x) | | |
281 | ||
282 | // = cot(B) + | ----------------------- - cot(B) | | |
283 | // | tan(B) + tan(x) | | |
284 | // \ / | |
285 | // | |
286 | // [tan(B) + cot(B)] * tan(x) | |
287 | // = cot(B) - ---------------------------- | |
288 | // tan(B) + tan(x) | |
289 | // | |
290 | // (1/[sin(B)*cos(B)]) * tan(x) | |
291 | // = cot(B) - -------------------------------- | |
292 | // tan(B) + tan(x) | |
293 | // | |
294 | // | |
295 | // Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that | |
296 | // are needed are the same set of values needed in the previous | |
297 | // case. | |
298 | // | |
299 | // Finally, we can put all the ingredients together as follows: | |
300 | // | |
301 | // Arg = N * pi/2 + r + c ...accurately | |
302 | // | |
303 | // tan(Arg) = tan(r) + correction if N is even; | |
304 | // = -cot(r) + correction otherwise. | |
305 | // | |
306 | // For Cases 2 and 4, | |
307 | // | |
308 | // Case 2: | |
309 | // tan(Arg) = tan(r + c) = r + c + r^3/3 N even | |
310 | // = -cot(r + c) = -1/(r+c) + r/3 N odd | |
311 | // Case 4: | |
312 | // tan(Arg) = tan(r + c) = r + c + r^3/3 + 2r^5/15 N even | |
313 | // = -cot(r + c) = -1/(r+c) + r/3 + r^3/45 N odd | |
314 | // | |
315 | // | |
316 | // For Cases 1 and 3, | |
317 | // | |
318 | // Case small_r: |r| < 2^(-2) | |
319 | // | |
320 | // tan(Arg) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19 | |
321 | // + c*(1 + r^2) N even | |
322 | // | |
323 | // = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13 | |
324 | // + Q1_1*c N odd | |
325 | // | |
326 | // Case normal_r: 2^(-2) <= |r| <= pi/4 | |
327 | // | |
328 | // tan(Arg) = tan(r) + c * sec^2(r) N even | |
329 | // = -cot(r) + c * csc^2(r) otherwise | |
330 | // | |
331 | // For N even, | |
332 | // | |
333 | // tan(Arg) = tan(r) + c*sec^2(r) | |
334 | // = tan( sgn_r * (B+x) ) + c * sec^2(|r|) | |
335 | // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(|r|) ) | |
336 | // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(B) ) | |
337 | // | |
338 | // since B approximates |r| to 2^(-6) in relative accuracy. | |
339 | // | |
340 | // / (1/[sin(B)*cos(B)]) * tan(x) | |
341 | // tan(Arg) = sgn_r * | tan(B) + -------------------------------- | |
342 | // \ cot(B) - tan(x) | |
343 | // \ | |
344 | // + CORR | | |
345 | ||
346 | // / | |
347 | // where | |
348 | // | |
349 | // CORR = sgn_r*c*tan(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)). | |
350 | // | |
351 | // For N odd, | |
352 | // | |
353 | // tan(Arg) = -cot(r) + c*csc^2(r) | |
354 | // = -cot( sgn_r * (B+x) ) + c * csc^2(|r|) | |
355 | // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(|r|) ) | |
356 | // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(B) ) | |
357 | // | |
358 | // since B approximates |r| to 2^(-6) in relative accuracy. | |
359 | // | |
360 | // / (1/[sin(B)*cos(B)]) * tan(x) | |
361 | // tan(Arg) = sgn_r * | -cot(B) + -------------------------------- | |
362 | // \ tan(B) + tan(x) | |
363 | // \ | |
364 | // + CORR | | |
365 | ||
366 | // / | |
367 | // where | |
368 | // | |
369 | // CORR = sgn_r*c*cot(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)). | |
370 | // | |
371 | // | |
372 | // The actual algorithm prescribes how all the mathematical formulas | |
373 | // are calculated. | |
374 | // | |
375 | // | |
376 | // 2. Algorithmic Description | |
377 | // ========================== | |
378 | // | |
379 | // 2.1 Computation for Cases 2 and 4. | |
380 | // ---------------------------------- | |
381 | // | |
382 | // For Case 2, we use two-term polynomials. | |
383 | // | |
384 | // For N even, | |
385 | // | |
386 | // rsq := r * r | |
387 | // Poly := c + r * rsq * P1_1 | |
388 | // Result := r + Poly ...in user-defined rounding | |
389 | // | |
390 | // For N odd, | |
391 | // S_hi := -frcpa(r) ...8 bits | |
392 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits | |
393 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits | |
394 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits | |
395 | // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c ) | |
396 | // ...S_hi + S_lo is -1/(r+c) to extra precision | |
397 | // S_lo := S_lo + Q1_1*r | |
398 | // | |
399 | // Result := S_hi + S_lo ...in user-defined rounding | |
400 | // | |
401 | // For Case 4, we use three-term polynomials | |
402 | // | |
403 | // For N even, | |
404 | // | |
405 | // rsq := r * r | |
406 | // Poly := c + r * rsq * (P1_1 + rsq * P1_2) | |
407 | // Result := r + Poly ...in user-defined rounding | |
408 | // | |
409 | // For N odd, | |
410 | // S_hi := -frcpa(r) ...8 bits | |
411 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits | |
412 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits | |
413 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits | |
414 | // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c ) | |
415 | // ...S_hi + S_lo is -1/(r+c) to extra precision | |
416 | // rsq := r * r | |
417 | // P := Q1_1 + rsq*Q1_2 | |
418 | // S_lo := S_lo + r*P | |
419 | // | |
420 | // Result := S_hi + S_lo ...in user-defined rounding | |
421 | // | |
422 | // | |
423 | // Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are | |
424 | // the same as those used in the small_r case of Cases 1 and 3 | |
425 | // below. | |
426 | // | |
427 | // | |
428 | // 2.2 Computation for Cases 1 and 3. | |
429 | // ---------------------------------- | |
430 | // This is further divided into the case of small_r, | |
431 | // where |r| < 2^(-2), and the case of normal_r, where |r| lies between | |
432 | // 2^(-2) and pi/4. | |
433 | // | |
434 | // Algorithm for the case of small_r | |
435 | // --------------------------------- | |
436 | // | |
437 | // For N even, | |
438 | // rsq := r * r | |
439 | // Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3)) | |
440 | // r_to_the_8 := rsq * rsq | |
441 | // r_to_the_8 := r_to_the_8 * r_to_the_8 | |
442 | // Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9)) | |
443 | // CORR := c * ( 1 + rsq ) | |
444 | // Poly := Poly1 + r_to_the_8*Poly2 | |
445 | // Poly := r*Poly + CORR | |
446 | // Result := r + Poly ...in user-defined rounding | |
447 | // ...note that Poly1 and r_to_the_8 can be computed in parallel | |
448 | // ...with Poly2 (Poly1 is intentionally set to be much | |
449 | // ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden) | |
450 | // | |
451 | // For N odd, | |
452 | // S_hi := -frcpa(r) ...8 bits | |
453 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits | |
454 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits | |
455 | // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits | |
456 | // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c ) | |
457 | // ...S_hi + S_lo is -1/(r+c) to extra precision | |
458 | // S_lo := S_lo + Q1_1*c | |
459 | // | |
460 | // ...S_hi and S_lo are computed in parallel with | |
461 | // ...the following | |
462 | // rsq := r*r | |
463 | // P := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7)) | |
464 | // | |
465 | // Poly := r*P + S_lo | |
466 | // Result := S_hi + Poly ...in user-defined rounding | |
467 | // | |
468 | // | |
469 | // Algorithm for the case of normal_r | |
470 | // ---------------------------------- | |
471 | // | |
472 | // Here, we first consider the computation of tan( r + c ). As | |
473 | // presented in the previous section, | |
474 | // | |
475 | // tan( r + c ) = tan(r) + c * sec^2(r) | |
476 | // = sgn_r * [ tan(B+x) + CORR ] | |
477 | // CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)] | |
478 | // | |
479 | // because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits. | |
480 | // | |
481 | // tan( r + c ) = | |
482 | // / (1/[sin(B)*cos(B)]) * tan(x) | |
483 | // sgn_r * | tan(B) + -------------------------------- + | |
484 | // \ cot(B) - tan(x) | |
485 | // \ | |
486 | // CORR | | |
487 | ||
488 | // / | |
489 | // | |
490 | // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are | |
491 | // calculated beforehand and stored in a table. Specifically, | |
492 | // the table values are | |
493 | // | |
494 | // tan(B) as T_hi + T_lo; | |
495 | // cot(B) as C_hi + C_lo; | |
496 | // 1/[sin(B)*cos(B)] as SC_inv | |
497 | // | |
498 | // T_hi, C_hi are in double-precision memory format; | |
499 | // T_lo, C_lo are in single-precision memory format; | |
500 | // SC_inv is in extended-precision memory format. | |
501 | // | |
502 | // The value of tan(x) will be approximated by a short polynomial of | |
503 | // the form | |
504 | // | |
505 | // tan(x) as x + x * P, where | |
506 | // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3)) | |
507 | // | |
508 | // Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x) | |
509 | // to a relative accuracy better than 2^(-20). Thus, a good | |
510 | // initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative | |
511 | // division is: | |
512 | // | |
513 | // 1/(cot(B) - tan(x)) is approximately | |
514 | // 1/(cot(B) - x) is | |
515 | // tan(B)/(1 - x*tan(B)) is approximately | |
516 | // T_hi / ( 1 - T_hi * x ) is approximately | |
517 | // | |
518 | // T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ] | |
519 | // | |
520 | // The calculation of tan(r+c) therefore proceed as follows: | |
521 | // | |
522 | // Tx := T_hi * x | |
523 | // xsq := x * x | |
524 | // | |
525 | // V_hi := T_hi*(1 + Tx*(1 + Tx)) | |
526 | // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3)) | |
527 | // ...V_hi serves as an initial guess of 1/(cot(B) - tan(x)) | |
528 | // ...good to about 20 bits of accuracy | |
529 | // | |
530 | // tanx := x + x*P | |
531 | // D := C_hi - tanx | |
532 | // ...D is a double precision denominator: cot(B) - tan(x) | |
533 | // | |
534 | // V_hi := V_hi + V_hi*(1 - V_hi*D) | |
535 | // ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits | |
536 | // | |
537 | // V_lo := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ] | |
538 | // - V_hi*C_lo ) ...observe all order | |
539 | // ...V_hi + V_lo approximates 1/(cot(B) - tan(x)) | |
540 | // ...to extra accuracy | |
541 | // | |
542 | // ... SC_inv(B) * (x + x*P) | |
543 | // ... tan(B) + ------------------------- + CORR | |
544 | // ... cot(B) - (x + x*P) | |
545 | // ... | |
546 | // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR | |
547 | // ... | |
548 | // | |
549 | // Sx := SC_inv * x | |
550 | // CORR := sgn_r * c * SC_inv * T_hi | |
551 | // | |
552 | // ...put the ingredients together to compute | |
553 | // ... SC_inv(B) * (x + x*P) | |
554 | // ... tan(B) + ------------------------- + CORR | |
555 | // ... cot(B) - (x + x*P) | |
556 | // ... | |
557 | // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR | |
558 | // ... | |
559 | // ... = T_hi + T_lo + CORR + | |
560 | // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo) | |
561 | // | |
562 | // CORR := CORR + T_lo | |
563 | // tail := V_lo + P*(V_hi + V_lo) | |
564 | // tail := Sx * tail + CORR | |
565 | // tail := Sx * V_hi + tail | |
566 | // T_hi := sgn_r * T_hi | |
567 | // | |
568 | // ...T_hi + sgn_r*tail now approximate | |
569 | // ...sgn_r*(tan(B+x) + CORR) accurately | |
570 | // | |
571 | // Result := T_hi + sgn_r*tail ...in user-defined | |
572 | // ...rounding control | |
573 | // ...It is crucial that independent paths be fully | |
574 | // ...exploited for performance's sake. | |
575 | // | |
576 | // | |
577 | // Next, we consider the computation of -cot( r + c ). As | |
578 | // presented in the previous section, | |
579 | // | |
580 | // -cot( r + c ) = -cot(r) + c * csc^2(r) | |
581 | // = sgn_r * [ -cot(B+x) + CORR ] | |
582 | // CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)] | |
583 | // | |
584 | // because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits. | |
585 | // | |
586 | // -cot( r + c ) = | |
587 | // / (1/[sin(B)*cos(B)]) * tan(x) | |
588 | // sgn_r * | -cot(B) + -------------------------------- + | |
589 | // \ tan(B) + tan(x) | |
590 | // \ | |
591 | // CORR | | |
592 | ||
593 | // / | |
594 | // | |
595 | // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are | |
596 | // calculated beforehand and stored in a table. Specifically, | |
597 | // the table values are | |
598 | // | |
599 | // tan(B) as T_hi + T_lo; | |
600 | // cot(B) as C_hi + C_lo; | |
601 | // 1/[sin(B)*cos(B)] as SC_inv | |
602 | // | |
603 | // T_hi, C_hi are in double-precision memory format; | |
604 | // T_lo, C_lo are in single-precision memory format; | |
605 | // SC_inv is in extended-precision memory format. | |
606 | // | |
607 | // The value of tan(x) will be approximated by a short polynomial of | |
608 | // the form | |
609 | // | |
610 | // tan(x) as x + x * P, where | |
611 | // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3)) | |
612 | // | |
613 | // Because |x| <= 2^(-7), tan(B) + x approximates tan(B) + tan(x) | |
614 | // to a relative accuracy better than 2^(-18). Thus, a good | |
615 | // initial guess of 1/( tan(B) + tan(x) ) to initiate the iterative | |
616 | // division is: | |
617 | // | |
618 | // 1/(tan(B) + tan(x)) is approximately | |
619 | // 1/(tan(B) + x) is | |
620 | // cot(B)/(1 + x*cot(B)) is approximately | |
621 | // C_hi / ( 1 + C_hi * x ) is approximately | |
622 | // | |
623 | // C_hi * [ 1 - (C_hi * x) + (C_hi * x)^2 ] | |
624 | // | |
625 | // The calculation of -cot(r+c) therefore proceed as follows: | |
626 | // | |
627 | // Cx := C_hi * x | |
628 | // xsq := x * x | |
629 | // | |
630 | // V_hi := C_hi*(1 - Cx*(1 - Cx)) | |
631 | // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3)) | |
632 | // ...V_hi serves as an initial guess of 1/(tan(B) + tan(x)) | |
633 | // ...good to about 18 bits of accuracy | |
634 | // | |
635 | // tanx := x + x*P | |
636 | // D := T_hi + tanx | |
637 | // ...D is a double precision denominator: tan(B) + tan(x) | |
638 | // | |
639 | // V_hi := V_hi + V_hi*(1 - V_hi*D) | |
640 | // ....V_hi approximates 1/(tan(B)+tan(x)) to 40 bits | |
641 | // | |
642 | // V_lo := V_hi * ( [ (1 - V_hi*T_hi) - V_hi*tanx ] | |
643 | // - V_hi*T_lo ) ...observe all order | |
644 | // ...V_hi + V_lo approximates 1/(tan(B) + tan(x)) | |
645 | // ...to extra accuracy | |
646 | // | |
647 | // ... SC_inv(B) * (x + x*P) | |
648 | // ... -cot(B) + ------------------------- + CORR | |
649 | // ... tan(B) + (x + x*P) | |
650 | // ... | |
651 | // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR | |
652 | // ... | |
653 | // | |
654 | // Sx := SC_inv * x | |
655 | // CORR := sgn_r * c * SC_inv * C_hi | |
656 | // | |
657 | // ...put the ingredients together to compute | |
658 | // ... SC_inv(B) * (x + x*P) | |
659 | // ... -cot(B) + ------------------------- + CORR | |
660 | // ... tan(B) + (x + x*P) | |
661 | // ... | |
662 | // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR | |
663 | // ... | |
664 | // ... =-C_hi - C_lo + CORR + | |
665 | // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo) | |
666 | // | |
667 | // CORR := CORR - C_lo | |
668 | // tail := V_lo + P*(V_hi + V_lo) | |
669 | // tail := Sx * tail + CORR | |
670 | // tail := Sx * V_hi + tail | |
671 | // C_hi := -sgn_r * C_hi | |
672 | // | |
673 | // ...C_hi + sgn_r*tail now approximates | |
674 | // ...sgn_r*(-cot(B+x) + CORR) accurately | |
675 | // | |
676 | // Result := C_hi + sgn_r*tail in user-defined rounding control | |
677 | // ...It is crucial that independent paths be fully | |
678 | // ...exploited for performance's sake. | |
679 | // | |
680 | // 3. Implementation Notes | |
681 | // ======================= | |
682 | // | |
683 | // Table entries T_hi, T_lo; C_hi, C_lo; SC_inv | |
684 | // | |
685 | // Recall that 2^(-2) <= |r| <= pi/4; | |
686 | // | |
687 | // r = sgn_r * 2^k * 1.b_1 b_2 ... b_63 | |
688 | // | |
689 | // and | |
690 | // | |
691 | // B = 2^k * 1.b_1 b_2 b_3 b_4 b_5 1 | |
692 | // | |
693 | // Thus, for k = -2, possible values of B are | |
694 | // | |
695 | // B = 2^(-2) * ( 1 + index/32 + 1/64 ), | |
696 | // index ranges from 0 to 31 | |
697 | // | |
698 | // For k = -1, however, since |r| <= pi/4 = 0.78... | |
699 | // possible values of B are | |
700 | // | |
701 | // B = 2^(-1) * ( 1 + index/32 + 1/64 ) | |
702 | // index ranges from 0 to 19. | |
703 | // | |
704 | // | |
705 | ||
706 | RODATA | |
707 | .align 16 | |
708 | ||
709 | LOCAL_OBJECT_START(TANL_BASE_CONSTANTS) | |
710 | ||
711 | tanl_table_1: | |
712 | data8 0xA2F9836E4E44152A, 0x00003FFE // two_by_pi | |
713 | data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0 | |
714 | data8 0xC90FDAA22168C235, 0x00003FFF // P_1 | |
715 | data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2 | |
716 | data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3 | |
717 | LOCAL_OBJECT_END(TANL_BASE_CONSTANTS) | |
718 | ||
719 | LOCAL_OBJECT_START(tanl_table_2) | |
720 | data8 0xC90FDAA22168C234, 0x00003FFE // PI_BY_4 | |
721 | data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0 | |
722 | data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1 | |
723 | data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2 | |
724 | data4 0x3E800000 // two**-2 | |
725 | data4 0xBE800000 // -two**-2 | |
726 | data4 0x00000000 // pad | |
727 | data4 0x00000000 // pad | |
728 | LOCAL_OBJECT_END(tanl_table_2) | |
729 | ||
730 | LOCAL_OBJECT_START(tanl_table_p1) | |
731 | data8 0xAAAAAAAAAAAAAABD, 0x00003FFD // P1_1 | |
732 | data8 0x8888888888882E6A, 0x00003FFC // P1_2 | |
733 | data8 0xDD0DD0DD0F0177B6, 0x00003FFA // P1_3 | |
734 | data8 0xB327A440646B8C6D, 0x00003FF9 // P1_4 | |
735 | data8 0x91371B251D5F7D20, 0x00003FF8 // P1_5 | |
736 | data8 0xEB69A5F161C67914, 0x00003FF6 // P1_6 | |
737 | data8 0xBEDD37BE019318D2, 0x00003FF5 // P1_7 | |
738 | data8 0x9979B1463C794015, 0x00003FF4 // P1_8 | |
739 | data8 0x8EBD21A38C6EB58A, 0x00003FF3 // P1_9 | |
740 | LOCAL_OBJECT_END(tanl_table_p1) | |
741 | ||
742 | LOCAL_OBJECT_START(tanl_table_q1) | |
743 | data8 0xAAAAAAAAAAAAAAB4, 0x00003FFD // Q1_1 | |
744 | data8 0xB60B60B60B5FC93E, 0x00003FF9 // Q1_2 | |
745 | data8 0x8AB355E00C9BBFBF, 0x00003FF6 // Q1_3 | |
746 | data8 0xDDEBBC89CBEE3D4C, 0x00003FF2 // Q1_4 | |
747 | data8 0xB3548A685F80BBB6, 0x00003FEF // Q1_5 | |
748 | data8 0x913625604CED5BF1, 0x00003FEC // Q1_6 | |
749 | data8 0xF189D95A8EE92A83, 0x00003FE8 // Q1_7 | |
750 | LOCAL_OBJECT_END(tanl_table_q1) | |
751 | ||
752 | LOCAL_OBJECT_START(tanl_table_p2) | |
753 | data8 0xAAAAAAAAAAAB362F, 0x00003FFD // P2_1 | |
754 | data8 0x88888886E97A6097, 0x00003FFC // P2_2 | |
755 | data8 0xDD108EE025E716A1, 0x00003FFA // P2_3 | |
756 | LOCAL_OBJECT_END(tanl_table_p2) | |
757 | ||
758 | LOCAL_OBJECT_START(tanl_table_tm2) | |
759 | // | |
760 | // Entries T_hi double-precision memory format | |
761 | // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64) | |
762 | // Entries T_lo single-precision memory format | |
763 | // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64) | |
764 | // | |
765 | data8 0x3FD09BC362400794 | |
766 | data4 0x23A05C32, 0x00000000 | |
767 | data8 0x3FD124A9DFFBC074 | |
768 | data4 0x240078B2, 0x00000000 | |
769 | data8 0x3FD1AE235BD4920F | |
770 | data4 0x23826B8E, 0x00000000 | |
771 | data8 0x3FD2383515E2701D | |
772 | data4 0x22D31154, 0x00000000 | |
773 | data8 0x3FD2C2E463739C2D | |
774 | data4 0x2265C9E2, 0x00000000 | |
775 | data8 0x3FD34E36AFEEA48B | |
776 | data4 0x245C05EB, 0x00000000 | |
777 | data8 0x3FD3DA317DBB35D1 | |
778 | data4 0x24749F2D, 0x00000000 | |
779 | data8 0x3FD466DA67321619 | |
780 | data4 0x2462CECE, 0x00000000 | |
781 | data8 0x3FD4F4371F94A4D5 | |
782 | data4 0x246D0DF1, 0x00000000 | |
783 | data8 0x3FD5824D740C3E6D | |
784 | data4 0x240A85B5, 0x00000000 | |
785 | data8 0x3FD611234CB1E73D | |
786 | data4 0x23F96E33, 0x00000000 | |
787 | data8 0x3FD6A0BEAD9EA64B | |
788 | data4 0x247C5393, 0x00000000 | |
789 | data8 0x3FD73125B804FD01 | |
790 | data4 0x241F3B29, 0x00000000 | |
791 | data8 0x3FD7C25EAB53EE83 | |
792 | data4 0x2479989B, 0x00000000 | |
793 | data8 0x3FD8546FE6640EED | |
794 | data4 0x23B343BC, 0x00000000 | |
795 | data8 0x3FD8E75FE8AF1892 | |
796 | data4 0x241454D1, 0x00000000 | |
797 | data8 0x3FD97B3553928BDA | |
798 | data4 0x238613D9, 0x00000000 | |
799 | data8 0x3FDA0FF6EB9DE4DE | |
800 | data4 0x22859FA7, 0x00000000 | |
801 | data8 0x3FDAA5AB99ECF92D | |
802 | data4 0x237A6D06, 0x00000000 | |
803 | data8 0x3FDB3C5A6D8F1796 | |
804 | data4 0x23952F6C, 0x00000000 | |
805 | data8 0x3FDBD40A9CFB8BE4 | |
806 | data4 0x2280FC95, 0x00000000 | |
807 | data8 0x3FDC6CC387943100 | |
808 | data4 0x245D2EC0, 0x00000000 | |
809 | data8 0x3FDD068CB736C500 | |
810 | data4 0x23C4AD7D, 0x00000000 | |
811 | data8 0x3FDDA16DE1DDBC31 | |
812 | data4 0x23D076E6, 0x00000000 | |
813 | data8 0x3FDE3D6EEB515A93 | |
814 | data4 0x244809A6, 0x00000000 | |
815 | data8 0x3FDEDA97E6E9E5F1 | |
816 | data4 0x220856C8, 0x00000000 | |
817 | data8 0x3FDF78F11963CE69 | |
818 | data4 0x244BE993, 0x00000000 | |
819 | data8 0x3FE00C417D635BCE | |
820 | data4 0x23D21799, 0x00000000 | |
821 | data8 0x3FE05CAB1C302CD3 | |
822 | data4 0x248A1B1D, 0x00000000 | |
823 | data8 0x3FE0ADB9DB6A1FA0 | |
824 | data4 0x23D53E33, 0x00000000 | |
825 | data8 0x3FE0FF724A20BA81 | |
826 | data4 0x24DB9ED5, 0x00000000 | |
827 | data8 0x3FE151D9153FA6F5 | |
828 | data4 0x24E9E451, 0x00000000 | |
829 | LOCAL_OBJECT_END(tanl_table_tm2) | |
830 | ||
831 | LOCAL_OBJECT_START(tanl_table_tm1) | |
832 | // | |
833 | // Entries T_hi double-precision memory format | |
834 | // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64) | |
835 | // Entries T_lo single-precision memory format | |
836 | // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64) | |
837 | // | |
838 | data8 0x3FE1CEC4BA1BE39E | |
839 | data4 0x24B60F9E, 0x00000000 | |
840 | data8 0x3FE277E45ABD9B2D | |
841 | data4 0x248C2474, 0x00000000 | |
842 | data8 0x3FE324180272B110 | |
843 | data4 0x247B8311, 0x00000000 | |
844 | data8 0x3FE3D38B890E2DF0 | |
845 | data4 0x24C55751, 0x00000000 | |
846 | data8 0x3FE4866D46236871 | |
847 | data4 0x24E5BC34, 0x00000000 | |
848 | data8 0x3FE53CEE45E044B0 | |
849 | data4 0x24001BA4, 0x00000000 | |
850 | data8 0x3FE5F74282EC06E4 | |
851 | data4 0x24B973DC, 0x00000000 | |
852 | data8 0x3FE6B5A125DF43F9 | |
853 | data4 0x24895440, 0x00000000 | |
854 | data8 0x3FE77844CAFD348C | |
855 | data4 0x240021CA, 0x00000000 | |
856 | data8 0x3FE83F6BCEED6B92 | |
857 | data4 0x24C45372, 0x00000000 | |
858 | data8 0x3FE90B58A34F3665 | |
859 | data4 0x240DAD33, 0x00000000 | |
860 | data8 0x3FE9DC522C1E56B4 | |
861 | data4 0x24F846CE, 0x00000000 | |
862 | data8 0x3FEAB2A427041578 | |
863 | data4 0x2323FB6E, 0x00000000 | |
864 | data8 0x3FEB8E9F9DD8C373 | |
865 | data4 0x24B3090B, 0x00000000 | |
866 | data8 0x3FEC709B65C9AA7B | |
867 | data4 0x2449F611, 0x00000000 | |
868 | data8 0x3FED58F4ACCF8435 | |
869 | data4 0x23616A7E, 0x00000000 | |
870 | data8 0x3FEE480F97635082 | |
871 | data4 0x24C2FEAE, 0x00000000 | |
872 | data8 0x3FEF3E57F0ACC544 | |
873 | data4 0x242CE964, 0x00000000 | |
874 | data8 0x3FF01E20F7E06E4B | |
875 | data4 0x2480D3EE, 0x00000000 | |
876 | data8 0x3FF0A1258A798A69 | |
877 | data4 0x24DB8967, 0x00000000 | |
878 | LOCAL_OBJECT_END(tanl_table_tm1) | |
879 | ||
880 | LOCAL_OBJECT_START(tanl_table_cm2) | |
881 | // | |
882 | // Entries C_hi double-precision memory format | |
883 | // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64) | |
884 | // Entries C_lo single-precision memory format | |
885 | // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64) | |
886 | // | |
887 | data8 0x400ED3E2E63EFBD0 | |
888 | data4 0x259D94D4, 0x00000000 | |
889 | data8 0x400DDDB4C515DAB5 | |
890 | data4 0x245F0537, 0x00000000 | |
891 | data8 0x400CF57ABE19A79F | |
892 | data4 0x25D4EA9F, 0x00000000 | |
893 | data8 0x400C1A06D15298ED | |
894 | data4 0x24AE40A0, 0x00000000 | |
895 | data8 0x400B4A4C164B2708 | |
896 | data4 0x25A5AAB6, 0x00000000 | |
897 | data8 0x400A855A5285B068 | |
898 | data4 0x25524F18, 0x00000000 | |
899 | data8 0x4009CA5A3FFA549F | |
900 | data4 0x24C999C0, 0x00000000 | |
901 | data8 0x4009188A646AF623 | |
902 | data4 0x254FD801, 0x00000000 | |
903 | data8 0x40086F3C6084D0E7 | |
904 | data4 0x2560F5FD, 0x00000000 | |
905 | data8 0x4007CDD2A29A76EE | |
906 | data4 0x255B9D19, 0x00000000 | |
907 | data8 0x400733BE6C8ECA95 | |
908 | data4 0x25CB021B, 0x00000000 | |
909 | data8 0x4006A07E1F8DDC52 | |
910 | data4 0x24AB4722, 0x00000000 | |
911 | data8 0x4006139BC298AD58 | |
912 | data4 0x252764E2, 0x00000000 | |
913 | data8 0x40058CABBAD7164B | |
914 | data4 0x24DAF5DB, 0x00000000 | |
915 | data8 0x40050B4BAE31A5D3 | |
916 | data4 0x25EA20F4, 0x00000000 | |
917 | data8 0x40048F2189F85A8A | |
918 | data4 0x2583A3E8, 0x00000000 | |
919 | data8 0x400417DAA862380D | |
920 | data4 0x25DCC4CC, 0x00000000 | |
921 | data8 0x4003A52B1088FCFE | |
922 | data4 0x2430A492, 0x00000000 | |
923 | data8 0x400336CCCD3527D5 | |
924 | data4 0x255F77CF, 0x00000000 | |
925 | data8 0x4002CC7F5760766D | |
926 | data4 0x25DA0BDA, 0x00000000 | |
927 | data8 0x4002660711CE02E3 | |
928 | data4 0x256FF4A2, 0x00000000 | |
929 | data8 0x4002032CD37BBE04 | |
930 | data4 0x25208AED, 0x00000000 | |
931 | data8 0x4001A3BD7F050775 | |
932 | data4 0x24B72DD6, 0x00000000 | |
933 | data8 0x40014789A554848A | |
934 | data4 0x24AB4DAA, 0x00000000 | |
935 | data8 0x4000EE65323E81B7 | |
936 | data4 0x2584C440, 0x00000000 | |
937 | data8 0x4000982721CF1293 | |
938 | data4 0x25C9428D, 0x00000000 | |
939 | data8 0x400044A93D415EEB | |
940 | data4 0x25DC8482, 0x00000000 | |
941 | data8 0x3FFFE78FBD72C577 | |
942 | data4 0x257F5070, 0x00000000 | |
943 | data8 0x3FFF4AC375EFD28E | |
944 | data4 0x23EBBF7A, 0x00000000 | |
945 | data8 0x3FFEB2AF60B52DDE | |
946 | data4 0x22EECA07, 0x00000000 | |
947 | data8 0x3FFE1F1935204180 | |
948 | data4 0x24191079, 0x00000000 | |
949 | data8 0x3FFD8FCA54F7E60A | |
950 | data4 0x248D3058, 0x00000000 | |
951 | LOCAL_OBJECT_END(tanl_table_cm2) | |
952 | ||
953 | LOCAL_OBJECT_START(tanl_table_cm1) | |
954 | // | |
955 | // Entries C_hi double-precision memory format | |
956 | // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64) | |
957 | // Entries C_lo single-precision memory format | |
958 | // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64) | |
959 | // | |
960 | data8 0x3FFCC06A79F6FADE | |
961 | data4 0x239C7886, 0x00000000 | |
962 | data8 0x3FFBB91F891662A6 | |
963 | data4 0x250BD191, 0x00000000 | |
964 | data8 0x3FFABFB6529F155D | |
965 | data4 0x256CC3E6, 0x00000000 | |
966 | data8 0x3FF9D3002E964AE9 | |
967 | data4 0x250843E3, 0x00000000 | |
968 | data8 0x3FF8F1EF89DCB383 | |
969 | data4 0x2277C87E, 0x00000000 | |
970 | data8 0x3FF81B937C87DBD6 | |
971 | data4 0x256DA6CF, 0x00000000 | |
972 | data8 0x3FF74F141042EDE4 | |
973 | data4 0x2573D28A, 0x00000000 | |
974 | data8 0x3FF68BAF1784B360 | |
975 | data4 0x242E489A, 0x00000000 | |
976 | data8 0x3FF5D0B57C923C4C | |
977 | data4 0x2532D940, 0x00000000 | |
978 | data8 0x3FF51D88F418EF20 | |
979 | data4 0x253C7DD6, 0x00000000 | |
980 | data8 0x3FF4719A02F88DAE | |
981 | data4 0x23DB59BF, 0x00000000 | |
982 | data8 0x3FF3CC6649DA0788 | |
983 | data4 0x252B4756, 0x00000000 | |
984 | data8 0x3FF32D770B980DB8 | |
985 | data4 0x23FE585F, 0x00000000 | |
986 | data8 0x3FF2945FE56C987A | |
987 | data4 0x25378A63, 0x00000000 | |
988 | data8 0x3FF200BDB16523F6 | |
989 | data4 0x247BB2E0, 0x00000000 | |
990 | data8 0x3FF172358CE27778 | |
991 | data4 0x24446538, 0x00000000 | |
992 | data8 0x3FF0E873FDEFE692 | |
993 | data4 0x2514638F, 0x00000000 | |
994 | data8 0x3FF0632C33154062 | |
995 | data4 0x24A7FC27, 0x00000000 | |
996 | data8 0x3FEFC42EB3EF115F | |
997 | data4 0x248FD0FE, 0x00000000 | |
998 | data8 0x3FEEC9E8135D26F6 | |
999 | data4 0x2385C719, 0x00000000 | |
1000 | LOCAL_OBJECT_END(tanl_table_cm1) | |
1001 | ||
1002 | LOCAL_OBJECT_START(tanl_table_scim2) | |
1003 | // | |
1004 | // Entries SC_inv in Swapped IEEE format (extended) | |
1005 | // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64) | |
1006 | // | |
1007 | data8 0x839D6D4A1BF30C9E, 0x00004001 | |
1008 | data8 0x80092804554B0EB0, 0x00004001 | |
1009 | data8 0xF959F94CA1CF0DE9, 0x00004000 | |
1010 | data8 0xF3086BA077378677, 0x00004000 | |
1011 | data8 0xED154515CCD4723C, 0x00004000 | |
1012 | data8 0xE77909441C27CF25, 0x00004000 | |
1013 | data8 0xE22D037D8DDACB88, 0x00004000 | |
1014 | data8 0xDD2B2D8A89C73522, 0x00004000 | |
1015 | data8 0xD86E1A23BB2C1171, 0x00004000 | |
1016 | data8 0xD3F0E288DFF5E0F9, 0x00004000 | |
1017 | data8 0xCFAF16B1283BEBD5, 0x00004000 | |
1018 | data8 0xCBA4AFAA0D88DD53, 0x00004000 | |
1019 | data8 0xC7CE03CCCA67C43D, 0x00004000 | |
1020 | data8 0xC427BC820CA0DDB0, 0x00004000 | |
1021 | data8 0xC0AECD57F13D8CAB, 0x00004000 | |
1022 | data8 0xBD606C3871ECE6B1, 0x00004000 | |
1023 | data8 0xBA3A0A96A44C4929, 0x00004000 | |
1024 | data8 0xB7394F6FE5CCCEC1, 0x00004000 | |
1025 | data8 0xB45C12039637D8BC, 0x00004000 | |
1026 | data8 0xB1A0552892CB051B, 0x00004000 | |
1027 | data8 0xAF04432B6BA2FFD0, 0x00004000 | |
1028 | data8 0xAC862A237221235F, 0x00004000 | |
1029 | data8 0xAA2478AF5F00A9D1, 0x00004000 | |
1030 | data8 0xA7DDBB0C81E082BF, 0x00004000 | |
1031 | data8 0xA5B0987D45684FEE, 0x00004000 | |
1032 | data8 0xA39BD0F5627A8F53, 0x00004000 | |
1033 | data8 0xA19E3B036EC5C8B0, 0x00004000 | |
1034 | data8 0x9FB6C1F091CD7C66, 0x00004000 | |
1035 | data8 0x9DE464101FA3DF8A, 0x00004000 | |
1036 | data8 0x9C263139A8F6B888, 0x00004000 | |
1037 | data8 0x9A7B4968C27B0450, 0x00004000 | |
1038 | data8 0x98E2DB7E5EE614EE, 0x00004000 | |
1039 | LOCAL_OBJECT_END(tanl_table_scim2) | |
1040 | ||
1041 | LOCAL_OBJECT_START(tanl_table_scim1) | |
1042 | // | |
1043 | // Entries SC_inv in Swapped IEEE format (extended) | |
1044 | // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64) | |
1045 | // | |
1046 | data8 0x969F335C13B2B5BA, 0x00004000 | |
1047 | data8 0x93D446D9D4C0F548, 0x00004000 | |
1048 | data8 0x9147094F61B798AF, 0x00004000 | |
1049 | data8 0x8EF317CC758787AC, 0x00004000 | |
1050 | data8 0x8CD498B3B99EEFDB, 0x00004000 | |
1051 | data8 0x8AE82A7DDFF8BC37, 0x00004000 | |
1052 | data8 0x892AD546E3C55D42, 0x00004000 | |
1053 | data8 0x8799FEA9D15573C1, 0x00004000 | |
1054 | data8 0x86335F88435A4B4C, 0x00004000 | |
1055 | data8 0x84F4FB6E3E93A87B, 0x00004000 | |
1056 | data8 0x83DD195280A382FB, 0x00004000 | |
1057 | data8 0x82EA3D7FA4CB8C9E, 0x00004000 | |
1058 | data8 0x821B247C6861D0A8, 0x00004000 | |
1059 | data8 0x816EBED163E8D244, 0x00004000 | |
1060 | data8 0x80E42D9127E4CFC6, 0x00004000 | |
1061 | data8 0x807ABF8D28E64AFD, 0x00004000 | |
1062 | data8 0x8031EF26863B4FD8, 0x00004000 | |
1063 | data8 0x800960ADAE8C11FD, 0x00004000 | |
1064 | data8 0x8000E1475FDBEC21, 0x00004000 | |
1065 | data8 0x80186650A07791FA, 0x00004000 | |
1066 | LOCAL_OBJECT_END(tanl_table_scim1) | |
1067 | ||
1068 | Arg = f8 | |
1069 | Save_Norm_Arg = f8 // For input to reduction routine | |
1070 | Result = f8 | |
1071 | r = f8 // For output from reduction routine | |
1072 | c = f9 // For output from reduction routine | |
1073 | U_2 = f10 | |
1074 | rsq = f11 | |
1075 | C_hi = f12 | |
1076 | C_lo = f13 | |
1077 | T_hi = f14 | |
1078 | T_lo = f15 | |
1079 | ||
1080 | d_1 = f33 | |
1081 | N_0 = f34 | |
1082 | tail = f35 | |
1083 | tanx = f36 | |
1084 | Cx = f37 | |
1085 | Sx = f38 | |
1086 | sgn_r = f39 | |
1087 | CORR = f40 | |
1088 | P = f41 | |
1089 | D = f42 | |
1090 | ArgPrime = f43 | |
1091 | P_0 = f44 | |
1092 | ||
1093 | P2_1 = f45 | |
1094 | P2_2 = f46 | |
1095 | P2_3 = f47 | |
1096 | ||
1097 | P1_1 = f45 | |
1098 | P1_2 = f46 | |
1099 | P1_3 = f47 | |
1100 | ||
1101 | P1_4 = f48 | |
1102 | P1_5 = f49 | |
1103 | P1_6 = f50 | |
1104 | P1_7 = f51 | |
1105 | P1_8 = f52 | |
1106 | P1_9 = f53 | |
1107 | ||
1108 | x = f56 | |
1109 | xsq = f57 | |
1110 | Tx = f58 | |
1111 | Tx1 = f59 | |
1112 | Set = f60 | |
1113 | poly1 = f61 | |
1114 | poly2 = f62 | |
1115 | Poly = f63 | |
1116 | Poly1 = f64 | |
1117 | Poly2 = f65 | |
1118 | r_to_the_8 = f66 | |
1119 | B = f67 | |
1120 | SC_inv = f68 | |
1121 | Pos_r = f69 | |
1122 | N_0_fix = f70 | |
1123 | d_2 = f71 | |
1124 | PI_BY_4 = f72 | |
1125 | TWO_TO_NEG14 = f74 | |
1126 | TWO_TO_NEG33 = f75 | |
1127 | NEGTWO_TO_NEG14 = f76 | |
1128 | NEGTWO_TO_NEG33 = f77 | |
1129 | two_by_PI = f78 | |
1130 | N = f79 | |
1131 | N_fix = f80 | |
1132 | P_1 = f81 | |
1133 | P_2 = f82 | |
1134 | P_3 = f83 | |
1135 | s_val = f84 | |
1136 | w = f85 | |
1137 | B_mask1 = f86 | |
1138 | B_mask2 = f87 | |
1139 | w2 = f88 | |
1140 | A = f89 | |
1141 | a = f90 | |
1142 | t = f91 | |
1143 | U_1 = f92 | |
1144 | NEGTWO_TO_NEG2 = f93 | |
1145 | TWO_TO_NEG2 = f94 | |
1146 | Q1_1 = f95 | |
1147 | Q1_2 = f96 | |
1148 | Q1_3 = f97 | |
1149 | Q1_4 = f98 | |
1150 | Q1_5 = f99 | |
1151 | Q1_6 = f100 | |
1152 | Q1_7 = f101 | |
1153 | Q1_8 = f102 | |
1154 | S_hi = f103 | |
1155 | S_lo = f104 | |
1156 | V_hi = f105 | |
1157 | V_lo = f106 | |
1158 | U_hi = f107 | |
1159 | U_lo = f108 | |
1160 | U_hiabs = f109 | |
1161 | V_hiabs = f110 | |
1162 | V = f111 | |
1163 | Inv_P_0 = f112 | |
1164 | ||
1165 | FR_inv_pi_2to63 = f113 | |
1166 | FR_rshf_2to64 = f114 | |
1167 | FR_2tom64 = f115 | |
1168 | FR_rshf = f116 | |
1169 | Norm_Arg = f117 | |
1170 | Abs_Arg = f118 | |
1171 | TWO_TO_NEG65 = f119 | |
1172 | fp_tmp = f120 | |
1173 | mOne = f121 | |
1174 | ||
1175 | GR_SAVE_B0 = r33 | |
1176 | GR_SAVE_GP = r34 | |
1177 | GR_SAVE_PFS = r35 | |
1178 | table_base = r36 | |
1179 | table_ptr1 = r37 | |
1180 | table_ptr2 = r38 | |
1181 | table_ptr3 = r39 | |
1182 | lookup = r40 | |
1183 | N_fix_gr = r41 | |
1184 | GR_exp_2tom2 = r42 | |
1185 | GR_exp_2tom65 = r43 | |
1186 | exp_r = r44 | |
1187 | sig_r = r45 | |
1188 | bmask1 = r46 | |
1189 | table_offset = r47 | |
1190 | bmask2 = r48 | |
1191 | gr_tmp = r49 | |
1192 | cot_flag = r50 | |
1193 | ||
1194 | GR_sig_inv_pi = r51 | |
1195 | GR_rshf_2to64 = r52 | |
1196 | GR_exp_2tom64 = r53 | |
1197 | GR_rshf = r54 | |
1198 | GR_exp_2_to_63 = r55 | |
1199 | GR_exp_2_to_24 = r56 | |
1200 | GR_signexp_x = r57 | |
1201 | GR_exp_x = r58 | |
1202 | GR_exp_mask = r59 | |
1203 | GR_exp_2tom14 = r60 | |
1204 | GR_exp_m2tom14 = r61 | |
1205 | GR_exp_2tom33 = r62 | |
1206 | GR_exp_m2tom33 = r63 | |
1207 | ||
1208 | GR_SAVE_B0 = r64 | |
1209 | GR_SAVE_PFS = r65 | |
1210 | GR_SAVE_GP = r66 | |
1211 | ||
1212 | GR_Parameter_X = r67 | |
1213 | GR_Parameter_Y = r68 | |
1214 | GR_Parameter_RESULT = r69 | |
1215 | GR_Parameter_Tag = r70 | |
1216 | ||
1217 | ||
1218 | .section .text | |
1219 | .global __libm_tanl# | |
1220 | .global __libm_cotl# | |
1221 | ||
1222 | .proc __libm_cotl# | |
1223 | __libm_cotl: | |
1224 | .endp __libm_cotl# | |
1225 | LOCAL_LIBM_ENTRY(cotl) | |
1226 | ||
1227 | { .mlx | |
1228 | alloc r32 = ar.pfs, 0,35,4,0 | |
1229 | movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi | |
1230 | } | |
1231 | { .mlx | |
1232 | mov GR_exp_mask = 0x1ffff // Exponent mask | |
1233 | movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) | |
1234 | } | |
1235 | ;; | |
1236 | ||
1237 | // Check for NatVals, Infs , NaNs, and Zeros | |
1238 | { .mfi | |
1239 | getf.exp GR_signexp_x = Arg // Get sign and exponent of x | |
1240 | fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero | |
1241 | mov cot_flag = 0x1 | |
1242 | } | |
1243 | { .mfb | |
1244 | addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr | |
1245 | fnorm.s1 Norm_Arg = Arg // Normalize x | |
1246 | br.cond.sptk COMMON_PATH | |
1247 | };; | |
1248 | ||
1249 | LOCAL_LIBM_END(cotl) | |
1250 | ||
1251 | ||
1252 | .proc __libm_tanl# | |
1253 | __libm_tanl: | |
1254 | .endp __libm_tanl# | |
1255 | GLOBAL_IEEE754_ENTRY(tanl) | |
1256 | ||
1257 | { .mlx | |
1258 | alloc r32 = ar.pfs, 0,35,4,0 | |
1259 | movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi | |
1260 | } | |
1261 | { .mlx | |
1262 | mov GR_exp_mask = 0x1ffff // Exponent mask | |
1263 | movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) | |
1264 | } | |
1265 | ;; | |
1266 | ||
1267 | // Check for NatVals, Infs , NaNs, and Zeros | |
1268 | { .mfi | |
1269 | getf.exp GR_signexp_x = Arg // Get sign and exponent of x | |
1270 | fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero | |
1271 | mov cot_flag = 0x0 | |
1272 | } | |
1273 | { .mfi | |
1274 | addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr | |
1275 | fnorm.s1 Norm_Arg = Arg // Normalize x | |
1276 | nop.i 0 | |
1277 | };; | |
1278 | ||
1279 | // Common path for both tanl and cotl | |
1280 | COMMON_PATH: | |
1281 | { .mfi | |
1282 | setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63 | |
1283 | fclass.m p9, p0 = Arg, 0x0b // Test x denormal | |
1284 | mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N | |
1285 | } | |
1286 | { .mlx | |
1287 | setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64) | |
1288 | movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63 | |
1289 | } | |
1290 | ;; | |
1291 | ||
1292 | // Check for everything - if false, then must be pseudo-zero or pseudo-nan. | |
1293 | // Branch out to deal with special values. | |
1294 | { .mfi | |
1295 | addl gr_tmp = -1,r0 | |
1296 | fclass.nm p7,p0 = Arg, 0x1FF // Test x unsupported | |
1297 | mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63 | |
1298 | } | |
1299 | { .mfb | |
1300 | ld8 table_base = [table_base] // Get pointer to constant table | |
1301 | fms.s1 mOne = f0, f0, f1 | |
1302 | (p6) br.cond.spnt TANL_SPECIAL // Branch if x natval, nan, inf, zero | |
1303 | } | |
1304 | ;; | |
1305 | ||
1306 | { .mmb | |
1307 | setf.sig fp_tmp = gr_tmp // Make a constant so fmpy produces inexact | |
1308 | mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24 | |
1309 | (p9) br.cond.spnt TANL_DENORMAL // Branch if x denormal | |
1310 | } | |
1311 | ;; | |
1312 | ||
1313 | TANL_COMMON: | |
1314 | // Return to here if x denormal | |
1315 | // | |
1316 | // Do fcmp to generate Denormal exception | |
1317 | // - can't do FNORM (will generate Underflow when U is unmasked!) | |
1318 | // Branch out to deal with unsupporteds values. | |
1319 | { .mfi | |
1320 | setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float | |
1321 | fcmp.eq.s0 p0, p6 = Arg, f1 // Dummy to flag denormals | |
1322 | add table_ptr1 = 0, table_base // Point to tanl_table_1 | |
1323 | } | |
1324 | { .mib | |
1325 | setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63 | |
1326 | add table_ptr2 = 80, table_base // Point to tanl_table_2 | |
1327 | (p7) br.cond.spnt TANL_UNSUPPORTED // Branch if x unsupported type | |
1328 | } | |
1329 | ;; | |
1330 | ||
1331 | { .mfi | |
1332 | and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x | |
1333 | fmpy.s1 Save_Norm_Arg = Norm_Arg, f1 // Save x if large arg reduction | |
1334 | dep.z bmask1 = 0x7c, 56, 8 // Form mask to get 5 msb of r | |
1335 | // bmask1 = 0x7c00000000000000 | |
1336 | } | |
1337 | ;; | |
1338 | ||
1339 | // | |
1340 | // Decide about the paths to take: | |
1341 | // Set PR_6 if |Arg| >= 2**63 | |
1342 | // Set PR_9 if |Arg| < 2**24 - CASE 1 OR 2 | |
1343 | // OTHERWISE Set PR_8 - CASE 3 OR 4 | |
1344 | // | |
1345 | // Branch out if the magnitude of the input argument is >= 2^63 | |
1346 | // - do this branch before the next. | |
1347 | { .mfi | |
1348 | ldfe two_by_PI = [table_ptr1],16 // Load 2/pi | |
1349 | nop.f 999 | |
1350 | dep.z bmask2 = 0x41, 57, 7 // Form mask to OR to produce B | |
1351 | // bmask2 = 0x8200000000000000 | |
1352 | } | |
1353 | { .mib | |
1354 | ldfe PI_BY_4 = [table_ptr2],16 // Load pi/4 | |
1355 | cmp.ge p6,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63 | |
1356 | (p6) br.cond.spnt TANL_ARG_TOO_LARGE // Branch if |x| >= 2^63 | |
1357 | } | |
1358 | ;; | |
1359 | ||
1360 | { .mmi | |
1361 | ldfe P_0 = [table_ptr1],16 // Load P_0 | |
1362 | ldfe Inv_P_0 = [table_ptr2],16 // Load Inv_P_0 | |
1363 | nop.i 999 | |
1364 | } | |
1365 | ;; | |
1366 | ||
1367 | { .mfi | |
1368 | ldfe P_1 = [table_ptr1],16 // Load P_1 | |
1369 | fmerge.s Abs_Arg = f0, Norm_Arg // Get |x| | |
1370 | mov GR_exp_m2tom33 = 0x2ffff - 33 // Form signexp of -2^-33 | |
1371 | } | |
1372 | { .mfi | |
1373 | ldfe d_1 = [table_ptr2],16 // Load d_1 for 2^24 <= |x| < 2^63 | |
1374 | nop.f 999 | |
1375 | mov GR_exp_2tom33 = 0xffff - 33 // Form signexp of 2^-33 | |
1376 | } | |
1377 | ;; | |
1378 | ||
1379 | { .mmi | |
1380 | ldfe P_2 = [table_ptr1],16 // Load P_2 | |
1381 | ldfe d_2 = [table_ptr2],16 // Load d_2 for 2^24 <= |x| < 2^63 | |
1382 | cmp.ge p8,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24 | |
1383 | } | |
1384 | ;; | |
1385 | ||
1386 | // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits | |
1387 | // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24 | |
1388 | { .mfb | |
1389 | ldfe P_3 = [table_ptr1],16 // Load P_3 | |
1390 | fma.s1 N_fix = Norm_Arg, FR_inv_pi_2to63, FR_rshf_2to64 | |
1391 | (p8) br.cond.spnt TANL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63 | |
1392 | } | |
1393 | ;; | |
1394 | ||
1395 | // Here if 0 < |x| < 2^24 | |
1396 | // ARGUMENT REDUCTION CODE - CASE 1 and 2 | |
1397 | // | |
1398 | { .mmf | |
1399 | setf.exp TWO_TO_NEG33 = GR_exp_2tom33 // Form 2^-33 | |
1400 | setf.exp NEGTWO_TO_NEG33 = GR_exp_m2tom33 // Form -2^-33 | |
1401 | fmerge.s r = Norm_Arg,Norm_Arg // Assume r=x, ok if |x| < pi/4 | |
1402 | } | |
1403 | ;; | |
1404 | ||
1405 | // | |
1406 | // If |Arg| < pi/4, set PR_8, else pi/4 <=|Arg| < 2^24 - set PR_9. | |
1407 | // | |
1408 | // Case 2: Convert integer N_fix back to normalized floating-point value. | |
1409 | { .mfi | |
1410 | getf.sig sig_r = Norm_Arg // Get sig_r if 1/4 <= |x| < pi/4 | |
1411 | fcmp.lt.s1 p8,p9= Abs_Arg,PI_BY_4 // Test |x| < pi/4 | |
1412 | mov GR_exp_2tom2 = 0xffff - 2 // Form signexp of 2^-2 | |
1413 | } | |
1414 | { .mfi | |
1415 | ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] // Load 2^-2, -2^-2 | |
1416 | fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated | |
1417 | mov N_fix_gr = r0 // Assume N=0, ok if |x| < pi/4 | |
1418 | } | |
1419 | ;; | |
1420 | ||
1421 | // | |
1422 | // Case 1: Is |r| < 2**(-2). | |
1423 | // Arg is the same as r in this case. | |
1424 | // r = Arg | |
1425 | // c = 0 | |
1426 | // | |
1427 | // Case 2: Place integer part of N in GP register. | |
1428 | { .mfi | |
1429 | (p9) getf.sig N_fix_gr = N_fix | |
1430 | fmerge.s c = f0, f0 // Assume c=0, ok if |x| < pi/4 | |
1431 | cmp.lt p10, p0 = GR_exp_x, GR_exp_2tom2 // Test if |x| < 1/4 | |
1432 | } | |
1433 | ;; | |
1434 | ||
1435 | { .mfi | |
1436 | setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r | |
1437 | nop.f 999 | |
1438 | mov exp_r = GR_exp_x // Get exp_r if 1/4 <= |x| < pi/4 | |
1439 | } | |
1440 | { .mbb | |
1441 | setf.sig B_mask2 = bmask2 // Form mask to form B from r | |
1442 | (p10) br.cond.spnt TANL_SMALL_R // Branch if 0 < |x| < 1/4 | |
1443 | (p8) br.cond.spnt TANL_NORMAL_R // Branch if 1/4 <= |x| < pi/4 | |
1444 | } | |
1445 | ;; | |
1446 | ||
1447 | // Here if pi/4 <= |x| < 2^24 | |
1448 | // | |
1449 | // Case 1: PR_3 is only affected when PR_1 is set. | |
1450 | // | |
1451 | // | |
1452 | // Case 2: w = N * P_2 | |
1453 | // Case 2: s_val = -N * P_1 + Arg | |
1454 | // | |
1455 | ||
1456 | { .mfi | |
1457 | nop.m 999 | |
1458 | fnma.s1 s_val = N, P_1, Norm_Arg | |
1459 | nop.i 999 | |
1460 | } | |
1461 | { .mfi | |
1462 | nop.m 999 | |
1463 | fmpy.s1 w = N, P_2 // w = N * P_2 for |s| >= 2^-33 | |
1464 | nop.i 999 | |
1465 | } | |
1466 | ;; | |
1467 | ||
1468 | // Case 2_reduce: w = N * P_3 (change sign) | |
1469 | { .mfi | |
1470 | nop.m 999 | |
1471 | fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-33 | |
1472 | nop.i 999 | |
1473 | } | |
1474 | ;; | |
1475 | ||
1476 | // Case 1_reduce: r = s + w (change sign) | |
1477 | { .mfi | |
1478 | nop.m 999 | |
1479 | fsub.s1 r = s_val, w // r = s_val - w for |s| >= 2^-33 | |
1480 | nop.i 999 | |
1481 | } | |
1482 | ;; | |
1483 | ||
1484 | // Case 2_reduce: U_1 = N * P_2 + w | |
1485 | { .mfi | |
1486 | nop.m 999 | |
1487 | fma.s1 U_1 = N, P_2, w2 // U_1 = N * P_2 + w for |s| < 2^-33 | |
1488 | nop.i 999 | |
1489 | } | |
1490 | ;; | |
1491 | ||
1492 | // | |
1493 | // Decide between case_1 and case_2 reduce: | |
1494 | // Case 1_reduce: |s| >= 2**(-33) | |
1495 | // Case 2_reduce: |s| < 2**(-33) | |
1496 | // | |
1497 | { .mfi | |
1498 | nop.m 999 | |
1499 | fcmp.lt.s1 p9, p8 = s_val, TWO_TO_NEG33 | |
1500 | nop.i 999 | |
1501 | } | |
1502 | ;; | |
1503 | ||
1504 | { .mfi | |
1505 | nop.m 999 | |
1506 | (p9) fcmp.gt.s1 p9, p8 = s_val, NEGTWO_TO_NEG33 | |
1507 | nop.i 999 | |
1508 | } | |
1509 | ;; | |
1510 | ||
1511 | // Case 1_reduce: c = s - r | |
1512 | { .mfi | |
1513 | nop.m 999 | |
1514 | fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-33 | |
1515 | nop.i 999 | |
1516 | } | |
1517 | ;; | |
1518 | ||
1519 | // Case 2_reduce: r is complete here - continue to calculate c . | |
1520 | // r = s - U_1 | |
1521 | { .mfi | |
1522 | nop.m 999 | |
1523 | (p9) fsub.s1 r = s_val, U_1 | |
1524 | nop.i 999 | |
1525 | } | |
1526 | { .mfi | |
1527 | nop.m 999 | |
1528 | (p9) fms.s1 U_2 = N, P_2, U_1 | |
1529 | nop.i 999 | |
1530 | } | |
1531 | ;; | |
1532 | ||
1533 | // | |
1534 | // Case 1_reduce: Is |r| < 2**(-2), if so set PR_10 | |
1535 | // else set PR_13. | |
1536 | // | |
1537 | ||
1538 | { .mfi | |
1539 | nop.m 999 | |
1540 | fand B = B_mask1, r | |
1541 | nop.i 999 | |
1542 | } | |
1543 | { .mfi | |
1544 | nop.m 999 | |
1545 | (p8) fcmp.lt.unc.s1 p10, p13 = r, TWO_TO_NEG2 | |
1546 | nop.i 999 | |
1547 | } | |
1548 | ;; | |
1549 | ||
1550 | { .mfi | |
1551 | (p8) getf.sig sig_r = r // Get signif of r if |s| >= 2^-33 | |
1552 | nop.f 999 | |
1553 | nop.i 999 | |
1554 | } | |
1555 | ;; | |
1556 | ||
1557 | { .mfi | |
1558 | (p8) getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33 | |
1559 | (p10) fcmp.gt.s1 p10, p13 = r, NEGTWO_TO_NEG2 | |
1560 | nop.i 999 | |
1561 | } | |
1562 | ;; | |
1563 | ||
1564 | // Case 1_reduce: c is complete here. | |
1565 | // Case 1: Branch to SMALL_R or NORMAL_R. | |
1566 | // c = c + w (w has not been negated.) | |
1567 | { .mfi | |
1568 | nop.m 999 | |
1569 | (p8) fsub.s1 c = c, w // c = c - w for |s| >= 2^-33 | |
1570 | nop.i 999 | |
1571 | } | |
1572 | { .mbb | |
1573 | nop.m 999 | |
1574 | (p10) br.cond.spnt TANL_SMALL_R // Branch if pi/4 < |x| < 2^24 and |r|<1/4 | |
1575 | (p13) br.cond.sptk TANL_NORMAL_R_A // Branch if pi/4 < |x| < 2^24 and |r|>=1/4 | |
1576 | } | |
1577 | ;; | |
1578 | ||
1579 | ||
1580 | // Here if pi/4 < |x| < 2^24 and |s| < 2^-33 | |
1581 | // | |
1582 | // Is i_1 = lsb of N_fix_gr even or odd? | |
1583 | // if i_1 == 0, set p11, else set p12. | |
1584 | // | |
1585 | { .mfi | |
1586 | nop.m 999 | |
1587 | fsub.s1 s_val = s_val, r | |
1588 | add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl) | |
1589 | } | |
1590 | { .mfi | |
1591 | nop.m 999 | |
1592 | // | |
1593 | // Case 2_reduce: | |
1594 | // U_2 = N * P_2 - U_1 | |
1595 | // Not needed until later. | |
1596 | // | |
1597 | fadd.s1 U_2 = U_2, w2 | |
1598 | // | |
1599 | // Case 2_reduce: | |
1600 | // s = s - r | |
1601 | // U_2 = U_2 + w | |
1602 | // | |
1603 | nop.i 999 | |
1604 | } | |
1605 | ;; | |
1606 | ||
1607 | // | |
1608 | // Case 2_reduce: | |
1609 | // c = c - U_2 | |
1610 | // c is complete here | |
1611 | // Argument reduction ends here. | |
1612 | // | |
1613 | { .mfi | |
1614 | nop.m 999 | |
1615 | fmpy.s1 rsq = r, r | |
1616 | tbit.z p11, p12 = N_fix_gr, 0 ;; // Set p11 if N even, p12 if odd | |
1617 | } | |
1618 | ||
1619 | { .mfi | |
1620 | nop.m 999 | |
1621 | (p12) frcpa.s1 S_hi,p0 = f1, r | |
1622 | nop.i 999 | |
1623 | } | |
1624 | { .mfi | |
1625 | nop.m 999 | |
1626 | fsub.s1 c = s_val, U_1 | |
1627 | nop.i 999 | |
1628 | } | |
1629 | ;; | |
1630 | ||
1631 | { .mmi | |
1632 | add table_ptr1 = 160, table_base ;; // Point to tanl_table_p1 | |
1633 | ldfe P1_1 = [table_ptr1],144 | |
1634 | nop.i 999 ;; | |
1635 | } | |
1636 | // | |
1637 | // Load P1_1 and point to Q1_1 . | |
1638 | // | |
1639 | { .mfi | |
1640 | ldfe Q1_1 = [table_ptr1] | |
1641 | // | |
1642 | // N even: rsq = r * Z | |
1643 | // N odd: S_hi = frcpa(r) | |
1644 | // | |
1645 | (p12) fmerge.ns S_hi = S_hi, S_hi | |
1646 | nop.i 999 | |
1647 | } | |
1648 | { .mfi | |
1649 | nop.m 999 | |
1650 | // | |
1651 | // Case 2_reduce: | |
1652 | // c = s - U_1 | |
1653 | // | |
1654 | (p9) fsub.s1 c = c, U_2 | |
1655 | nop.i 999 ;; | |
1656 | } | |
1657 | { .mfi | |
1658 | nop.m 999 | |
1659 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
1660 | nop.i 999 ;; | |
1661 | } | |
1662 | { .mfi | |
1663 | nop.m 999 | |
1664 | // | |
1665 | // N odd: Change sign of S_hi | |
1666 | // | |
1667 | (p11) fmpy.s1 rsq = rsq, P1_1 | |
1668 | nop.i 999 ;; | |
1669 | } | |
1670 | { .mfi | |
1671 | nop.m 999 | |
1672 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
1673 | nop.i 999 ;; | |
1674 | } | |
1675 | { .mfi | |
1676 | nop.m 999 | |
1677 | // | |
1678 | // N even: rsq = rsq * P1_1 | |
1679 | // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary | |
1680 | // | |
1681 | (p11) fma.s1 Poly = r, rsq, c | |
1682 | nop.i 999 ;; | |
1683 | } | |
1684 | { .mfi | |
1685 | nop.m 999 | |
1686 | // | |
1687 | // N even: Poly = c + r * rsq | |
1688 | // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary | |
1689 | // | |
1690 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
1691 | (p11) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl | |
1692 | } | |
1693 | { .mfi | |
1694 | nop.m 999 | |
1695 | // | |
1696 | // N even: Result = Poly + r | |
1697 | // N odd: poly1 = 1.0 + S_hi * r 32 bits partial | |
1698 | // | |
1699 | (p14) fadd.s0 Result = r, Poly // for tanl | |
1700 | nop.i 999 | |
1701 | } | |
1702 | { .mfi | |
1703 | nop.m 999 | |
1704 | (p15) fms.s0 Result = r, mOne, Poly // for cotl | |
1705 | nop.i 999 | |
1706 | } | |
1707 | ;; | |
1708 | ||
1709 | { .mfi | |
1710 | nop.m 999 | |
1711 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
1712 | nop.i 999 ;; | |
1713 | } | |
1714 | { .mfi | |
1715 | nop.m 999 | |
1716 | // | |
1717 | // N even: Result1 = Result + r | |
1718 | // N odd: S_hi = S_hi * poly1 + S_hi 32 bits | |
1719 | // | |
1720 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
1721 | nop.i 999 ;; | |
1722 | } | |
1723 | { .mfi | |
1724 | nop.m 999 | |
1725 | // | |
1726 | // N odd: poly1 = S_hi * r + 1.0 64 bits partial | |
1727 | // | |
1728 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
1729 | nop.i 999 ;; | |
1730 | } | |
1731 | { .mfi | |
1732 | nop.m 999 | |
1733 | // | |
1734 | // N odd: poly1 = S_hi * poly + 1.0 64 bits | |
1735 | // | |
1736 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
1737 | nop.i 999 ;; | |
1738 | } | |
1739 | { .mfi | |
1740 | nop.m 999 | |
1741 | // | |
1742 | // N odd: poly1 = S_hi * r + 1.0 | |
1743 | // | |
1744 | (p12) fma.s1 poly1 = S_hi, c, poly1 | |
1745 | nop.i 999 ;; | |
1746 | } | |
1747 | { .mfi | |
1748 | nop.m 999 | |
1749 | // | |
1750 | // N odd: poly1 = S_hi * c + poly1 | |
1751 | // | |
1752 | (p12) fmpy.s1 S_lo = S_hi, poly1 | |
1753 | nop.i 999 ;; | |
1754 | } | |
1755 | { .mfi | |
1756 | nop.m 999 | |
1757 | // | |
1758 | // N odd: S_lo = S_hi * poly1 | |
1759 | // | |
1760 | (p12) fma.s1 S_lo = Q1_1, r, S_lo | |
1761 | (p12) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl | |
1762 | } | |
1763 | { .mfi | |
1764 | nop.m 999 | |
1765 | // | |
1766 | // N odd: Result = S_hi + S_lo | |
1767 | // | |
1768 | fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact | |
1769 | nop.i 999 ;; | |
1770 | } | |
1771 | { .mfi | |
1772 | nop.m 999 | |
1773 | // | |
1774 | // N odd: S_lo = S_lo + Q1_1 * r | |
1775 | // | |
1776 | (p14) fadd.s0 Result = S_hi, S_lo // for tanl | |
1777 | nop.i 999 | |
1778 | } | |
1779 | { .mfb | |
1780 | nop.m 999 | |
1781 | (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl | |
1782 | br.ret.sptk b0 ;; // Exit for pi/4 <= |x| < 2^24 and |s| < 2^-33 | |
1783 | } | |
1784 | ||
1785 | ||
1786 | TANL_LARGER_ARG: | |
1787 | // Here if 2^24 <= |x| < 2^63 | |
1788 | // | |
1789 | // ARGUMENT REDUCTION CODE - CASE 3 and 4 | |
1790 | // | |
1791 | ||
1792 | { .mmf | |
1793 | mov GR_exp_2tom14 = 0xffff - 14 // Form signexp of 2^-14 | |
1794 | mov GR_exp_m2tom14 = 0x2ffff - 14 // Form signexp of -2^-14 | |
1795 | fmpy.s1 N_0 = Norm_Arg, Inv_P_0 | |
1796 | } | |
1797 | ;; | |
1798 | ||
1799 | { .mmi | |
1800 | setf.exp TWO_TO_NEG14 = GR_exp_2tom14 // Form 2^-14 | |
1801 | setf.exp NEGTWO_TO_NEG14 = GR_exp_m2tom14// Form -2^-14 | |
1802 | nop.i 999 | |
1803 | } | |
1804 | ;; | |
1805 | ||
1806 | ||
1807 | // | |
1808 | // Adjust table_ptr1 to beginning of table. | |
1809 | // N_0 = Arg * Inv_P_0 | |
1810 | // | |
1811 | { .mmi | |
1812 | add table_ptr2 = 144, table_base ;; // Point to 2^-2 | |
1813 | ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] | |
1814 | nop.i 999 | |
1815 | } | |
1816 | ;; | |
1817 | ||
1818 | // | |
1819 | // N_0_fix = integer part of N_0 . | |
1820 | // | |
1821 | // | |
1822 | // Make N_0 the integer part. | |
1823 | // | |
1824 | { .mfi | |
1825 | nop.m 999 | |
1826 | fcvt.fx.s1 N_0_fix = N_0 | |
1827 | nop.i 999 ;; | |
1828 | } | |
1829 | { .mfi | |
1830 | setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r | |
1831 | fcvt.xf N_0 = N_0_fix | |
1832 | nop.i 999 ;; | |
1833 | } | |
1834 | { .mfi | |
1835 | setf.sig B_mask2 = bmask2 // Form mask to form B from r | |
1836 | fnma.s1 ArgPrime = N_0, P_0, Norm_Arg | |
1837 | nop.i 999 | |
1838 | } | |
1839 | { .mfi | |
1840 | nop.m 999 | |
1841 | fmpy.s1 w = N_0, d_1 | |
1842 | nop.i 999 ;; | |
1843 | } | |
1844 | // | |
1845 | // ArgPrime = -N_0 * P_0 + Arg | |
1846 | // w = N_0 * d_1 | |
1847 | // | |
1848 | // | |
1849 | // N = ArgPrime * 2/pi | |
1850 | // | |
1851 | // fcvt.fx.s1 N_fix = N | |
1852 | // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits | |
1853 | // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24 | |
1854 | { .mfi | |
1855 | nop.m 999 | |
1856 | fma.s1 N_fix = ArgPrime, FR_inv_pi_2to63, FR_rshf_2to64 | |
1857 | ||
1858 | nop.i 999 ;; | |
1859 | } | |
1860 | // Convert integer N_fix back to normalized floating-point value. | |
1861 | { .mfi | |
1862 | nop.m 999 | |
1863 | fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated | |
1864 | nop.i 999 | |
1865 | } | |
1866 | ;; | |
1867 | ||
1868 | // | |
1869 | // N is the integer part of the reduced-reduced argument. | |
1870 | // Put the integer in a GP register. | |
1871 | // | |
1872 | { .mfi | |
1873 | getf.sig N_fix_gr = N_fix | |
1874 | nop.f 999 | |
1875 | nop.i 999 | |
1876 | } | |
1877 | ;; | |
1878 | ||
1879 | // | |
1880 | // s_val = -N*P_1 + ArgPrime | |
1881 | // w = -N*P_2 + w | |
1882 | // | |
1883 | { .mfi | |
1884 | nop.m 999 | |
1885 | fnma.s1 s_val = N, P_1, ArgPrime | |
1886 | nop.i 999 | |
1887 | } | |
1888 | { .mfi | |
1889 | nop.m 999 | |
1890 | fnma.s1 w = N, P_2, w | |
1891 | nop.i 999 | |
1892 | } | |
1893 | ;; | |
1894 | ||
1895 | // Case 4: V_hi = N * P_2 | |
1896 | // Case 4: U_hi = N_0 * d_1 | |
1897 | { .mfi | |
1898 | nop.m 999 | |
1899 | fmpy.s1 V_hi = N, P_2 // V_hi = N * P_2 for |s| < 2^-14 | |
1900 | nop.i 999 | |
1901 | } | |
1902 | { .mfi | |
1903 | nop.m 999 | |
1904 | fmpy.s1 U_hi = N_0, d_1 // U_hi = N_0 * d_1 for |s| < 2^-14 | |
1905 | nop.i 999 | |
1906 | } | |
1907 | ;; | |
1908 | ||
1909 | // Case 3: r = s_val + w (Z complete) | |
1910 | // Case 4: w = N * P_3 | |
1911 | { .mfi | |
1912 | nop.m 999 | |
1913 | fadd.s1 r = s_val, w // r = s_val + w for |s| >= 2^-14 | |
1914 | nop.i 999 | |
1915 | } | |
1916 | { .mfi | |
1917 | nop.m 999 | |
1918 | fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-14 | |
1919 | nop.i 999 | |
1920 | } | |
1921 | ;; | |
1922 | ||
1923 | // Case 4: A = U_hi + V_hi | |
1924 | // Note: Worry about switched sign of V_hi, so subtract instead of add. | |
1925 | // Case 4: V_lo = -N * P_2 - V_hi (U_hi is in place of V_hi in writeup) | |
1926 | // Note: the (-) is still missing for V_hi. | |
1927 | { .mfi | |
1928 | nop.m 999 | |
1929 | fsub.s1 A = U_hi, V_hi // A = U_hi - V_hi for |s| < 2^-14 | |
1930 | nop.i 999 | |
1931 | } | |
1932 | { .mfi | |
1933 | nop.m 999 | |
1934 | fnma.s1 V_lo = N, P_2, V_hi // V_lo = V_hi - N * P_2 for |s| < 2^-14 | |
1935 | nop.i 999 | |
1936 | } | |
1937 | ;; | |
1938 | ||
1939 | // Decide between case 3 and 4: | |
1940 | // Case 3: |s| >= 2**(-14) Set p10 | |
1941 | // Case 4: |s| < 2**(-14) Set p11 | |
1942 | // | |
1943 | // Case 4: U_lo = N_0 * d_1 - U_hi | |
1944 | { .mfi | |
1945 | nop.m 999 | |
1946 | fms.s1 U_lo = N_0, d_1, U_hi // U_lo = N_0*d_1 - U_hi for |s| < 2^-14 | |
1947 | nop.i 999 | |
1948 | } | |
1949 | { .mfi | |
1950 | nop.m 999 | |
1951 | fcmp.lt.s1 p11, p10 = s_val, TWO_TO_NEG14 | |
1952 | nop.i 999 | |
1953 | } | |
1954 | ;; | |
1955 | ||
1956 | // Case 4: We need abs of both U_hi and V_hi - dont | |
1957 | // worry about switched sign of V_hi. | |
1958 | { .mfi | |
1959 | nop.m 999 | |
1960 | fabs V_hiabs = V_hi // |V_hi| for |s| < 2^-14 | |
1961 | nop.i 999 | |
1962 | } | |
1963 | { .mfi | |
1964 | nop.m 999 | |
1965 | (p11) fcmp.gt.s1 p11, p10 = s_val, NEGTWO_TO_NEG14 | |
1966 | nop.i 999 | |
1967 | } | |
1968 | ;; | |
1969 | ||
1970 | // Case 3: c = s_val - r | |
1971 | { .mfi | |
1972 | nop.m 999 | |
1973 | fabs U_hiabs = U_hi // |U_hi| for |s| < 2^-14 | |
1974 | nop.i 999 | |
1975 | } | |
1976 | { .mfi | |
1977 | nop.m 999 | |
1978 | fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-14 | |
1979 | nop.i 999 | |
1980 | } | |
1981 | ;; | |
1982 | ||
1983 | // For Case 3, |s| >= 2^-14, determine if |r| < 1/4 | |
1984 | // | |
1985 | // Case 4: C_hi = s_val + A | |
1986 | // | |
1987 | { .mfi | |
1988 | nop.m 999 | |
1989 | (p11) fadd.s1 C_hi = s_val, A // C_hi = s_val + A for |s| < 2^-14 | |
1990 | nop.i 999 | |
1991 | } | |
1992 | { .mfi | |
1993 | nop.m 999 | |
1994 | (p10) fcmp.lt.unc.s1 p14, p15 = r, TWO_TO_NEG2 | |
1995 | nop.i 999 | |
1996 | } | |
1997 | ;; | |
1998 | ||
1999 | { .mfi | |
2000 | getf.sig sig_r = r // Get signif of r if |s| >= 2^-33 | |
2001 | fand B = B_mask1, r | |
2002 | nop.i 999 | |
2003 | } | |
2004 | ;; | |
2005 | ||
2006 | // Case 4: t = U_lo + V_lo | |
2007 | { .mfi | |
2008 | getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33 | |
2009 | (p11) fadd.s1 t = U_lo, V_lo // t = U_lo + V_lo for |s| < 2^-14 | |
2010 | nop.i 999 | |
2011 | } | |
2012 | { .mfi | |
2013 | nop.m 999 | |
2014 | (p14) fcmp.gt.s1 p14, p15 = r, NEGTWO_TO_NEG2 | |
2015 | nop.i 999 | |
2016 | } | |
2017 | ;; | |
2018 | ||
2019 | // Case 3: c = (s - r) + w (c complete) | |
2020 | { .mfi | |
2021 | nop.m 999 | |
2022 | (p10) fadd.s1 c = c, w // c = c + w for |s| >= 2^-14 | |
2023 | nop.i 999 | |
2024 | } | |
2025 | { .mbb | |
2026 | nop.m 999 | |
2027 | (p14) br.cond.spnt TANL_SMALL_R // Branch if 2^24 <= |x| < 2^63 and |r|< 1/4 | |
2028 | (p15) br.cond.sptk TANL_NORMAL_R_A // Branch if 2^24 <= |x| < 2^63 and |r|>=1/4 | |
2029 | } | |
2030 | ;; | |
2031 | ||
2032 | ||
2033 | // Here if 2^24 <= |x| < 2^63 and |s| < 2^-14 >>>>>>> Case 4. | |
2034 | // | |
2035 | // Case 4: Set P_12 if U_hiabs >= V_hiabs | |
2036 | // Case 4: w = w + N_0 * d_2 | |
2037 | // Note: the (-) is now incorporated in w . | |
2038 | { .mfi | |
2039 | add table_ptr1 = 160, table_base // Point to tanl_table_p1 | |
2040 | fcmp.ge.unc.s1 p12, p13 = U_hiabs, V_hiabs | |
2041 | nop.i 999 | |
2042 | } | |
2043 | { .mfi | |
2044 | nop.m 999 | |
2045 | fms.s1 w2 = N_0, d_2, w2 | |
2046 | nop.i 999 | |
2047 | } | |
2048 | ;; | |
2049 | ||
2050 | // Case 4: C_lo = s_val - C_hi | |
2051 | { .mfi | |
2052 | ldfe P1_1 = [table_ptr1], 16 // Load P1_1 | |
2053 | fsub.s1 C_lo = s_val, C_hi | |
2054 | nop.i 999 | |
2055 | } | |
2056 | ;; | |
2057 | ||
2058 | // | |
2059 | // Case 4: a = U_hi - A | |
2060 | // a = V_hi - A (do an add to account for missing (-) on V_hi | |
2061 | // | |
2062 | { .mfi | |
2063 | ldfe P1_2 = [table_ptr1], 128 // Load P1_2 | |
2064 | (p12) fsub.s1 a = U_hi, A | |
2065 | nop.i 999 | |
2066 | } | |
2067 | { .mfi | |
2068 | nop.m 999 | |
2069 | (p13) fadd.s1 a = V_hi, A | |
2070 | nop.i 999 | |
2071 | } | |
2072 | ;; | |
2073 | ||
2074 | // Case 4: t = U_lo + V_lo + w | |
2075 | { .mfi | |
2076 | ldfe Q1_1 = [table_ptr1], 16 // Load Q1_1 | |
2077 | fadd.s1 t = t, w2 | |
2078 | nop.i 999 | |
2079 | } | |
2080 | ;; | |
2081 | ||
2082 | // Case 4: a = (U_hi - A) + V_hi | |
2083 | // a = (V_hi - A) + U_hi | |
2084 | // In each case account for negative missing form V_hi . | |
2085 | // | |
2086 | { .mfi | |
2087 | ldfe Q1_2 = [table_ptr1], 16 // Load Q1_2 | |
2088 | (p12) fsub.s1 a = a, V_hi | |
2089 | nop.i 999 | |
2090 | } | |
2091 | { .mfi | |
2092 | nop.m 999 | |
2093 | (p13) fsub.s1 a = U_hi, a | |
2094 | nop.i 999 | |
2095 | } | |
2096 | ;; | |
2097 | ||
2098 | // | |
2099 | // Case 4: C_lo = (s_val - C_hi) + A | |
2100 | // | |
2101 | { .mfi | |
2102 | nop.m 999 | |
2103 | fadd.s1 C_lo = C_lo, A | |
2104 | nop.i 999 ;; | |
2105 | } | |
2106 | // | |
2107 | // Case 4: t = t + a | |
2108 | // | |
2109 | { .mfi | |
2110 | nop.m 999 | |
2111 | fadd.s1 t = t, a | |
2112 | nop.i 999 | |
2113 | } | |
2114 | ;; | |
2115 | ||
2116 | // Case 4: C_lo = C_lo + t | |
2117 | // Case 4: r = C_hi + C_lo | |
2118 | { .mfi | |
2119 | nop.m 999 | |
2120 | fadd.s1 C_lo = C_lo, t | |
2121 | nop.i 999 | |
2122 | } | |
2123 | ;; | |
2124 | ||
2125 | { .mfi | |
2126 | nop.m 999 | |
2127 | fadd.s1 r = C_hi, C_lo | |
2128 | nop.i 999 | |
2129 | } | |
2130 | ;; | |
2131 | ||
2132 | // | |
2133 | // Case 4: c = C_hi - r | |
2134 | // | |
2135 | { .mfi | |
2136 | nop.m 999 | |
2137 | fsub.s1 c = C_hi, r | |
2138 | nop.i 999 | |
2139 | } | |
2140 | { .mfi | |
2141 | nop.m 999 | |
2142 | fmpy.s1 rsq = r, r | |
2143 | add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl) | |
2144 | } | |
2145 | ;; | |
2146 | ||
2147 | // Case 4: c = c + C_lo finished. | |
2148 | // | |
2149 | // Is i_1 = lsb of N_fix_gr even or odd? | |
2150 | // if i_1 == 0, set PR_11, else set PR_12. | |
2151 | // | |
2152 | { .mfi | |
2153 | nop.m 999 | |
2154 | fadd.s1 c = c , C_lo | |
2155 | tbit.z p11, p12 = N_fix_gr, 0 | |
2156 | } | |
2157 | ;; | |
2158 | ||
2159 | // r and c have been computed. | |
2160 | { .mfi | |
2161 | nop.m 999 | |
2162 | (p12) frcpa.s1 S_hi, p0 = f1, r | |
2163 | nop.i 999 | |
2164 | } | |
2165 | { .mfi | |
2166 | nop.m 999 | |
2167 | // | |
2168 | // N odd: Change sign of S_hi | |
2169 | // | |
2170 | (p11) fma.s1 Poly = rsq, P1_2, P1_1 | |
2171 | nop.i 999 ;; | |
2172 | } | |
2173 | { .mfi | |
2174 | nop.m 999 | |
2175 | (p12) fma.s1 P = rsq, Q1_2, Q1_1 | |
2176 | nop.i 999 | |
2177 | } | |
2178 | { .mfi | |
2179 | nop.m 999 | |
2180 | // | |
2181 | // N odd: Result = S_hi + S_lo (User supplied rounding mode for C1) | |
2182 | // | |
2183 | fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact | |
2184 | nop.i 999 ;; | |
2185 | } | |
2186 | { .mfi | |
2187 | nop.m 999 | |
2188 | // | |
2189 | // N even: rsq = r * r | |
2190 | // N odd: S_hi = frcpa(r) | |
2191 | // | |
2192 | (p12) fmerge.ns S_hi = S_hi, S_hi | |
2193 | nop.i 999 | |
2194 | } | |
2195 | { .mfi | |
2196 | nop.m 999 | |
2197 | // | |
2198 | // N even: rsq = rsq * P1_2 + P1_1 | |
2199 | // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary | |
2200 | // | |
2201 | (p11) fmpy.s1 Poly = rsq, Poly | |
2202 | nop.i 999 ;; | |
2203 | } | |
2204 | { .mfi | |
2205 | nop.m 999 | |
2206 | (p12) fma.s1 poly1 = S_hi, r,f1 | |
2207 | (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl | |
2208 | } | |
2209 | { .mfi | |
2210 | nop.m 999 | |
2211 | // | |
2212 | // N even: Poly = Poly * rsq | |
2213 | // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary | |
2214 | // | |
2215 | (p11) fma.s1 Poly = r, Poly, c | |
2216 | nop.i 999 ;; | |
2217 | } | |
2218 | { .mfi | |
2219 | nop.m 999 | |
2220 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2221 | nop.i 999 | |
2222 | } | |
2223 | { .mfi | |
2224 | nop.m 999 | |
2225 | // | |
2226 | // N odd: S_hi = S_hi * poly1 + S_hi 32 bits | |
2227 | // | |
2228 | (p14) fadd.s0 Result = r, Poly // for tanl | |
2229 | nop.i 999 ;; | |
2230 | } | |
2231 | ||
2232 | .pred.rel "mutex",p15,p12 | |
2233 | { .mfi | |
2234 | nop.m 999 | |
2235 | (p15) fms.s0 Result = r, mOne, Poly // for cotl | |
2236 | nop.i 999 | |
2237 | } | |
2238 | { .mfi | |
2239 | nop.m 999 | |
2240 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2241 | nop.i 999 ;; | |
2242 | } | |
2243 | { .mfi | |
2244 | nop.m 999 | |
2245 | // | |
2246 | // N even: Poly = Poly * r + c | |
2247 | // N odd: poly1 = 1.0 + S_hi * r 32 bits partial | |
2248 | // | |
2249 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2250 | nop.i 999 ;; | |
2251 | } | |
2252 | { .mfi | |
2253 | nop.m 999 | |
2254 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2255 | nop.i 999 ;; | |
2256 | } | |
2257 | { .mfi | |
2258 | nop.m 999 | |
2259 | // | |
2260 | // N even: Result = Poly + r (Rounding mode S0) | |
2261 | // N odd: poly1 = S_hi * r + 1.0 64 bits partial | |
2262 | // | |
2263 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2264 | nop.i 999 ;; | |
2265 | } | |
2266 | { .mfi | |
2267 | nop.m 999 | |
2268 | // | |
2269 | // N odd: poly1 = S_hi * poly + S_hi 64 bits | |
2270 | // | |
2271 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2272 | nop.i 999 ;; | |
2273 | } | |
2274 | { .mfi | |
2275 | nop.m 999 | |
2276 | // | |
2277 | // N odd: poly1 = S_hi * r + 1.0 | |
2278 | // | |
2279 | (p12) fma.s1 poly1 = S_hi, c, poly1 | |
2280 | nop.i 999 ;; | |
2281 | } | |
2282 | { .mfi | |
2283 | nop.m 999 | |
2284 | // | |
2285 | // N odd: poly1 = S_hi * c + poly1 | |
2286 | // | |
2287 | (p12) fmpy.s1 S_lo = S_hi, poly1 | |
2288 | nop.i 999 ;; | |
2289 | } | |
2290 | { .mfi | |
2291 | nop.m 999 | |
2292 | // | |
2293 | // N odd: S_lo = S_hi * poly1 | |
2294 | // | |
2295 | (p12) fma.s1 S_lo = P, r, S_lo | |
2296 | (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl | |
2297 | } | |
2298 | ||
2299 | { .mfi | |
2300 | nop.m 999 | |
2301 | (p14) fadd.s0 Result = S_hi, S_lo // for tanl | |
2302 | nop.i 999 | |
2303 | } | |
2304 | { .mfb | |
2305 | nop.m 999 | |
2306 | // | |
2307 | // N odd: S_lo = S_lo + r * P | |
2308 | // | |
2309 | (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl | |
2310 | br.ret.sptk b0 ;; // Exit for 2^24 <= |x| < 2^63 and |s| < 2^-14 | |
2311 | } | |
2312 | ||
2313 | ||
2314 | TANL_SMALL_R: | |
2315 | // Here if |r| < 1/4 | |
2316 | // r and c have been computed. | |
2317 | // ***************************************************************** | |
2318 | // ***************************************************************** | |
2319 | // ***************************************************************** | |
2320 | // N odd: S_hi = frcpa(r) | |
2321 | // Get [i_1] - lsb of N_fix_gr. Set p11 if N even, p12 if N odd. | |
2322 | // N even: rsq = r * r | |
2323 | { .mfi | |
2324 | add table_ptr1 = 160, table_base // Point to tanl_table_p1 | |
2325 | frcpa.s1 S_hi, p0 = f1, r // S_hi for N odd | |
2326 | add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl) | |
2327 | } | |
2328 | { .mfi | |
2329 | add table_ptr2 = 400, table_base // Point to Q1_7 | |
2330 | fmpy.s1 rsq = r, r | |
2331 | nop.i 999 | |
2332 | } | |
2333 | ;; | |
2334 | ||
2335 | { .mmi | |
2336 | ldfe P1_1 = [table_ptr1], 16 | |
2337 | ;; | |
2338 | ldfe P1_2 = [table_ptr1], 16 | |
2339 | tbit.z p11, p12 = N_fix_gr, 0 | |
2340 | } | |
2341 | ;; | |
2342 | ||
2343 | ||
2344 | { .mfi | |
2345 | ldfe P1_3 = [table_ptr1], 96 | |
2346 | nop.f 999 | |
2347 | nop.i 999 | |
2348 | } | |
2349 | ;; | |
2350 | ||
2351 | { .mfi | |
2352 | (p11) ldfe P1_9 = [table_ptr1], -16 | |
2353 | (p12) fmerge.ns S_hi = S_hi, S_hi | |
2354 | nop.i 999 | |
2355 | } | |
2356 | { .mfi | |
2357 | nop.m 999 | |
2358 | (p11) fmpy.s1 r_to_the_8 = rsq, rsq | |
2359 | nop.i 999 | |
2360 | } | |
2361 | ;; | |
2362 | ||
2363 | // | |
2364 | // N even: Poly2 = P1_7 + Poly2 * rsq | |
2365 | // N odd: poly2 = Q1_5 + poly2 * rsq | |
2366 | // | |
2367 | { .mfi | |
2368 | (p11) ldfe P1_8 = [table_ptr1], -16 | |
2369 | (p11) fadd.s1 CORR = rsq, f1 | |
2370 | nop.i 999 | |
2371 | } | |
2372 | ;; | |
2373 | ||
2374 | // | |
2375 | // N even: Poly1 = P1_2 + P1_3 * rsq | |
2376 | // N odd: poly1 = 1.0 + S_hi * r | |
2377 | // 16 bits partial account for necessary (-1) | |
2378 | // | |
2379 | { .mmi | |
2380 | (p11) ldfe P1_7 = [table_ptr1], -16 | |
2381 | ;; | |
2382 | (p11) ldfe P1_6 = [table_ptr1], -16 | |
2383 | nop.i 999 | |
2384 | } | |
2385 | ;; | |
2386 | ||
2387 | // | |
2388 | // N even: Poly1 = P1_1 + Poly1 * rsq | |
2389 | // N odd: S_hi = S_hi + S_hi * poly1) 16 bits account for necessary | |
2390 | // | |
2391 | // | |
2392 | // N even: Poly2 = P1_5 + Poly2 * rsq | |
2393 | // N odd: poly2 = Q1_3 + poly2 * rsq | |
2394 | // | |
2395 | { .mfi | |
2396 | (p11) ldfe P1_5 = [table_ptr1], -16 | |
2397 | (p11) fmpy.s1 r_to_the_8 = r_to_the_8, r_to_the_8 | |
2398 | nop.i 999 | |
2399 | } | |
2400 | { .mfi | |
2401 | nop.m 999 | |
2402 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2403 | nop.i 999 | |
2404 | } | |
2405 | ;; | |
2406 | ||
2407 | // | |
2408 | // N even: Poly1 = Poly1 * rsq | |
2409 | // N odd: poly1 = 1.0 + S_hi * r 32 bits partial | |
2410 | // | |
2411 | ||
2412 | // | |
2413 | // N even: CORR = CORR * c | |
2414 | // N odd: S_hi = S_hi * poly1 + S_hi 32 bits | |
2415 | // | |
2416 | ||
2417 | // | |
2418 | // N even: Poly2 = P1_6 + Poly2 * rsq | |
2419 | // N odd: poly2 = Q1_4 + poly2 * rsq | |
2420 | // | |
2421 | ||
2422 | { .mmf | |
2423 | (p11) ldfe P1_4 = [table_ptr1], -16 | |
2424 | nop.m 999 | |
2425 | (p11) fmpy.s1 CORR = CORR, c | |
2426 | } | |
2427 | ;; | |
2428 | ||
2429 | { .mfi | |
2430 | nop.m 999 | |
2431 | (p11) fma.s1 Poly1 = P1_3, rsq, P1_2 | |
2432 | nop.i 999 ;; | |
2433 | } | |
2434 | { .mfi | |
2435 | (p12) ldfe Q1_7 = [table_ptr2], -16 | |
2436 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2437 | nop.i 999 ;; | |
2438 | } | |
2439 | { .mfi | |
2440 | (p12) ldfe Q1_6 = [table_ptr2], -16 | |
2441 | (p11) fma.s1 Poly2 = P1_9, rsq, P1_8 | |
2442 | nop.i 999 ;; | |
2443 | } | |
2444 | { .mmi | |
2445 | (p12) ldfe Q1_5 = [table_ptr2], -16 ;; | |
2446 | (p12) ldfe Q1_4 = [table_ptr2], -16 | |
2447 | nop.i 999 ;; | |
2448 | } | |
2449 | { .mfi | |
2450 | (p12) ldfe Q1_3 = [table_ptr2], -16 | |
2451 | // | |
2452 | // N even: Poly2 = P1_8 + P1_9 * rsq | |
2453 | // N odd: poly2 = Q1_6 + Q1_7 * rsq | |
2454 | // | |
2455 | (p11) fma.s1 Poly1 = Poly1, rsq, P1_1 | |
2456 | nop.i 999 ;; | |
2457 | } | |
2458 | { .mfi | |
2459 | (p12) ldfe Q1_2 = [table_ptr2], -16 | |
2460 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2461 | nop.i 999 ;; | |
2462 | } | |
2463 | { .mfi | |
2464 | (p12) ldfe Q1_1 = [table_ptr2], -16 | |
2465 | (p11) fma.s1 Poly2 = Poly2, rsq, P1_7 | |
2466 | nop.i 999 ;; | |
2467 | } | |
2468 | { .mfi | |
2469 | nop.m 999 | |
2470 | // | |
2471 | // N even: CORR = rsq + 1 | |
2472 | // N even: r_to_the_8 = rsq * rsq | |
2473 | // | |
2474 | (p11) fmpy.s1 Poly1 = Poly1, rsq | |
2475 | nop.i 999 ;; | |
2476 | } | |
2477 | { .mfi | |
2478 | nop.m 999 | |
2479 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2480 | nop.i 999 | |
2481 | } | |
2482 | { .mfi | |
2483 | nop.m 999 | |
2484 | (p12) fma.s1 poly2 = Q1_7, rsq, Q1_6 | |
2485 | nop.i 999 ;; | |
2486 | } | |
2487 | { .mfi | |
2488 | nop.m 999 | |
2489 | (p11) fma.s1 Poly2 = Poly2, rsq, P1_6 | |
2490 | nop.i 999 ;; | |
2491 | } | |
2492 | { .mfi | |
2493 | nop.m 999 | |
2494 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2495 | nop.i 999 | |
2496 | } | |
2497 | { .mfi | |
2498 | nop.m 999 | |
2499 | (p12) fma.s1 poly2 = poly2, rsq, Q1_5 | |
2500 | nop.i 999 ;; | |
2501 | } | |
2502 | { .mfi | |
2503 | nop.m 999 | |
2504 | (p11) fma.s1 Poly2= Poly2, rsq, P1_5 | |
2505 | nop.i 999 ;; | |
2506 | } | |
2507 | { .mfi | |
2508 | nop.m 999 | |
2509 | (p12) fma.s1 S_hi = S_hi, poly1, S_hi | |
2510 | nop.i 999 | |
2511 | } | |
2512 | { .mfi | |
2513 | nop.m 999 | |
2514 | (p12) fma.s1 poly2 = poly2, rsq, Q1_4 | |
2515 | nop.i 999 ;; | |
2516 | } | |
2517 | { .mfi | |
2518 | nop.m 999 | |
2519 | // | |
2520 | // N even: r_to_the_8 = r_to_the_8 * r_to_the_8 | |
2521 | // N odd: poly1 = S_hi * r + 1.0 64 bits partial | |
2522 | // | |
2523 | (p11) fma.s1 Poly2 = Poly2, rsq, P1_4 | |
2524 | nop.i 999 ;; | |
2525 | } | |
2526 | { .mfi | |
2527 | nop.m 999 | |
2528 | // | |
2529 | // N even: Poly = CORR + Poly * r | |
2530 | // N odd: P = Q1_1 + poly2 * rsq | |
2531 | // | |
2532 | (p12) fma.s1 poly1 = S_hi, r, f1 | |
2533 | nop.i 999 | |
2534 | } | |
2535 | { .mfi | |
2536 | nop.m 999 | |
2537 | (p12) fma.s1 poly2 = poly2, rsq, Q1_3 | |
2538 | nop.i 999 ;; | |
2539 | } | |
2540 | { .mfi | |
2541 | nop.m 999 | |
2542 | // | |
2543 | // N even: Poly2 = P1_4 + Poly2 * rsq | |
2544 | // N odd: poly2 = Q1_2 + poly2 * rsq | |
2545 | // | |
2546 | (p11) fma.s1 Poly = Poly2, r_to_the_8, Poly1 | |
2547 | nop.i 999 ;; | |
2548 | } | |
2549 | { .mfi | |
2550 | nop.m 999 | |
2551 | (p12) fma.s1 poly1 = S_hi, c, poly1 | |
2552 | nop.i 999 | |
2553 | } | |
2554 | { .mfi | |
2555 | nop.m 999 | |
2556 | (p12) fma.s1 poly2 = poly2, rsq, Q1_2 | |
2557 | nop.i 999 ;; | |
2558 | } | |
2559 | ||
2560 | { .mfi | |
2561 | nop.m 999 | |
2562 | // | |
2563 | // N even: Poly = Poly1 + Poly2 * r_to_the_8 | |
2564 | // N odd: S_hi = S_hi * poly1 + S_hi 64 bits | |
2565 | // | |
2566 | (p11) fma.s1 Poly = Poly, r, CORR | |
2567 | nop.i 999 ;; | |
2568 | } | |
2569 | { .mfi | |
2570 | nop.m 999 | |
2571 | // | |
2572 | // N even: Result = r + Poly (User supplied rounding mode) | |
2573 | // N odd: poly1 = S_hi * c + poly1 | |
2574 | // | |
2575 | (p12) fmpy.s1 S_lo = S_hi, poly1 | |
2576 | (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl | |
2577 | } | |
2578 | { .mfi | |
2579 | nop.m 999 | |
2580 | (p12) fma.s1 P = poly2, rsq, Q1_1 | |
2581 | nop.i 999 ;; | |
2582 | } | |
2583 | { .mfi | |
2584 | nop.m 999 | |
2585 | // | |
2586 | // N odd: poly1 = S_hi * r + 1.0 | |
2587 | // | |
2588 | // | |
2589 | // N odd: S_lo = S_hi * poly1 | |
2590 | // | |
2591 | (p14) fadd.s0 Result = Poly, r // for tanl | |
2592 | nop.i 999 | |
2593 | } | |
2594 | { .mfi | |
2595 | nop.m 999 | |
2596 | (p15) fms.s0 Result = Poly, mOne, r // for cotl | |
2597 | nop.i 999 ;; | |
2598 | } | |
2599 | ||
2600 | { .mfi | |
2601 | nop.m 999 | |
2602 | // | |
2603 | // N odd: S_lo = Q1_1 * c + S_lo | |
2604 | // | |
2605 | (p12) fma.s1 S_lo = Q1_1, c, S_lo | |
2606 | nop.i 999 | |
2607 | } | |
2608 | { .mfi | |
2609 | nop.m 999 | |
2610 | fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact | |
2611 | nop.i 999 ;; | |
2612 | } | |
2613 | { .mfi | |
2614 | nop.m 999 | |
2615 | // | |
2616 | // N odd: Result = S_lo + r * P | |
2617 | // | |
2618 | (p12) fma.s1 Result = P, r, S_lo | |
2619 | (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl | |
2620 | } | |
2621 | ||
2622 | // | |
2623 | // N odd: Result = Result + S_hi (user supplied rounding mode) | |
2624 | // | |
2625 | { .mfi | |
2626 | nop.m 999 | |
2627 | (p14) fadd.s0 Result = Result, S_hi // for tanl | |
2628 | nop.i 999 | |
2629 | } | |
2630 | { .mfb | |
2631 | nop.m 999 | |
2632 | (p15) fms.s0 Result = Result, mOne, S_hi // for cotl | |
2633 | br.ret.sptk b0 ;; // Exit |r| < 1/4 path | |
2634 | } | |
2635 | ||
2636 | ||
2637 | TANL_NORMAL_R: | |
2638 | // Here if 1/4 <= |x| < pi/4 or if |x| >= 2^63 and |r| >= 1/4 | |
2639 | // ******************************************************************* | |
2640 | // ******************************************************************* | |
2641 | // ******************************************************************* | |
2642 | // | |
2643 | // r and c have been computed. | |
2644 | // | |
2645 | { .mfi | |
2646 | nop.m 999 | |
2647 | fand B = B_mask1, r | |
2648 | nop.i 999 | |
2649 | } | |
2650 | ;; | |
2651 | ||
2652 | TANL_NORMAL_R_A: | |
2653 | // Enter here if pi/4 <= |x| < 2^63 and |r| >= 1/4 | |
2654 | // Get the 5 bits or r for the lookup. 1.xxxxx .... | |
2655 | { .mmi | |
2656 | add table_ptr1 = 416, table_base // Point to tanl_table_p2 | |
2657 | mov GR_exp_2tom65 = 0xffff - 65 // Scaling constant for B | |
2658 | extr.u lookup = sig_r, 58, 5 | |
2659 | } | |
2660 | ;; | |
2661 | ||
2662 | { .mmi | |
2663 | ldfe P2_1 = [table_ptr1], 16 | |
2664 | setf.exp TWO_TO_NEG65 = GR_exp_2tom65 // 2^-65 for scaling B if exp_r=-2 | |
2665 | add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl) | |
2666 | } | |
2667 | ;; | |
2668 | ||
2669 | .pred.rel "mutex",p11,p12 | |
2670 | // B = 2^63 * 1.xxxxx 100...0 | |
2671 | { .mfi | |
2672 | ldfe P2_2 = [table_ptr1], 16 | |
2673 | for B = B_mask2, B | |
2674 | mov table_offset = 512 // Assume table offset is 512 | |
2675 | } | |
2676 | ;; | |
2677 | ||
2678 | { .mfi | |
2679 | ldfe P2_3 = [table_ptr1], 16 | |
2680 | fmerge.s Pos_r = f1, r | |
2681 | tbit.nz p8,p9 = exp_r, 0 | |
2682 | } | |
2683 | ;; | |
2684 | ||
2685 | // Is B = 2** -2 or B= 2** -1? If 2**-1, then | |
2686 | // we want an offset of 512 for table addressing. | |
2687 | { .mii | |
2688 | add table_ptr2 = 1296, table_base // Point to tanl_table_cm2 | |
2689 | (p9) shladd table_offset = lookup, 4, table_offset | |
2690 | (p8) shladd table_offset = lookup, 4, r0 | |
2691 | } | |
2692 | ;; | |
2693 | ||
2694 | { .mmi | |
2695 | add table_ptr1 = table_ptr1, table_offset // Point to T_hi | |
2696 | add table_ptr2 = table_ptr2, table_offset // Point to C_hi | |
2697 | add table_ptr3 = 2128, table_base // Point to tanl_table_scim2 | |
2698 | } | |
2699 | ;; | |
2700 | ||
2701 | { .mmi | |
2702 | ldfd T_hi = [table_ptr1], 8 // Load T_hi | |
2703 | ;; | |
2704 | ldfd C_hi = [table_ptr2], 8 // Load C_hi | |
2705 | add table_ptr3 = table_ptr3, table_offset // Point to SC_inv | |
2706 | } | |
2707 | ;; | |
2708 | ||
2709 | // | |
2710 | // x = |r| - B | |
2711 | // | |
2712 | // Convert B so it has the same exponent as Pos_r before subtracting | |
2713 | { .mfi | |
2714 | ldfs T_lo = [table_ptr1] // Load T_lo | |
2715 | (p9) fnma.s1 x = B, FR_2tom64, Pos_r | |
2716 | nop.i 999 | |
2717 | } | |
2718 | { .mfi | |
2719 | nop.m 999 | |
2720 | (p8) fnma.s1 x = B, TWO_TO_NEG65, Pos_r | |
2721 | nop.i 999 | |
2722 | } | |
2723 | ;; | |
2724 | ||
2725 | { .mfi | |
2726 | ldfs C_lo = [table_ptr2] // Load C_lo | |
2727 | nop.f 999 | |
2728 | nop.i 999 | |
2729 | } | |
2730 | ;; | |
2731 | ||
2732 | { .mfi | |
2733 | ldfe SC_inv = [table_ptr3] // Load SC_inv | |
2734 | fmerge.s sgn_r = r, f1 | |
2735 | tbit.z p11, p12 = N_fix_gr, 0 // p11 if N even, p12 if odd | |
2736 | ||
2737 | } | |
2738 | ;; | |
2739 | ||
2740 | // | |
2741 | // xsq = x * x | |
2742 | // N even: Tx = T_hi * x | |
2743 | // | |
2744 | // N even: Tx1 = Tx + 1 | |
2745 | // N odd: Cx1 = 1 - Cx | |
2746 | // | |
2747 | ||
2748 | { .mfi | |
2749 | nop.m 999 | |
2750 | fmpy.s1 xsq = x, x | |
2751 | nop.i 999 | |
2752 | } | |
2753 | { .mfi | |
2754 | nop.m 999 | |
2755 | (p11) fmpy.s1 Tx = T_hi, x | |
2756 | nop.i 999 | |
2757 | } | |
2758 | ;; | |
2759 | ||
2760 | // | |
2761 | // N odd: Cx = C_hi * x | |
2762 | // | |
2763 | { .mfi | |
2764 | nop.m 999 | |
2765 | (p12) fmpy.s1 Cx = C_hi, x | |
2766 | nop.i 999 | |
2767 | } | |
2768 | ;; | |
2769 | // | |
2770 | // N even and odd: P = P2_3 + P2_2 * xsq | |
2771 | // | |
2772 | { .mfi | |
2773 | nop.m 999 | |
2774 | fma.s1 P = P2_3, xsq, P2_2 | |
2775 | nop.i 999 | |
2776 | } | |
2777 | { .mfi | |
2778 | nop.m 999 | |
2779 | (p11) fadd.s1 Tx1 = Tx, f1 | |
2780 | nop.i 999 ;; | |
2781 | } | |
2782 | { .mfi | |
2783 | nop.m 999 | |
2784 | // | |
2785 | // N even: D = C_hi - tanx | |
2786 | // N odd: D = T_hi + tanx | |
2787 | // | |
2788 | (p11) fmpy.s1 CORR = SC_inv, T_hi | |
2789 | nop.i 999 | |
2790 | } | |
2791 | { .mfi | |
2792 | nop.m 999 | |
2793 | fmpy.s1 Sx = SC_inv, x | |
2794 | nop.i 999 ;; | |
2795 | } | |
2796 | { .mfi | |
2797 | nop.m 999 | |
2798 | (p12) fmpy.s1 CORR = SC_inv, C_hi | |
2799 | nop.i 999 ;; | |
2800 | } | |
2801 | { .mfi | |
2802 | nop.m 999 | |
2803 | (p12) fsub.s1 V_hi = f1, Cx | |
2804 | nop.i 999 ;; | |
2805 | } | |
2806 | { .mfi | |
2807 | nop.m 999 | |
2808 | fma.s1 P = P, xsq, P2_1 | |
2809 | nop.i 999 | |
2810 | } | |
2811 | { .mfi | |
2812 | nop.m 999 | |
2813 | // | |
2814 | // N even and odd: P = P2_1 + P * xsq | |
2815 | // | |
2816 | (p11) fma.s1 V_hi = Tx, Tx1, f1 | |
2817 | nop.i 999 ;; | |
2818 | } | |
2819 | { .mfi | |
2820 | nop.m 999 | |
2821 | // | |
2822 | // N even: Result = sgn_r * tail + T_hi (user rounding mode for C1) | |
2823 | // N odd: Result = sgn_r * tail + C_hi (user rounding mode for C1) | |
2824 | // | |
2825 | fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact | |
2826 | nop.i 999 ;; | |
2827 | } | |
2828 | { .mfi | |
2829 | nop.m 999 | |
2830 | fmpy.s1 CORR = CORR, c | |
2831 | nop.i 999 ;; | |
2832 | } | |
2833 | { .mfi | |
2834 | nop.m 999 | |
2835 | (p12) fnma.s1 V_hi = Cx,V_hi,f1 | |
2836 | nop.i 999 ;; | |
2837 | } | |
2838 | { .mfi | |
2839 | nop.m 999 | |
2840 | // | |
2841 | // N even: V_hi = Tx * Tx1 + 1 | |
2842 | // N odd: Cx1 = 1 - Cx * Cx1 | |
2843 | // | |
2844 | fmpy.s1 P = P, xsq | |
2845 | nop.i 999 | |
2846 | } | |
2847 | { .mfi | |
2848 | nop.m 999 | |
2849 | // | |
2850 | // N even and odd: P = P * xsq | |
2851 | // | |
2852 | (p11) fmpy.s1 V_hi = V_hi, T_hi | |
2853 | nop.i 999 ;; | |
2854 | } | |
2855 | { .mfi | |
2856 | nop.m 999 | |
2857 | // | |
2858 | // N even and odd: tail = P * tail + V_lo | |
2859 | // | |
2860 | (p11) fmpy.s1 T_hi = sgn_r, T_hi | |
2861 | nop.i 999 ;; | |
2862 | } | |
2863 | { .mfi | |
2864 | nop.m 999 | |
2865 | fmpy.s1 CORR = CORR, sgn_r | |
2866 | nop.i 999 ;; | |
2867 | } | |
2868 | { .mfi | |
2869 | nop.m 999 | |
2870 | (p12) fmpy.s1 V_hi = V_hi,C_hi | |
2871 | nop.i 999 ;; | |
2872 | } | |
2873 | { .mfi | |
2874 | nop.m 999 | |
2875 | // | |
2876 | // N even: V_hi = T_hi * V_hi | |
2877 | // N odd: V_hi = C_hi * V_hi | |
2878 | // | |
2879 | fma.s1 tanx = P, x, x | |
2880 | nop.i 999 | |
2881 | } | |
2882 | { .mfi | |
2883 | nop.m 999 | |
2884 | (p12) fnmpy.s1 C_hi = sgn_r, C_hi | |
2885 | nop.i 999 ;; | |
2886 | } | |
2887 | { .mfi | |
2888 | nop.m 999 | |
2889 | // | |
2890 | // N even: V_lo = 1 - V_hi + C_hi | |
2891 | // N odd: V_lo = 1 - V_hi + T_hi | |
2892 | // | |
2893 | (p11) fadd.s1 CORR = CORR, T_lo | |
2894 | nop.i 999 | |
2895 | } | |
2896 | { .mfi | |
2897 | nop.m 999 | |
2898 | (p12) fsub.s1 CORR = CORR, C_lo | |
2899 | nop.i 999 ;; | |
2900 | } | |
2901 | { .mfi | |
2902 | nop.m 999 | |
2903 | // | |
2904 | // N even and odd: tanx = x + x * P | |
2905 | // N even and odd: Sx = SC_inv * x | |
2906 | // | |
2907 | (p11) fsub.s1 D = C_hi, tanx | |
2908 | nop.i 999 | |
2909 | } | |
2910 | { .mfi | |
2911 | nop.m 999 | |
2912 | (p12) fadd.s1 D = T_hi, tanx | |
2913 | nop.i 999 ;; | |
2914 | } | |
2915 | { .mfi | |
2916 | nop.m 999 | |
2917 | // | |
2918 | // N odd: CORR = SC_inv * C_hi | |
2919 | // N even: CORR = SC_inv * T_hi | |
2920 | // | |
2921 | fnma.s1 D = V_hi, D, f1 | |
2922 | nop.i 999 ;; | |
2923 | } | |
2924 | { .mfi | |
2925 | nop.m 999 | |
2926 | // | |
2927 | // N even and odd: D = 1 - V_hi * D | |
2928 | // N even and odd: CORR = CORR * c | |
2929 | // | |
2930 | fma.s1 V_hi = V_hi, D, V_hi | |
2931 | nop.i 999 ;; | |
2932 | } | |
2933 | { .mfi | |
2934 | nop.m 999 | |
2935 | // | |
2936 | // N even and odd: V_hi = V_hi + V_hi * D | |
2937 | // N even and odd: CORR = sgn_r * CORR | |
2938 | // | |
2939 | (p11) fnma.s1 V_lo = V_hi, C_hi, f1 | |
2940 | nop.i 999 | |
2941 | } | |
2942 | { .mfi | |
2943 | nop.m 999 | |
2944 | (p12) fnma.s1 V_lo = V_hi, T_hi, f1 | |
2945 | nop.i 999 ;; | |
2946 | } | |
2947 | { .mfi | |
2948 | nop.m 999 | |
2949 | // | |
2950 | // N even: CORR = COOR + T_lo | |
2951 | // N odd: CORR = CORR - C_lo | |
2952 | // | |
2953 | (p11) fma.s1 V_lo = tanx, V_hi, V_lo | |
2954 | tbit.nz p15, p0 = cot_flag, 0 // p15=1 if we compute cotl | |
2955 | } | |
2956 | { .mfi | |
2957 | nop.m 999 | |
2958 | (p12) fnma.s1 V_lo = tanx, V_hi, V_lo | |
2959 | nop.i 999 ;; | |
2960 | } | |
2961 | ||
2962 | { .mfi | |
2963 | nop.m 999 | |
2964 | (p15) fms.s1 T_hi = f0, f0, T_hi // to correct result's sign for cotl | |
2965 | nop.i 999 | |
2966 | } | |
2967 | { .mfi | |
2968 | nop.m 999 | |
2969 | (p15) fms.s1 C_hi = f0, f0, C_hi // to correct result's sign for cotl | |
2970 | nop.i 999 | |
2971 | };; | |
2972 | ||
2973 | { .mfi | |
2974 | nop.m 999 | |
2975 | (p15) fms.s1 sgn_r = f0, f0, sgn_r // to correct result's sign for cotl | |
2976 | nop.i 999 | |
2977 | };; | |
2978 | ||
2979 | { .mfi | |
2980 | nop.m 999 | |
2981 | // | |
2982 | // N even: V_lo = V_lo + V_hi * tanx | |
2983 | // N odd: V_lo = V_lo - V_hi * tanx | |
2984 | // | |
2985 | (p11) fnma.s1 V_lo = C_lo, V_hi, V_lo | |
2986 | nop.i 999 | |
2987 | } | |
2988 | { .mfi | |
2989 | nop.m 999 | |
2990 | (p12) fnma.s1 V_lo = T_lo, V_hi, V_lo | |
2991 | nop.i 999 ;; | |
2992 | } | |
2993 | { .mfi | |
2994 | nop.m 999 | |
2995 | // | |
2996 | // N even: V_lo = V_lo - V_hi * C_lo | |
2997 | // N odd: V_lo = V_lo - V_hi * T_lo | |
2998 | // | |
2999 | fmpy.s1 V_lo = V_hi, V_lo | |
3000 | nop.i 999 ;; | |
3001 | } | |
3002 | { .mfi | |
3003 | nop.m 999 | |
3004 | // | |
3005 | // N even and odd: V_lo = V_lo * V_hi | |
3006 | // | |
3007 | fadd.s1 tail = V_hi, V_lo | |
3008 | nop.i 999 ;; | |
3009 | } | |
3010 | { .mfi | |
3011 | nop.m 999 | |
3012 | // | |
3013 | // N even and odd: tail = V_hi + V_lo | |
3014 | // | |
3015 | fma.s1 tail = tail, P, V_lo | |
3016 | nop.i 999 ;; | |
3017 | } | |
3018 | { .mfi | |
3019 | nop.m 999 | |
3020 | // | |
3021 | // N even: T_hi = sgn_r * T_hi | |
3022 | // N odd : C_hi = -sgn_r * C_hi | |
3023 | // | |
3024 | fma.s1 tail = tail, Sx, CORR | |
3025 | nop.i 999 ;; | |
3026 | } | |
3027 | { .mfi | |
3028 | nop.m 999 | |
3029 | // | |
3030 | // N even and odd: tail = Sx * tail + CORR | |
3031 | // | |
3032 | fma.s1 tail = V_hi, Sx, tail | |
3033 | nop.i 999 ;; | |
3034 | } | |
3035 | { .mfi | |
3036 | nop.m 999 | |
3037 | // | |
3038 | // N even an odd: tail = Sx * V_hi + tail | |
3039 | // | |
3040 | (p11) fma.s0 Result = sgn_r, tail, T_hi | |
3041 | nop.i 999 | |
3042 | } | |
3043 | { .mfb | |
3044 | nop.m 999 | |
3045 | (p12) fma.s0 Result = sgn_r, tail, C_hi | |
3046 | br.ret.sptk b0 ;; // Exit for 1/4 <= |r| < pi/4 | |
3047 | } | |
3048 | ||
3049 | TANL_DENORMAL: | |
3050 | // Here if x denormal | |
3051 | { .mfb | |
3052 | getf.exp GR_signexp_x = Norm_Arg // Get sign and exponent of x | |
3053 | nop.f 999 | |
3054 | br.cond.sptk TANL_COMMON // Return to common code | |
3055 | } | |
3056 | ;; | |
3057 | ||
3058 | ||
3059 | TANL_SPECIAL: | |
3060 | TANL_UNSUPPORTED: | |
3061 | // | |
3062 | // Code for NaNs, Unsupporteds, Infs, or +/- zero ? | |
3063 | // Invalid raised for Infs and SNaNs. | |
3064 | // | |
3065 | ||
3066 | { .mfi | |
3067 | nop.m 999 | |
3068 | fmerge.s f10 = f8, f8 // Save input for error call | |
3069 | tbit.nz p6, p7 = cot_flag, 0 // p6=1 if we compute cotl | |
3070 | } | |
3071 | ;; | |
3072 | ||
3073 | { .mfi | |
3074 | nop.m 999 | |
3075 | (p6) fclass.m p6, p7 = f8, 0x7 // Test for zero (cotl only) | |
3076 | nop.i 999 | |
3077 | } | |
3078 | ;; | |
3079 | ||
3080 | .pred.rel "mutex", p6, p7 | |
3081 | { .mfi | |
3082 | (p6) mov GR_Parameter_Tag = 225 // (cotl) | |
3083 | (p6) frcpa.s0 f8, p0 = f1, f8 // cotl(+-0) = +-Inf | |
3084 | nop.i 999 | |
3085 | } | |
3086 | { .mfb | |
3087 | nop.m 999 | |
3088 | (p7) fmpy.s0 f8 = f8, f0 | |
3089 | (p7) br.ret.sptk b0 | |
3090 | } | |
3091 | ;; | |
3092 | ||
3093 | GLOBAL_IEEE754_END(tanl) | |
3094 | ||
3095 | ||
3096 | LOCAL_LIBM_ENTRY(__libm_error_region) | |
3097 | .prologue | |
3098 | ||
3099 | // (1) | |
3100 | { .mfi | |
3101 | add GR_Parameter_Y=-32,sp // Parameter 2 value | |
3102 | nop.f 0 | |
3103 | .save ar.pfs,GR_SAVE_PFS | |
3104 | mov GR_SAVE_PFS=ar.pfs // Save ar.pfs | |
3105 | } | |
3106 | { .mfi | |
3107 | .fframe 64 | |
3108 | add sp=-64,sp // Create new stack | |
3109 | nop.f 0 | |
3110 | mov GR_SAVE_GP=gp // Save gp | |
3111 | };; | |
3112 | ||
3113 | // (2) | |
3114 | { .mmi | |
3115 | stfe [GR_Parameter_Y] = f1,16 // STORE Parameter 2 on stack | |
3116 | add GR_Parameter_X = 16,sp // Parameter 1 address | |
3117 | .save b0, GR_SAVE_B0 | |
3118 | mov GR_SAVE_B0=b0 // Save b0 | |
3119 | };; | |
3120 | ||
3121 | .body | |
3122 | // (3) | |
3123 | { .mib | |
3124 | stfe [GR_Parameter_X] = f10 // STORE Parameter 1 on stack | |
3125 | add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address | |
3126 | nop.b 0 | |
3127 | } | |
3128 | { .mib | |
3129 | stfe [GR_Parameter_Y] = f8 // STORE Parameter 3 on stack | |
3130 | add GR_Parameter_Y = -16,GR_Parameter_Y | |
3131 | br.call.sptk b0=__libm_error_support# // Call error handling function | |
3132 | };; | |
3133 | { .mmi | |
3134 | nop.m 0 | |
3135 | nop.m 0 | |
3136 | add GR_Parameter_RESULT = 48,sp | |
3137 | };; | |
3138 | ||
3139 | // (4) | |
3140 | { .mmi | |
3141 | ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack | |
3142 | .restore sp | |
3143 | add sp = 64,sp // Restore stack pointer | |
3144 | mov b0 = GR_SAVE_B0 // Restore return address | |
3145 | };; | |
3146 | { .mib | |
3147 | mov gp = GR_SAVE_GP // Restore gp | |
3148 | mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs | |
3149 | br.ret.sptk b0 // Return | |
3150 | };; | |
3151 | ||
3152 | LOCAL_LIBM_END(__libm_error_region) | |
3153 | ||
3154 | .type __libm_error_support#,@function | |
3155 | .global __libm_error_support# | |
3156 | ||
3157 | ||
3158 | // ******************************************************************* | |
3159 | // ******************************************************************* | |
3160 | // ******************************************************************* | |
3161 | // | |
3162 | // Special Code to handle very large argument case. | |
3163 | // Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63 | |
3164 | // The interface is custom: | |
3165 | // On input: | |
3166 | // (Arg or x) is in f8 | |
3167 | // On output: | |
3168 | // r is in f8 | |
3169 | // c is in f9 | |
3170 | // N is in r8 | |
3171 | // We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We | |
3172 | // use this to eliminate save/restore of key fp registers in this calling | |
3173 | // function. | |
3174 | // | |
3175 | // ******************************************************************* | |
3176 | // ******************************************************************* | |
3177 | // ******************************************************************* | |
3178 | ||
3179 | LOCAL_LIBM_ENTRY(__libm_callout) | |
3180 | TANL_ARG_TOO_LARGE: | |
3181 | .prologue | |
3182 | { .mfi | |
3183 | add table_ptr2 = 144, table_base // Point to 2^-2 | |
3184 | nop.f 999 | |
3185 | .save ar.pfs,GR_SAVE_PFS | |
3186 | mov GR_SAVE_PFS=ar.pfs // Save ar.pfs | |
3187 | } | |
3188 | ;; | |
3189 | ||
3190 | // Load 2^-2, -2^-2 | |
3191 | { .mmi | |
3192 | ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] | |
3193 | setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r | |
3194 | .save b0, GR_SAVE_B0 | |
3195 | mov GR_SAVE_B0=b0 // Save b0 | |
3196 | };; | |
3197 | ||
3198 | .body | |
3199 | // | |
3200 | // Call argument reduction with x in f8 | |
3201 | // Returns with N in r8, r in f8, c in f9 | |
3202 | // Assumes f71-127 are preserved across the call | |
3203 | // | |
3204 | { .mib | |
3205 | setf.sig B_mask2 = bmask2 // Form mask to form B from r | |
3206 | mov GR_SAVE_GP=gp // Save gp | |
3207 | br.call.sptk b0=__libm_pi_by_2_reduce# | |
3208 | } | |
3209 | ;; | |
3210 | ||
3211 | // | |
3212 | // Is |r| < 2**(-2) | |
3213 | // | |
3214 | { .mfi | |
3215 | getf.sig sig_r = r // Extract significand of r | |
3216 | fcmp.lt.s1 p6, p0 = r, TWO_TO_NEG2 | |
3217 | mov gp = GR_SAVE_GP // Restore gp | |
3218 | } | |
3219 | ;; | |
3220 | ||
3221 | { .mfi | |
3222 | getf.exp exp_r = r // Extract signexp of r | |
3223 | nop.f 999 | |
3224 | mov b0 = GR_SAVE_B0 // Restore return address | |
3225 | } | |
3226 | ;; | |
3227 | ||
3228 | // | |
3229 | // Get N_fix_gr | |
3230 | // | |
3231 | { .mfi | |
3232 | mov N_fix_gr = r8 | |
3233 | (p6) fcmp.gt.unc.s1 p6, p0 = r, NEGTWO_TO_NEG2 | |
3234 | mov ar.pfs = GR_SAVE_PFS // Restore pfs | |
3235 | } | |
3236 | ;; | |
3237 | ||
3238 | { .mbb | |
3239 | nop.m 999 | |
3240 | (p6) br.cond.spnt TANL_SMALL_R // Branch if |r| < 1/4 | |
3241 | br.cond.sptk TANL_NORMAL_R // Branch if 1/4 <= |r| < pi/4 | |
3242 | } | |
3243 | ;; | |
3244 | ||
3245 | LOCAL_LIBM_END(__libm_callout) | |
3246 | ||
3247 | .type __libm_pi_by_2_reduce#,@function | |
3248 | .global __libm_pi_by_2_reduce# |