| blob | 28d44c1850fda288e50e748a644e0c2bd4a1b7d7 |
1 .file "atanl.s"
3 // Copyright (C) 2000, 2001, Intel Corporation
4 // All rights reserved.
5 //
6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
7 // and Ping Tak Peter Tang of the Computational Software Lab, 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 //
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
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
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
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 //
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://developer.intel.com/opensource.
39 //
40 //
41 // *********************************************************************
42 //
43 // History
44 // 2/02/00 (hand-optimized)
45 // 4/04/00 Unwind support added
46 // 8/15/00 Bundle added after call to __libm_error_support to properly
47 // set [the previously overwritten] GR_Parameter_RESULT.
48 //
49 // *********************************************************************
50 //
51 // Function: atanl(x) = inverse tangent(x), for double extended x values
52 // Function: atan2l(y,x) = atan(y/x), for double extended x values
53 //
54 // *********************************************************************
55 //
56 // Resources Used:
57 //
58 // Floating-Point Registers: f8 (Input and Return Value)
59 // f9-f15
60 // f32-f79
61 //
62 // General Purpose Registers:
63 // r32-r48
64 // r49,r50,r51,r52 (Arguments to error support for 0,0 case)
65 //
66 // Predicate Registers: p6-p15
67 //
68 // *********************************************************************
69 //
70 // IEEE Special Conditions:
71 //
72 // Denormal fault raised on denormal inputs
73 // Underflow exceptions may occur
74 // Special error handling for the y=0 and x=0 case
75 // Inexact raised when appropriate by algorithm
76 //
77 // atanl(SNaN) = QNaN
78 // atanl(QNaN) = QNaN
79 // atanl(+/-0) = +/- 0
80 // atanl(+/-Inf) = +/-pi/2
81 //
82 // atan2l(Any NaN for x or y) = QNaN
83 // atan2l(+/-0,x) = +/-0 for x > 0
84 // atan2l(+/-0,x) = +/-pi for x < 0
85 // atan2l(+/-0,+0) = +/-0
86 // atan2l(+/-0,-0) = +/-pi
87 // atan2l(y,+/-0) = pi/2 y > 0
88 // atan2l(y,+/-0) = -pi/2 y < 0
89 // atan2l(+/-y, Inf) = +/-0 for finite y > 0
90 // atan2l(+/-Inf, x) = +/-pi/2 for finite x
91 // atan2l(+/-y, -Inf) = +/-pi for finite y > 0
92 // atan2l(+/-Inf, Inf) = +/-pi/4
93 // atan2l(+/-Inf, -Inf) = +/-3pi/4
94 //
95 // *********************************************************************
96 //
97 // Mathematical Description
98 // ---------------------------
99 //
100 // The function ATANL( Arg_Y, Arg_X ) returns the "argument"
101 // or the "phase" of the complex number
102 //
103 // Arg_X + i Arg_Y
104 //
105 // or equivalently, the angle in radians from the positive
106 // x-axis to the line joining the origin and the point
107 // (Arg_X,Arg_Y)
108 //
109 //
110 // (Arg_X, Arg_Y) x
111 // \
112 // \
113 // \
114 // \
115 // \ angle between is ATANL(Arg_Y,Arg_X)
120 // \
121 // ------------------> X-axis
123 // Origin
124 //
125 // Moreover, this angle is reported in the range [-pi,pi] thus
126 //
127 // -pi <= ATANL( Arg_Y, Arg_X ) <= pi.
128 //
129 // From the geometry, it is easy to define ATANL when one of
130 // Arg_X or Arg_Y is +-0 or +-inf:
131 //
132 //
133 // \ Y |
134 // X \ | +0 | -0 | +inf | -inf | finite non-zero
135 // \ | | | | |
136 // ______________________________________________________
137 // | | | |
138 // +-0 | Invalid/ | pi/2 | -pi/2 | sign(Y)*pi/2
139 // | qNaN | | |
140 // --------------------------------------------------------
141 // | | | | |
142 // +inf | +0 | -0 | pi/4 | -pi/4 | sign(Y)*0
143 // --------------------------------------------------------
144 // | | | | |
145 // -inf | +pi | -pi | 3pi/4 | -3pi/4 | sign(Y)*pi
146 // --------------------------------------------------------
147 // finite | X>0? | pi/2 | -pi/2 | normal case
148 // non-zero| sign(Y)*0: | | |
149 // | sign(Y)*pi | | |
150 //
151 //
152 // One must take note that ATANL is NOT the arctangent of the
153 // value Arg_Y/Arg_X; but rather ATANL and arctan are related
154 // in a slightly more complicated way as follows:
155 //
156 // Let U := max(|Arg_X|, |Arg_Y|); V := min(|Arg_X|, |Arg_Y|);
157 // sign_X be the sign bit of Arg_X, i.e., sign_X is 0 or 1;
158 // s_X be the sign of Arg_X, i.e., s_X = (-1)^sign_X;
159 //
160 // sign_Y be the sign bit of Arg_Y, i.e., sign_Y is 0 or 1;
161 // s_Y be the sign of Arg_Y, i.e., s_Y = (-1)^sign_Y;
162 //
163 // swap be 0 if |Arg_X| >= |Arg_Y| and 1 otherwise.
164 //
165 // Then, ATANL(Arg_Y, Arg_X) =
166 //
167 // / arctan(V/U) \ sign_X = 0 & swap = 0
168 // | pi/2 - arctan(V/U) | sign_X = 0 & swap = 1
169 // s_Y * | |
170 // | pi - arctan(V/U) | sign_X = 1 & swap = 0
171 // \ pi/2 + arctan(V/U) / sign_X = 1 & swap = 1
172 //
173 //
174 // This relationship also suggest that the algorithm's major
175 // task is to calculate arctan(V/U) for 0 < V <= U; and the
176 // final Result is given by
177 //
178 // s_Y * { (P_hi + P_lo) + sigma * arctan(V/U) }
179 //
180 // where
181 //
182 // (P_hi,P_lo) represents M(sign_X,swap)*(pi/2) accurately
183 //
184 // M(sign_X,swap) = 0 for sign_X = 0 and swap = 0
185 // 1 for swap = 1
186 // 2 for sign_X = 1 and swap = 0
187 //
188 // and
189 //
190 // sigma = { (sign_X XOR swap) : -1.0 : 1.0 }
191 //
192 // = (-1) ^ ( sign_X XOR swap )
193 //
194 // Both (P_hi,P_lo) and sigma can be stored in a table and fetched
195 // using (sign_X,swap) as an index. (P_hi, P_lo) can be stored as a
196 // double-precision, and single-precision pair; and sigma can
197 // obviously be just a single-precision number.
198 //
199 // In the algorithm we propose, arctan(V/U) is calculated to high accuracy
200 // as A_hi + A_lo. Consequently, the Result ATANL( Arg_Y, Arg_X ) is
201 // given by
202 //
203 // s_Y*P_hi + s_Y*sigma*A_hi + s_Y*(sigma*A_lo + P_lo)
204 //
205 // We now discuss the calculation of arctan(V/U) for 0 < V <= U.
206 //
207 // For (V/U) < 2^(-3), we use a simple polynomial of the form
208 //
209 // z + z^3*(P_1 + z^2*(P_2 + z^2*(P_3 + ... + P_8)))
210 //
211 // where z = V/U.
212 //
213 // For the sake of accuracy, the first term "z" must approximate V/U to
214 // extra precision. For z^3 and higher power, a working precision
215 // approximation to V/U suffices. Thus, we obtain:
216 //
217 // z_hi + z_lo = V/U to extra precision and
218 // z = V/U to working precision
219 //
220 // The value arctan(V/U) is delivered as two pieces (A_hi, A_lo)
221 //
222 // (A_hi,A_lo) = (z_hi, z^3*(P_1 + ... + P_8) + z_lo).
223 //
224 //
225 // For 2^(-3) <= (V/U) <= 1, we use a table-driven approach.
226 // Consider
227 //
228 // (V/U) = 2^k * 1.b_1 b_2 .... b_63 b_64 b_65 ....
229 //
230 // Define
231 //
232 // z_hi = 2^k * 1.b_1 b_2 b_3 b_4 1
233 //
234 // then
235 // / \
236 // | (V/U) - z_hi |
238 // arctan(V/U) = arctan(z_hi) + acrtan| -------------- |
239 // | 1 + (V/U)*z_hi |
240 // \ /
241 //
242 // / \
243 // | V - z_hi*U |
245 // = arctan(z_hi) + acrtan| -------------- |
246 // | U + V*z_hi |
247 // \ /
248 //
249 // = arctan(z_hi) + acrtan( V' / U' )
250 //
251 //
252 // where
253 //
254 // V' = V - U*z_hi; U' = U + V*z_hi.
255 //
256 // Let
257 //
258 // w_hi + w_lo = V'/U' to extra precision and
259 // w = V'/U' to working precision
260 //
261 // then we can approximate arctan(V'/U') by
262 //
263 // arctan(V'/U') = w_hi + w_lo
264 // + w^3*(Q_1 + w^2*(Q_2 + w^2*(Q_3 + w^2*Q_4)))
265 //
266 // = w_hi + w_lo + poly
267 //
268 // Finally, arctan(z_hi) is calculated beforehand and stored in a table
269 // as Tbl_hi, Tbl_lo. Thus,
270 //
271 // (A_hi, A_lo) = (Tbl_hi, w_hi+(poly+(w_lo+Tbl_lo)))
272 //
273 // This completes the mathematical description.
274 //
275 //
276 // Algorithm
277 // -------------
278 //
279 // Step 0. Check for unsupported format.
280 //
281 // If
282 // ( expo(Arg_X) not zero AND msb(Arg_X) = 0 ) OR
283 // ( expo(Arg_Y) not zero AND msb(Arg_Y) = 0 )
284 //
285 // then one of the arguments is unsupported. Generate an
286 // invalid and return qNaN.
287 //
288 // Step 1. Initialize
289 //
290 // Normalize Arg_X and Arg_Y and set the following
291 //
292 // sign_X := sign_bit(Arg_X)
293 // s_Y := (sign_bit(Arg_Y)==0? 1.0 : -1.0)
294 // swap := (|Arg_X| >= |Arg_Y|? 0 : 1 )
295 // U := max( |Arg_X|, |Arg_Y| )
296 // V := min( |Arg_X|, |Arg_Y| )
297 //
298 // execute: frcap E, pred, V, U
299 // If pred is 0, go to Step 5 for special cases handling.
300 //
301 // Step 2. Decide on branch.
302 //
303 // Q := E * V
304 // If Q < 2^(-3) go to Step 4 for simple polynomial case.
305 //
306 // Step 3. Table-driven algorithm.
307 //
308 // Q is represented as
309 //
310 // 2^(-k) * 1.b_1 b_2 b_3 ... b_63; k = 0,-1,-2,-3
311 //
312 // and that if k = 0, b_1 = b_2 = b_3 = b_4 = 0.
313 //
314 // Define
315 //
316 // z_hi := 2^(-k) * 1.b_1 b_2 b_3 b_4 1
317 //
318 // (note that there are 49 possible values of z_hi).
319 //
320 // ...We now calculate V' and U'. While V' is representable
321 // ...as a 64-bit number because of cancellation, U' is
322 // ...not in general a 64-bit number. Obtaining U' accurately
323 // ...requires two working precision numbers
324 //
325 // U_prime_hi := U + V * z_hi ...WP approx. to U'
326 // U_prime_lo := ( U - U_prime_hi ) + V*z_hi ...observe order
327 // V_prime := V - U * z_hi ...this is exact
328 //
329 // C_hi := frcpa (1.0, U_prime_hi) ...C_hi approx 1/U'_hi
330 //
331 // loop 3 times
332 // C_hi := C_hi + C_hi*(1.0 - C_hi*U_prime_hi)
333 //
334 // ...at this point C_hi is (1/U_prime_hi) to roughly 64 bits
335 //
336 // w_hi := V_prime * C_hi ...w_hi is V_prime/U_prime to
337 // ...roughly working precision
338 //
339 // ...note that we want w_hi + w_lo to approximate
340 // ...V_prime/(U_prime_hi + U_prime_lo) to extra precision
341 // ...but for now, w_hi is good enough for the polynomial
342 // ...calculation.
343 //
344 // wsq := w_hi*w_hi
345 // poly := w_hi*wsq*(Q_1 + wsq*(Q_2 + wsq*(Q_3 + wsq*Q_4)))
346 //
347 // Fetch
348 // (Tbl_hi, Tbl_lo) = atan(z_hi) indexed by (k,b_1,b_2,b_3,b_4)
349 // ...Tbl_hi is a double-precision number
350 // ...Tbl_lo is a single-precision number
351 //
352 // (P_hi, P_lo) := M(sign_X,swap)*(Pi_by_2_hi, Pi_by_2_lo)
353 // ...as discussed previous. Again; the implementation can
354 // ...chose to fetch P_hi and P_lo from a table indexed by
355 // ...(sign_X, swap).
356 // ...P_hi is a double-precision number;
357 // ...P_lo is a single-precision number.
358 //
359 // ...calculate w_lo so that w_hi + w_lo is V'/U' accurately
360 // w_lo := ((V_prime - w_hi*U_prime_hi) -
361 // w_hi*U_prime_lo) * C_hi ...observe order
362 //
363 //
364 // ...Ready to deliver arctan(V'/U') as A_hi, A_lo
365 // A_hi := Tbl_hi
366 // A_lo := w_hi + (poly + (Tbl_lo + w_lo)) ...observe order
367 //
368 // ...Deliver final Result
369 // ...s_Y*P_hi + s_Y*sigma*A_hi + s_Y*(sigma*A_lo + P_lo)
370 //
371 // sigma := ( (sign_X XOR swap) ? -1.0 : 1.0 )
372 // ...sigma can be obtained by a table lookup using
373 // ...(sign_X,swap) as index and stored as single precision
374 // ...sigma should be calculated earlier
375 //
376 // P_hi := s_Y*P_hi
377 // A_hi := s_Y*A_hi
378 //
379 // Res_hi := P_hi + sigma*A_hi ...this is exact because
380 // ...both P_hi and Tbl_hi
381 // ...are double-precision
382 // ...and |Tbl_hi| > 2^(-4)
383 // ...P_hi is either 0 or
384 // ...between (1,4)
385 //
386 // Res_lo := sigma*A_lo + P_lo
387 //
388 // Return Res_hi + s_Y*Res_lo in user-defined rounding control
389 //
390 // Step 4. Simple polynomial case.
391 //
392 // ...E and Q are inherited from Step 2.
393 //
394 // A_hi := Q ...Q is inherited from Step 2 Q approx V/U
395 //
396 // loop 3 times
397 // E := E + E2(1.0 - E*U1
398 // ...at this point E approximates 1/U to roughly working precision
399 //
400 // z := V * E ...z approximates V/U to roughly working precision
401 // zsq := z * z
402 // z8 := zsq * zsq; z8 := z8 * z8
403 //
404 // poly1 := P_4 + zsq*(P_5 + zsq*(P_6 + zsq*(P_7 + zsq*P_8)))
405 // poly2 := zsq*(P_1 + zsq*(P_2 + zsq*P_3))
406 //
407 // poly := poly1 + z8*poly2
408 //
409 // z_lo := (V - A_hi*U)*E
410 //
411 // A_lo := z*poly + z_lo
412 // ...A_hi, A_lo approximate arctan(V/U) accurately
413 //
414 // (P_hi, P_lo) := M(sign_X,swap)*(Pi_by_2_hi, Pi_by_2_lo)
415 // ...one can store the M(sign_X,swap) as single precision
416 // ...values
417 //
418 // ...Deliver final Result
419 // ...s_Y*P_hi + s_Y*sigma*A_hi + s_Y*(sigma*A_lo + P_lo)
420 //
421 // sigma := ( (sign_X XOR swap) ? -1.0 : 1.0 )
422 // ...sigma can be obtained by a table lookup using
423 // ...(sign_X,swap) as index and stored as single precision
424 // ...sigma should be calculated earlier
425 //
426 // P_hi := s_Y*P_hi
427 // A_hi := s_Y*A_hi
428 //
429 // Res_hi := P_hi + sigma*A_hi ...need to compute
430 // ...P_hi + sigma*A_hi
431 // ...exactly
432 //
433 // tmp := (P_hi - Res_hi) + sigma*A_hi
434 //
435 // Res_lo := s_Y*(sigma*A_lo + P_lo) + tmp
436 //
437 // Return Res_hi + Res_lo in user-defined rounding control
438 //
439 // Step 5. Special Cases
440 //
441 // If pred is 0 where pred is obtained in
442 // frcap E, pred, V, U
443 //
444 // we are in one of those special cases of 0,+-inf or NaN
445 //
446 // If one of U and V is NaN, return U+V (which will generate
447 // invalid in case one is a signaling NaN). Otherwise,
448 // return the Result as described in the table
449 //
450 //
451 //
452 // \ Y |
453 // X \ | +0 | -0 | +inf | -inf | finite non-zero
454 // \ | | | | |
455 // ______________________________________________________
456 // | | | |
457 // +-0 | Invalid/ | pi/2 | -pi/2 | sign(Y)*pi/2
458 // | qNaN | | |
459 // --------------------------------------------------------
460 // | | | | |
461 // +inf | +0 | -0 | pi/4 | -pi/4 | sign(Y)*0
462 // --------------------------------------------------------
463 // | | | | |
464 // -inf | +pi | -pi | 3pi/4 | -3pi/4 | sign(Y)*pi
465 // --------------------------------------------------------
466 // finite | X>0? | pi/2 | -pi/2 |
467 // non-zero| sign(Y)*0: | | | N/A
468 // | sign(Y)*pi | | |
469 //
470 //
472 #include "libm_support.h"
474 ArgY_orig = f8
475 Result = f8
476 FR_RESULT = f8
477 ArgX_orig = f9
478 ArgX = f10
479 FR_X = f10
480 ArgY = f11
481 FR_Y = f11
482 s_Y = f12
483 U = f13
484 V = f14
485 E = f15
486 Q = f32
487 z_hi = f33
488 U_prime_hi = f34
489 U_prime_lo = f35
490 V_prime = f36
491 C_hi = f37
492 w_hi = f38
493 w_lo = f39
494 wsq = f40
495 poly = f41
496 Tbl_hi = f42
497 Tbl_lo = f43
498 P_hi = f44
499 P_lo = f45
500 A_hi = f46
501 A_lo = f47
502 sigma = f48
503 Res_hi = f49
504 Res_lo = f50
505 Z = f52
506 zsq = f53
507 z8 = f54
508 poly1 = f55
509 poly2 = f56
510 z_lo = f57
511 tmp = f58
512 P_1 = f59
513 Q_1 = f60
514 P_2 = f61
515 Q_2 = f62
516 P_3 = f63
517 Q_3 = f64
518 P_4 = f65
519 Q_4 = f66
520 P_5 = f67
521 P_6 = f68
522 P_7 = f69
523 P_8 = f70
524 TWO_TO_NEG3 = f71
525 U_hold = f72
526 C_hi_hold = f73
527 E_hold = f74
528 M = f75
529 ArgX_abs = f76
530 ArgY_abs = f77
531 Result_lo = f78
532 A_temp = f79
533 GR_SAVE_PFS = r33
534 GR_SAVE_B0 = r34
535 GR_SAVE_GP = r35
536 sign_X = r36
537 sign_Y = r37
538 swap = r38
539 table_ptr1 = r39
540 table_ptr2 = r40
541 k = r41
542 lookup = r42
543 exp_ArgX = r43
544 exp_ArgY = r44
545 exponent_Q = r45
546 significand_Q = r46
547 special = r47
548 special1 = r48
549 GR_Parameter_X = r49
550 GR_Parameter_Y = r50
551 GR_Parameter_RESULT = r51
552 GR_Parameter_TAG = r52
553 int_temp = r52
555 #ifdef _LIBC
556 .rodata
557 #else
558 .data
559 #endif
560 .align 64
562 Constants_atan:
563 ASM_TYPE_DIRECTIVE(Constants_atan,@object)
564 data4 0x54442D18, 0x3FF921FB, 0x248D3132, 0x3E000000
565 // double pi/2, single lo_pi/2, two**(-3)
566 data4 0xAAAAAAA3, 0xAAAAAAAA, 0x0000BFFD, 0x00000000 // P_1
567 data4 0xCCCC54B2, 0xCCCCCCCC, 0x00003FFC, 0x00000000 // P_2
568 data4 0x47E4D0C2, 0x92492492, 0x0000BFFC, 0x00000000 // P_3
569 data4 0x58870889, 0xE38E38E0, 0x00003FFB, 0x00000000 // P_4
570 data4 0x290149F8, 0xBA2E895B, 0x0000BFFB, 0x00000000 // P_5
571 data4 0x250F733D, 0x9D88E6D4, 0x00003FFB, 0x00000000 // P_6
572 data4 0xFB8745A0, 0x884E51FF, 0x0000BFFB, 0x00000000 // P_7
573 data4 0x394396BD, 0xE1C7412B, 0x00003FFA, 0x00000000 // P_8
574 data4 0xAAAAA52F, 0xAAAAAAAA, 0x0000BFFD, 0x00000000 // Q_1
575 data4 0xC75B60D3, 0xCCCCCCCC, 0x00003FFC, 0x00000000 // Q_2
576 data4 0x011F1940, 0x924923AD, 0x0000BFFC, 0x00000000 // Q_3
577 data4 0x2A5F89BD, 0xE36F716D, 0x00003FFB, 0x00000000 // Q_4
578 //
579 // Entries Tbl_hi (double precision)
580 // B = 1+Index/16+1/32 Index = 0
581 // Entries Tbl_lo (single precision)
582 // B = 1+Index/16+1/32 Index = 0
583 //
584 data4 0xA935BD8E, 0x3FE9A000, 0x23ACA08F, 0x00000000
585 //
586 // Entries Tbl_hi (double precision) Index = 0,1,...,15
587 // B = 2^(-1)*(1+Index/16+1/32)
588 // Entries Tbl_lo (single precision)
589 // Index = 0,1,...,15 B = 2^(-1)*(1+Index/16+1/32)
590 //
591 data4 0x7F175A34, 0x3FDE77EB, 0x238729EE, 0x00000000
592 data4 0x73C1A40B, 0x3FE0039C, 0x249334DB, 0x00000000
593 data4 0x5B5B43DA, 0x3FE0C614, 0x22CBA7D1, 0x00000000
594 data4 0x88BE7C13, 0x3FE1835A, 0x246310E7, 0x00000000
595 data4 0xE2CC9E6A, 0x3FE23B71, 0x236210E5, 0x00000000
596 data4 0x8406CBCA, 0x3FE2EE62, 0x2462EAF5, 0x00000000
597 data4 0x1CD41719, 0x3FE39C39, 0x24B73EF3, 0x00000000
598 data4 0x5B795B55, 0x3FE44506, 0x24C11260, 0x00000000
599 data4 0x5BB6EC04, 0x3FE4E8DE, 0x242519EE, 0x00000000
600 data4 0x1F732FBA, 0x3FE587D8, 0x24D4346C, 0x00000000
601 data4 0x115D7B8D, 0x3FE6220D, 0x24ED487B, 0x00000000
602 data4 0x920B3D98, 0x3FE6B798, 0x2495FF1E, 0x00000000
603 data4 0x8FBA8E0F, 0x3FE74897, 0x223D9531, 0x00000000
604 data4 0x289FA093, 0x3FE7D528, 0x242B0411, 0x00000000
605 data4 0x576CC2C5, 0x3FE85D69, 0x2335B374, 0x00000000
606 data4 0xA99CC05D, 0x3FE8E17A, 0x24C27CFB, 0x00000000
607 //
608 // Entries Tbl_hi (double precision) Index = 0,1,...,15
609 // B = 2^(-2)*(1+Index/16+1/32)
610 // Entries Tbl_lo (single precision)
611 // Index = 0,1,...,15 B = 2^(-2)*(1+Index/16+1/32)
612 //
613 data4 0x510665B5, 0x3FD025FA, 0x24263482, 0x00000000
614 data4 0x362431C9, 0x3FD1151A, 0x242C8DC9, 0x00000000
615 data4 0x67E47C95, 0x3FD20255, 0x245CF9BA, 0x00000000
616 data4 0x7A823CFE, 0x3FD2ED98, 0x235C892C, 0x00000000
617 data4 0x29271134, 0x3FD3D6D1, 0x2389BE52, 0x00000000
618 data4 0x586890E6, 0x3FD4BDEE, 0x24436471, 0x00000000
619 data4 0x175E0F4E, 0x3FD5A2E0, 0x2389DBD4, 0x00000000
620 data4 0x9F5FA6FD, 0x3FD68597, 0x2476D43F, 0x00000000
621 data4 0x52817501, 0x3FD76607, 0x24711774, 0x00000000
622 data4 0xB8DF95D7, 0x3FD84422, 0x23EBB501, 0x00000000
623 data4 0x7CD0C662, 0x3FD91FDE, 0x23883A0C, 0x00000000
624 data4 0x66168001, 0x3FD9F930, 0x240DF63F, 0x00000000
625 data4 0x5422058B, 0x3FDAD00F, 0x23FE261A, 0x00000000
626 data4 0x378624A5, 0x3FDBA473, 0x23A8CD0E, 0x00000000
627 data4 0x0AAD71F8, 0x3FDC7655, 0x2422D1D0, 0x00000000
628 data4 0xC9EC862B, 0x3FDD45AE, 0x2344A109, 0x00000000
629 //
630 // Entries Tbl_hi (double precision) Index = 0,1,...,15
631 // B = 2^(-3)*(1+Index/16+1/32)
632 // Entries Tbl_lo (single precision)
633 // Index = 0,1,...,15 B = 2^(-3)*(1+Index/16+1/32)
634 //
635 data4 0x84212B3D, 0x3FC068D5, 0x239874B6, 0x00000000
636 data4 0x41060850, 0x3FC16465, 0x2335E774, 0x00000000
637 data4 0x171A535C, 0x3FC25F6E, 0x233E36BE, 0x00000000
638 data4 0xEDEB99A3, 0x3FC359E8, 0x239680A3, 0x00000000
639 data4 0xC6092A9E, 0x3FC453CE, 0x230FB29E, 0x00000000
640 data4 0xBA11570A, 0x3FC54D18, 0x230C1418, 0x00000000
641 data4 0xFFB3AA73, 0x3FC645BF, 0x23F0564A, 0x00000000
642 data4 0xE8A7D201, 0x3FC73DBD, 0x23D4A5E1, 0x00000000
643 data4 0xE398EBC7, 0x3FC8350B, 0x23D4ADDA, 0x00000000
644 data4 0x7D050271, 0x3FC92BA3, 0x23BCB085, 0x00000000
645 data4 0x601081A5, 0x3FCA217E, 0x23BC841D, 0x00000000
646 data4 0x574D780B, 0x3FCB1696, 0x23CF4A8E, 0x00000000
647 data4 0x4D768466, 0x3FCC0AE5, 0x23BECC90, 0x00000000
648 data4 0x4E1D5395, 0x3FCCFE65, 0x2323DCD2, 0x00000000
649 data4 0x864C9D9D, 0x3FCDF110, 0x23F53F3A, 0x00000000
650 data4 0x451D980C, 0x3FCEE2E1, 0x23CCB11F, 0x00000000
652 data4 0x54442D18, 0x400921FB, 0x33145C07, 0x3CA1A626 // PI two doubles
653 data4 0x54442D18, 0x3FF921FB, 0x33145C07, 0x3C91A626 // PI_by_2 two dbles
654 data4 0x54442D18, 0x3FE921FB, 0x33145C07, 0x3C81A626 // PI_by_4 two dbles
655 data4 0x7F3321D2, 0x4002D97C, 0x4C9E8A0A, 0x3C9A7939 // 3PI_by_4 two dbles
656 ASM_SIZE_DIRECTIVE(Constants_atan)
659 .text
660 .proc atanl#
661 .global atanl#
662 .align 64
664 atanl:
665 { .mfb
666 nop.m 999
667 (p0) mov ArgX_orig = f1
668 (p0) br.cond.sptk atan2l ;;
669 }
670 .endp atanl
671 ASM_SIZE_DIRECTIVE(atanl)
673 .text
674 .proc atan2l#
675 .global atan2l#
676 #ifdef _LIBC
677 .proc __atan2l#
678 .global __atan2l#
679 .proc __ieee754_atan2l#
680 .global __ieee754_atan2l#
681 #endif
682 .align 64
685 atan2l:
686 #ifdef _LIBC
687 __atan2l:
688 __ieee754_atan2l:
689 #endif
690 { .mfi
691 alloc r32 = ar.pfs, 0, 17 , 4, 0
692 (p0) mov ArgY = ArgY_orig
693 }
694 { .mfi
695 nop.m 999
696 (p0) mov ArgX = ArgX_orig
697 nop.i 999
698 };;
699 { .mfi
700 nop.m 999
701 (p0) fclass.m.unc p7,p0 = ArgY_orig, 0x103
702 nop.i 999
703 }
704 { .mfi
705 nop.m 999
706 //
707 //
708 // Save original input args and load table ptr.
709 //
710 (p0) fclass.m.unc p6,p0 = ArgX_orig, 0x103
711 nop.i 999
712 };;
713 { .mfi
714 (p0) addl table_ptr1 = @ltoff(Constants_atan#), gp
715 (p0) fclass.m.unc p0,p9 = ArgY_orig, 0x1FF
716 nop.i 999 ;;
717 }
718 { .mfi
719 ld8 table_ptr1 = [table_ptr1]
720 (p0) fclass.m.unc p0,p8 = ArgX_orig, 0x1FF
721 nop.i 999
722 }
723 { .mfi
724 nop.m 999
725 (p0) fclass.m.unc p13,p0 = ArgY_orig, 0x0C3
726 nop.i 999 ;;
727 }
728 { .mfi
729 (p0) fclass.m.unc p12,p0 = ArgX_orig, 0x0C3
730 nop.i 999
731 }
734 //
735 // Check for NatVals.
736 // Check for everything - if false, then must be pseudo-zero
737 // or pseudo-nan (IA unsupporteds).
738 //
739 { .mib
740 nop.m 999
741 nop.i 999
742 (p6) br.cond.spnt L(ATANL_NATVAL) ;;
743 }
745 { .mib
746 nop.m 999
747 nop.i 999
748 (p7) br.cond.spnt L(ATANL_NATVAL) ;;
749 }
750 { .mib
751 (p0) ldfd P_hi = [table_ptr1],8
752 nop.i 999
753 (p8) br.cond.spnt L(ATANL_UNSUPPORTED) ;;
754 }
755 { .mbb
756 (p0) add table_ptr2 = 96, table_ptr1
757 (p9) br.cond.spnt L(ATANL_UNSUPPORTED)
758 //
759 // Load double precision high-order part of pi
760 //
761 (p12) br.cond.spnt L(ATANL_NAN) ;;
762 }
763 { .mfb
764 nop.m 999
765 (p0) fnorm.s1 ArgX = ArgX
766 (p13) br.cond.spnt L(ATANL_NAN) ;;
767 }
768 //
769 // Normalize the input argument.
770 // Branch out if NaN inputs
771 //
772 { .mmf
773 (p0) ldfs P_lo = [table_ptr1], 4
774 nop.m 999
775 (p0) fnorm.s1 ArgY = ArgY ;;
776 }
777 { .mmf
778 nop.m 999
779 (p0) ldfs TWO_TO_NEG3 = [table_ptr1], 180
780 //
781 // U = max(ArgX_abs,ArgY_abs)
782 // V = min(ArgX_abs,ArgY_abs)
783 // if PR1, swap = 0
784 // if PR2, swap = 1
785 //
786 (p0) mov M = f1 ;;
787 }
788 { .mfi
789 nop.m 999
790 //
791 // Get exp and sign of ArgX
792 // Get exp and sign of ArgY
793 // Load 2**(-3) and increment ptr to Q_4.
794 //
795 (p0) fmerge.s ArgX_abs = f1, ArgX
796 nop.i 999 ;;
797 }
798 //
799 // load single precision low-order part of pi = P_lo
800 //
801 { .mfi
802 (p0) getf.exp sign_X = ArgX
803 (p0) fmerge.s ArgY_abs = f1, ArgY
804 nop.i 999 ;;
805 }
806 { .mii
807 (p0) getf.exp sign_Y = ArgY
808 nop.i 999 ;;
809 (p0) shr sign_X = sign_X, 17 ;;
810 }
811 { .mii
812 nop.m 999
813 (p0) shr sign_Y = sign_Y, 17 ;;
814 (p0) cmp.eq.unc p8, p9 = 0x00000, sign_Y ;;
815 }
816 { .mfi
817 nop.m 999
818 //
819 // Is ArgX_abs >= ArgY_abs
820 // Is sign_Y == 0?
821 //
822 (p0) fmax.s1 U = ArgX_abs, ArgY_abs
823 nop.i 999
824 }
825 { .mfi
826 nop.m 999
827 //
828 // ArgX_abs = |ArgX|
829 // ArgY_abs = |ArgY|
830 // sign_X is sign bit of ArgX
831 // sign_Y is sign bit of ArgY
832 //
833 (p0) fcmp.ge.s1 p6, p7 = ArgX_abs, ArgY_abs
834 nop.i 999 ;;
835 }
836 { .mfi
837 nop.m 999
838 (p0) fmin.s1 V = ArgX_abs, ArgY_abs
839 nop.i 999 ;;
840 }
841 { .mfi
842 nop.m 999
843 (p8) fadd.s1 s_Y = f0, f1
844 (p6) cmp.eq.unc p10, p11 = 0x00000, sign_X
845 }
846 { .mii
847 (p6) add swap = r0, r0
848 nop.i 999 ;;
849 (p7) add swap = 1, r0
850 }
851 { .mfi
852 nop.m 999
853 //
854 // Let M = 1.0
855 // if p8, s_Y = 1.0
856 // if p9, s_Y = -1.0
857 //
858 (p10) fsub.s1 M = M, f1
859 nop.i 999 ;;
860 }
861 { .mfi
862 nop.m 999
863 (p9) fsub.s1 s_Y = f0, f1
864 nop.i 999 ;;
865 }
866 { .mfi
867 nop.m 999
868 (p0) frcpa.s1 E, p6 = V, U
869 nop.i 999 ;;
870 }
871 { .mbb
872 nop.m 999
873 //
874 // E = frcpa(V,U)
875 //
876 (p6) br.cond.sptk L(ATANL_STEP2)
877 (p0) br.cond.spnt L(ATANL_SPECIAL_HANDLING) ;;
878 }
879 L(ATANL_STEP2):
880 { .mfi
881 nop.m 999
882 (p0) fmpy.s1 Q = E, V
883 nop.i 999
884 }
885 { .mfi
886 nop.m 999
887 (p0) fcmp.eq.s0 p0, p9 = f1, ArgY_orig
888 nop.i 999 ;;
889 }
890 { .mfi
891 nop.m 999
892 //
893 // Is Q < 2**(-3)?
894 //
895 (p0) fcmp.eq.s0 p0, p8 = f1, ArgX_orig
896 nop.i 999
897 }
898 { .mfi
899 nop.m 999
900 (p11) fadd.s1 M = M, f1
901 nop.i 999 ;;
902 }
903 { .mlx
904 nop.m 999
905 // *************************************************
906 // ********************* STEP2 *********************
907 // *************************************************
908 (p0) movl special = 0x8400000000000000
909 }
910 { .mlx
911 nop.m 999
912 //
913 // lookup = b_1 b_2 b_3 B_4
914 //
915 (p0) movl special1 = 0x0000000000000100 ;;
916 }
917 { .mfi
918 nop.m 999
919 //
920 // Do fnorms to raise any denormal operand
921 // exceptions.
922 //
923 (p0) fmpy.s1 P_hi = M, P_hi
924 nop.i 999
925 }
926 { .mfi
927 nop.m 999
928 (p0) fmpy.s1 P_lo = M, P_lo
929 nop.i 999 ;;
930 }
931 { .mfi
932 nop.m 999
933 //
934 // Q = E * V
935 //
936 (p0) fcmp.lt.unc.s1 p6, p7 = Q, TWO_TO_NEG3
937 nop.i 999 ;;
938 }
939 { .mmb
940 (p0) getf.sig significand_Q = Q
941 (p0) getf.exp exponent_Q = Q
942 nop.b 999 ;;
943 }
944 { .mmi
945 nop.m 999 ;;
946 (p0) andcm k = 0x0003, exponent_Q
947 (p0) extr.u lookup = significand_Q, 59, 4 ;;
948 }
949 { .mib
950 nop.m 999
951 (p0) dep special = lookup, special, 59, 4
952 //
953 // Generate 1.b_1 b_2 b_3 b_4 1 0 0 0 ... 0
954 //
955 (p6) br.cond.spnt L(ATANL_POLY) ;;
956 }
957 { .mfi
958 (p0) cmp.eq.unc p8, p9 = 0x0000, k
959 (p0) fmpy.s1 P_hi = s_Y, P_hi
960 //
961 // We waited a few extra cycles so P_lo and P_hi could be calculated.
962 // Load the constant 256 for loading up table entries.
963 //
964 // *************************************************
965 // ******************** STEP3 **********************
966 // *************************************************
967 (p0) add table_ptr2 = 16, table_ptr1
968 }
969 //
970 // Let z_hi have exponent and sign of original Q
971 // Load the Tbl_hi(0) else, increment pointer.
972 //
973 { .mii
974 (p0) ldfe Q_4 = [table_ptr1], -16
975 (p0) xor swap = sign_X, swap ;;
976 (p9) sub k = k, r0, 1
977 }
978 { .mmi
979 (p0) setf.sig z_hi = special
980 (p0) ldfe Q_3 = [table_ptr1], -16
981 (p9) add table_ptr2 = 16, table_ptr2 ;;
982 }
983 //
984 // U_hold = U - U_prime_hi
985 // k = k * 256 - Result can be 0, 256, or 512.
986 //
987 { .mmb
988 (p0) ldfe Q_2 = [table_ptr1], -16
989 (p8) ldfd Tbl_hi = [table_ptr2], 8
990 nop.b 999 ;;
991 }
992 //
993 // U_prime_lo = U_hold + V * z_hi
994 // lookup -> lookup * 16 + k
995 //
996 { .mmi
997 (p0) ldfe Q_1 = [table_ptr1], -16 ;;
998 (p8) ldfs Tbl_lo = [table_ptr2], 8
999 //
1000 // U_prime_hi = U + V * z_hi
1001 // Load the Tbl_lo(0)
1002 //
1003 (p9) pmpy2.r k = k, special1 ;;
1004 }
1005 { .mii
1006 nop.m 999
1007 nop.i 999
1008 nop.i 999 ;;
1009 }
1010 { .mii
1011 nop.m 999
1012 nop.i 999
1013 nop.i 999 ;;
1014 }
1015 { .mii
1016 nop.m 999
1017 nop.i 999
1018 nop.i 999 ;;
1019 }
1020 { .mii
1021 nop.m 999
1022 nop.i 999 ;;
1023 (p9) shladd lookup = lookup, 0x0004, k ;;
1024 }
1025 { .mmi
1026 (p9) add table_ptr2 = table_ptr2, lookup ;;
1027 //
1028 // V_prime = V - U * z_hi
1029 //
1030 (p9) ldfd Tbl_hi = [table_ptr2], 8
1031 nop.i 999 ;;
1032 }
1033 { .mmf
1034 nop.m 999
1035 //
1036 // C_hi = frcpa(1,U_prime_hi)
1037 //
1038 (p9) ldfs Tbl_lo = [table_ptr2], 8
1039 //
1040 // z_hi = s exp 1.b_1 b_2 b_3 b_4 1 0 0 0 ... 0
1041 // Point to beginning of Tbl_hi entries - k = 0.
1042 //
1043 (p0) fmerge.se z_hi = Q, z_hi ;;
1044 }
1045 { .mfi
1046 nop.m 999
1047 (p0) fma.s1 U_prime_hi = V, z_hi, U
1048 nop.i 999
1049 }
1050 { .mfi
1051 nop.m 999
1052 (p0) fnma.s1 V_prime = U, z_hi, V
1053 nop.i 999 ;;
1054 }
1055 { .mfi
1056 nop.m 999
1057 (p0) mov A_hi = Tbl_hi
1058 nop.i 999 ;;
1059 }
1060 { .mfi
1061 nop.m 999
1062 (p0) fsub.s1 U_hold = U, U_prime_hi
1063 nop.i 999 ;;
1064 }
1065 { .mfi
1066 nop.m 999
1067 (p0) frcpa.s1 C_hi, p6 = f1, U_prime_hi
1068 nop.i 999 ;;
1069 }
1070 { .mfi
1071 (p0) cmp.eq.unc p7, p6 = 0x00000, swap
1072 (p0) fmpy.s1 A_hi = s_Y, A_hi
1073 nop.i 999 ;;
1074 }
1075 { .mfi
1076 nop.m 999
1077 //
1078 // poly = wsq * poly
1079 //
1080 (p7) fadd.s1 sigma = f0, f1
1081 nop.i 999 ;;
1082 }
1083 { .mfi
1084 nop.m 999
1085 (p0) fma.s1 U_prime_lo = z_hi, V, U_hold
1086 nop.i 999
1087 }
1088 { .mfi
1089 nop.m 999
1090 (p6) fsub.s1 sigma = f0, f1
1091 nop.i 999 ;;
1092 }
1093 { .mfi
1094 nop.m 999
1095 (p0) fnma.s1 C_hi_hold = C_hi, U_prime_hi, f1
1096 nop.i 999 ;;
1097 }
1098 { .mfi
1099 nop.m 999
1100 //
1101 // A_lo = A_lo + w_hi
1102 // A_hi = s_Y * A_hi
1103 //
1104 (p0) fma.s1 Res_hi = sigma, A_hi, P_hi
1105 nop.i 999 ;;
1106 }
1107 { .mfi
1108 nop.m 999
1109 //
1110 // C_hi_hold = 1 - C_hi * U_prime_hi (1)
1111 //
1112 (p0) fma.s1 C_hi = C_hi_hold, C_hi, C_hi
1113 nop.i 999 ;;
1114 }
1115 { .mfi
1116 nop.m 999
1117 //
1118 // C_hi = C_hi + C_hi * C_hi_hold (1)
1119 //
1120 (p0) fnma.s1 C_hi_hold = C_hi, U_prime_hi, f1
1121 nop.i 999 ;;
1122 }
1123 { .mfi
1124 nop.m 999
1125 //
1126 // C_hi_hold = 1 - C_hi * U_prime_hi (2)
1127 //
1128 (p0) fma.s1 C_hi = C_hi_hold, C_hi, C_hi
1129 nop.i 999 ;;
1130 }
1131 { .mfi
1132 nop.m 999
1133 //
1134 // C_hi = C_hi + C_hi * C_hi_hold (2)
1135 //
1136 (p0) fnma.s1 C_hi_hold = C_hi, U_prime_hi, f1
1137 nop.i 999 ;;
1138 }
1139 { .mfi
1140 nop.m 999
1141 //
1142 // C_hi_hold = 1 - C_hi * U_prime_hi (3)
1143 //
1144 (p0) fma.s1 C_hi = C_hi_hold, C_hi, C_hi
1145 nop.i 999 ;;
1146 }
1147 { .mfi
1148 nop.m 999
1149 //
1150 // C_hi = C_hi + C_hi * C_hi_hold (3)
1151 //
1152 (p0) fmpy.s1 w_hi = V_prime, C_hi
1153 nop.i 999 ;;
1154 }
1155 { .mfi
1156 nop.m 999
1157 //
1158 // w_hi = V_prime * C_hi
1159 //
1160 (p0) fmpy.s1 wsq = w_hi, w_hi
1161 nop.i 999
1162 }
1163 { .mfi
1164 nop.m 999
1165 (p0) fnma.s1 w_lo = w_hi, U_prime_hi, V_prime
1166 nop.i 999 ;;
1167 }
1168 { .mfi
1169 nop.m 999
1170 //
1171 // wsq = w_hi * w_hi
1172 // w_lo = = V_prime - w_hi * U_prime_hi
1173 //
1174 (p0) fma.s1 poly = wsq, Q_4, Q_3
1175 nop.i 999
1176 }
1177 { .mfi
1178 nop.m 999
1179 (p0) fnma.s1 w_lo = w_hi, U_prime_lo, w_lo
1180 nop.i 999 ;;
1181 }
1182 { .mfi
1183 nop.m 999
1184 //
1185 // poly = Q_3 + wsq * Q_4
1186 // w_lo = = w_lo - w_hi * U_prime_lo
1187 //
1188 (p0) fma.s1 poly = wsq, poly, Q_2
1189 nop.i 999
1190 }
1191 { .mfi
1192 nop.m 999
1193 (p0) fmpy.s1 w_lo = C_hi, w_lo
1194 nop.i 999 ;;
1195 }
1196 { .mfi
1197 nop.m 999
1198 //
1199 // poly = Q_2 + wsq * poly
1200 // w_lo = = w_lo * C_hi
1201 //
1202 (p0) fma.s1 poly = wsq, poly, Q_1
1203 nop.i 999
1204 }
1205 { .mfi
1206 nop.m 999
1207 (p0) fadd.s1 A_lo = Tbl_lo, w_lo
1208 nop.i 999 ;;
1209 }
1210 { .mfi
1211 nop.m 999
1212 //
1213 // Result = Res_hi + Res_lo * s_Y (User Supplied Rounding Mode)
1214 //
1215 (p0) fmpy.s0 Q_1 = Q_1, Q_1
1216 nop.i 999 ;;
1217 }
1218 { .mfi
1219 nop.m 999
1220 //
1221 // poly = Q_1 + wsq * poly
1222 // A_lo = Tbl_lo + w_lo
1223 // swap = xor(swap,sign_X)
1224 //
1225 (p0) fmpy.s1 poly = wsq, poly
1226 nop.i 999 ;;
1227 }
1228 { .mfi
1229 nop.m 999
1230 //
1231 // Is (swap) != 0 ?
1232 // poly = wsq * poly
1233 // A_hi = Tbl_hi
1234 //
1235 (p0) fmpy.s1 poly = w_hi, poly
1236 nop.i 999 ;;
1237 }
1238 { .mfi
1239 nop.m 999
1240 //
1241 // if (PR_1) sigma = -1.0
1242 // if (PR_2) sigma = 1.0
1243 //
1244 (p0) fadd.s1 A_lo = A_lo, poly
1245 nop.i 999 ;;
1246 }
1247 { .mfi
1248 nop.m 999
1249 //
1250 // P_hi = s_Y * P_hi
1251 // A_lo = A_lo + poly
1252 //
1253 (p0) fadd.s1 A_lo = A_lo, w_hi
1254 nop.i 999 ;;
1255 }
1256 { .mfi
1257 nop.m 999
1258 (p0) fma.s1 Res_lo = sigma, A_lo, P_lo
1259 nop.i 999 ;;
1260 }
1261 { .mfb
1262 nop.m 999
1263 //
1264 // Res_hi = P_hi + sigma * A_hi
1265 // Res_lo = P_lo + sigma * A_lo
1266 //
1267 (p0) fma.s0 Result = Res_lo, s_Y, Res_hi
1268 //
1269 // Raise inexact.
1270 //
1271 br.ret.sptk b0 ;;
1272 }
1273 //
1274 // poly1 = P_5 + zsq * poly1
1275 // poly2 = zsq * poly2
1276 //
1277 L(ATANL_POLY):
1278 { .mmf
1279 (p0) xor swap = sign_X, swap
1280 nop.m 999
1281 (p0) fnma.s1 E_hold = E, U, f1 ;;
1282 }
1283 { .mfi
1284 nop.m 999
1285 (p0) mov A_temp = Q
1286 //
1287 // poly1 = P_4 + zsq * poly1
1288 // swap = xor(swap,sign_X)
1289 //
1290 // sign_X gr_002
1291 // swap gr_004
1292 // poly1 = poly1 <== Done with poly1
1293 // poly1 = P_4 + zsq * poly1
1294 // swap = xor(swap,sign_X)
1295 //
1296 (p0) cmp.eq.unc p7, p6 = 0x00000, swap
1297 }
1298 { .mfi
1299 nop.m 999
1300 (p0) fmpy.s1 P_hi = s_Y, P_hi
1301 nop.i 999 ;;
1302 }
1303 { .mfi
1304 nop.m 999
1305 (p6) fsub.s1 sigma = f0, f1
1306 nop.i 999
1307 }
1308 { .mfi
1309 nop.m 999
1310 (p7) fadd.s1 sigma = f0, f1
1311 nop.i 999 ;;
1312 }
1314 // ***********************************************
1315 // ******************** STEP4 ********************
1316 // ***********************************************
1318 { .mmi
1319 nop.m 999
1320 (p0) addl table_ptr1 = @ltoff(Constants_atan#), gp
1321 nop.i 999
1322 }
1323 ;;
1325 { .mmi
1326 ld8 table_ptr1 = [table_ptr1]
1327 nop.m 999
1328 nop.i 999
1329 }
1330 ;;
1333 { .mfi
1334 nop.m 999
1335 (p0) fma.s1 E = E, E_hold, E
1336 //
1337 // Following:
1338 // Iterate 3 times E = E + E*(1.0 - E*U)
1339 // Also load P_8, P_7, P_6, P_5, P_4
1340 // E_hold = 1.0 - E * U (1)
1341 // A_temp = Q
1342 //
1343 (p0) add table_ptr1 = 128, table_ptr1 ;;
1344 }
1345 { .mmf
1346 nop.m 999
1347 //
1348 // E = E + E_hold*E (1)
1349 // Point to P_8.
1350 //
1351 (p0) ldfe P_8 = [table_ptr1], -16
1352 //
1353 // poly = z8*poly1 + poly2 (Typo in writeup)
1354 // Is (swap) != 0 ?
1355 //
1356 (p0) fnma.s1 z_lo = A_temp, U, V ;;
1357 }
1358 { .mmb
1359 nop.m 999
1360 //
1361 // E_hold = 1.0 - E * U (2)
1362 //
1363 (p0) ldfe P_7 = [table_ptr1], -16
1364 nop.b 999 ;;
1365 }
1366 { .mmb
1367 nop.m 999
1368 //
1369 // E = E + E_hold*E (2)
1370 //
1371 (p0) ldfe P_6 = [table_ptr1], -16
1372 nop.b 999 ;;
1373 }
1374 { .mmb
1375 nop.m 999
1376 //
1377 // E_hold = 1.0 - E * U (3)
1378 //
1379 (p0) ldfe P_5 = [table_ptr1], -16
1380 nop.b 999 ;;
1381 }
1382 { .mmf
1383 nop.m 999
1384 //
1385 // E = E + E_hold*E (3)
1386 //
1387 //
1388 // At this point E approximates 1/U to roughly working precision
1389 // z = V*E approximates V/U
1390 //
1391 (p0) ldfe P_4 = [table_ptr1], -16
1392 (p0) fnma.s1 E_hold = E, U, f1 ;;
1393 }
1394 { .mmb
1395 nop.m 999
1396 //
1397 // Z = V * E
1398 //
1399 (p0) ldfe P_3 = [table_ptr1], -16
1400 nop.b 999 ;;
1401 }
1402 { .mmb
1403 nop.m 999
1404 //
1405 // zsq = Z * Z
1406 //
1407 (p0) ldfe P_2 = [table_ptr1], -16
1408 nop.b 999 ;;
1409 }
1410 { .mmb
1411 nop.m 999
1412 //
1413 // z8 = zsq * zsq
1414 //
1415 (p0) ldfe P_1 = [table_ptr1], -16
1416 nop.b 999 ;;
1417 }
1418 { .mlx
1419 nop.m 999
1420 (p0) movl int_temp = 0x24005
1421 }
1422 { .mfi
1423 nop.m 999
1424 (p0) fma.s1 E = E, E_hold, E
1425 nop.i 999 ;;
1426 }
1427 { .mfi
1428 nop.m 999
1429 (p0) fnma.s1 E_hold = E, U, f1
1430 nop.i 999 ;;
1431 }
1432 { .mfi
1433 nop.m 999
1434 (p0) fma.s1 E = E, E_hold, E
1435 nop.i 999 ;;
1436 }
1437 { .mfi
1438 nop.m 999
1439 (p0) fmpy.s1 Z = V, E
1440 nop.i 999
1441 }
1442 { .mfi
1443 nop.m 999
1444 //
1445 // z_lo = V - A_temp * U
1446 // if (PR_2) sigma = 1.0
1447 //
1448 (p0) fmpy.s1 z_lo = z_lo, E
1449 nop.i 999 ;;
1450 }
1451 { .mfi
1452 nop.m 999
1453 (p0) fmpy.s1 zsq = Z, Z
1454 nop.i 999
1455 }
1456 { .mfi
1457 nop.m 999
1458 //
1459 // z_lo = z_lo * E
1460 // if (PR_1) sigma = -1.0
1461 //
1462 (p0) fadd.s1 A_hi = A_temp, z_lo
1463 nop.i 999 ;;
1464 }
1465 { .mfi
1466 nop.m 999
1467 //
1468 // z8 = z8 * z8
1469 //
1470 //
1471 // Now what we want to do is
1472 // poly1 = P_4 + zsq*(P_5 + zsq*(P_6 + zsq*(P_7 + zsq*P_8)))
1473 // poly2 = zsq*(P_1 + zsq*(P_2 + zsq*P_3))
1474 //
1475 (p0) fma.s1 poly1 = zsq, P_8, P_7
1476 nop.i 999
1477 }
1478 { .mfi
1479 nop.m 999
1480 (p0) fma.s1 poly2 = zsq, P_3, P_2
1481 nop.i 999 ;;
1482 }
1483 { .mfi
1484 nop.m 999
1485 (p0) fmpy.s1 z8 = zsq, zsq
1486 nop.i 999
1487 }
1488 { .mfi
1489 nop.m 999
1490 (p0) fsub.s1 A_temp = A_temp, A_hi
1491 nop.i 999 ;;
1492 }
1493 { .mfi
1494 nop.m 999
1495 //
1496 // A_lo = Z * poly + z_lo
1497 //
1498 (p0) fmerge.s tmp = A_hi, A_hi
1499 nop.i 999 ;;
1500 }
1501 { .mfi
1502 nop.m 999
1503 //
1504 // poly1 = P_7 + zsq * P_8
1505 // poly2 = P_2 + zsq * P_3
1506 //
1507 (p0) fma.s1 poly1 = zsq, poly1, P_6
1508 nop.i 999
1509 }
1510 { .mfi
1511 nop.m 999
1512 (p0) fma.s1 poly2 = zsq, poly2, P_1
1513 nop.i 999 ;;
1514 }
1515 { .mfi
1516 nop.m 999
1517 (p0) fmpy.s1 z8 = z8, z8
1518 nop.i 999
1519 }
1520 { .mfi
1521 nop.m 999
1522 (p0) fadd.s1 z_lo = A_temp, z_lo
1523 nop.i 999 ;;
1524 }
1525 { .mfi
1526 nop.m 999
1527 //
1528 // poly1 = P_6 + zsq * poly1
1529 // poly2 = P_2 + zsq * poly2
1530 //
1531 (p0) fma.s1 poly1 = zsq, poly1, P_5
1532 nop.i 999
1533 }
1534 { .mfi
1535 nop.m 999
1536 (p0) fmpy.s1 poly2 = poly2, zsq
1537 nop.i 999 ;;
1538 }
1539 { .mfi
1540 nop.m 999
1541 //
1542 // Result = Res_hi + Res_lo (User Supplied Rounding Mode)
1543 //
1544 (p0) fmpy.s1 P_5 = P_5, P_5
1545 nop.i 999 ;;
1546 }
1547 { .mfi
1548 nop.m 999
1549 (p0) fma.s1 poly1 = zsq, poly1, P_4
1550 nop.i 999 ;;
1551 }
1552 { .mfi
1553 nop.m 999
1554 (p0) fma.s1 poly = z8, poly1, poly2
1555 nop.i 999 ;;
1556 }
1557 { .mfi
1558 nop.m 999
1559 //
1560 // Fixup added to force inexact later -
1561 // A_hi = A_temp + z_lo
1562 // z_lo = (A_temp - A_hi) + z_lo
1563 //
1564 (p0) fma.s1 A_lo = Z, poly, z_lo
1565 nop.i 999 ;;
1566 }
1567 { .mfi
1568 nop.m 999
1569 (p0) fadd.s1 A_hi = tmp, A_lo
1570 nop.i 999 ;;
1571 }
1572 { .mfi
1573 nop.m 999
1574 (p0) fsub.s1 tmp = tmp, A_hi
1575 nop.i 999
1576 }
1577 { .mfi
1578 nop.m 999
1579 (p0) fmpy.s1 A_hi = s_Y, A_hi
1580 nop.i 999 ;;
1581 }
1582 { .mfi
1583 nop.m 999
1584 (p0) fadd.s1 A_lo = tmp, A_lo
1585 nop.i 999
1586 }
1587 { .mfi
1588 (p0) setf.exp tmp = int_temp
1589 //
1590 // P_hi = s_Y * P_hi
1591 // A_hi = s_Y * A_hi
1592 //
1593 (p0) fma.s1 Res_hi = sigma, A_hi, P_hi
1594 nop.i 999 ;;
1595 }
1596 { .mfi
1597 nop.m 999
1598 (p0) fclass.m.unc p6,p0 = A_lo, 0x007
1599 nop.i 999 ;;
1600 }
1601 { .mfi
1602 nop.m 999
1603 (p6) mov A_lo = tmp
1604 nop.i 999
1605 }
1606 { .mfi
1607 nop.m 999
1608 //
1609 // Res_hi = P_hi + sigma * A_hi
1610 //
1611 (p0) fsub.s1 tmp = P_hi, Res_hi
1612 nop.i 999 ;;
1613 }
1614 { .mfi
1615 nop.m 999
1616 //
1617 // tmp = P_hi - Res_hi
1618 //
1619 (p0) fma.s1 tmp = A_hi, sigma, tmp
1620 nop.i 999
1621 }
1622 { .mfi
1623 nop.m 999
1624 (p0) fma.s1 sigma = A_lo, sigma, P_lo
1625 nop.i 999 ;;
1626 }
1627 { .mfi
1628 nop.m 999
1629 //
1630 // tmp = sigma * A_hi + tmp
1631 // sigma = A_lo * sigma + P_lo
1632 //
1633 (p0) fma.s1 Res_lo = s_Y, sigma, tmp
1634 nop.i 999 ;;
1635 }
1636 { .mfb
1637 nop.m 999
1638 //
1639 // Res_lo = s_Y * sigma + tmp
1640 //
1641 (p0) fadd.s0 Result = Res_lo, Res_hi
1642 br.ret.sptk b0 ;;
1643 }
1644 L(ATANL_NATVAL):
1645 L(ATANL_UNSUPPORTED):
1646 L(ATANL_NAN):
1647 { .mfb
1648 nop.m 999
1649 (p0) fmpy.s0 Result = ArgX,ArgY
1650 (p0) br.ret.sptk b0 ;;
1651 }
1652 L(ATANL_SPECIAL_HANDLING):
1653 { .mfi
1654 nop.m 999
1655 (p0) fcmp.eq.s0 p0, p6 = f1, ArgY_orig
1656 nop.i 999
1657 }
1658 { .mfi
1659 nop.m 999
1660 (p0) fcmp.eq.s0 p0, p5 = f1, ArgX_orig
1661 nop.i 999 ;;
1662 }
1663 { .mfi
1664 nop.m 999
1665 (p0) fclass.m.unc p6, p7 = ArgY, 0x007
1666 nop.i 999
1667 }
1668 { .mlx
1669 nop.m 999
1670 (p0) movl special = 992
1671 }
1672 ;;
1675 { .mmi
1676 nop.m 999
1677 (p0) addl table_ptr1 = @ltoff(Constants_atan#), gp
1678 nop.i 999
1679 }
1680 ;;
1682 { .mmi
1683 ld8 table_ptr1 = [table_ptr1]
1684 nop.m 999
1685 nop.i 999
1686 }
1687 ;;
1690 { .mib
1691 (p0) add table_ptr1 = table_ptr1, special
1692 nop.i 999
1693 (p7) br.cond.spnt L(ATANL_ArgY_Not_ZERO) ;;
1694 }
1695 { .mmf
1696 (p0) ldfd Result = [table_ptr1], 8
1697 nop.m 999
1698 (p6) fclass.m.unc p14, p0 = ArgX, 0x035 ;;
1699 }
1700 { .mmf
1701 nop.m 999
1702 (p0) ldfd Result_lo = [table_ptr1], -8
1703 (p6) fclass.m.unc p15, p0 = ArgX, 0x036 ;;
1704 }
1705 { .mfi
1706 nop.m 999
1707 (p14) fmerge.s Result = ArgY, f0
1708 nop.i 999
1709 }
1710 { .mfi
1711 nop.m 999
1712 (p6) fclass.m.unc p13, p0 = ArgX, 0x007
1713 nop.i 999 ;;
1714 }
1715 { .mfi
1716 nop.m 999
1717 (p14) fmerge.s Result_lo = ArgY, f0
1718 nop.i 999 ;;
1719 }
1720 { .mfi
1721 (p13) mov GR_Parameter_TAG = 36
1722 nop.f 999
1723 nop.i 999 ;;
1724 }
1725 { .mfi
1726 nop.m 999
1727 //
1728 // Return sign_Y * 0 when ArgX > +0
1729 //
1730 (p15) fmerge.s Result = ArgY, Result
1731 nop.i 999 ;;
1732 }
1733 { .mfi
1734 nop.m 999
1735 (p15) fmerge.s Result_lo = ArgY, Result_lo
1736 nop.i 999 ;;
1737 }
1738 { .mfb
1739 nop.m 999
1740 //
1741 // Return sign_Y * 0 when ArgX < -0
1742 //
1743 (p0) fadd.s0 Result = Result, Result_lo
1744 (p13) br.cond.spnt __libm_error_region ;;
1745 }
1746 { .mib
1747 nop.m 999
1748 nop.i 999
1749 //
1750 // Call error support funciton for atan(0,0)
1751 //
1752 (p0) br.ret.sptk b0 ;;
1753 }
1754 L(ATANL_ArgY_Not_ZERO):
1755 { .mfi
1756 nop.m 999
1757 (p0) fclass.m.unc p9, p10 = ArgY, 0x023
1758 nop.i 999 ;;
1759 }
1760 { .mib
1761 nop.m 999
1762 nop.i 999
1763 (p10) br.cond.spnt L(ATANL_ArgY_Not_INF) ;;
1764 }
1765 { .mfi
1766 nop.m 999
1767 (p9) fclass.m.unc p6, p0 = ArgX, 0x017
1768 nop.i 999
1769 }
1770 { .mfi
1771 nop.m 999
1772 (p9) fclass.m.unc p7, p0 = ArgX, 0x021
1773 nop.i 999 ;;
1774 }
1775 { .mfi
1776 nop.m 999
1777 (p9) fclass.m.unc p8, p0 = ArgX, 0x022
1778 nop.i 999 ;;
1779 }
1780 { .mmi
1781 (p6) add table_ptr1 = 16, table_ptr1 ;;
1782 (p0) ldfd Result = [table_ptr1], 8
1783 nop.i 999 ;;
1784 }
1785 { .mfi
1786 (p0) ldfd Result_lo = [table_ptr1], -8
1787 nop.f 999
1788 nop.i 999 ;;
1789 }
1790 { .mfi
1791 nop.m 999
1792 (p6) fmerge.s Result = ArgY, Result
1793 nop.i 999 ;;
1794 }
1795 { .mfi
1796 nop.m 999
1797 (p6) fmerge.s Result_lo = ArgY, Result_lo
1798 nop.i 999 ;;
1799 }
1800 { .mfb
1801 nop.m 999
1802 (p6) fadd.s0 Result = Result, Result_lo
1803 (p6) br.ret.sptk b0 ;;
1804 }
1805 //
1806 // Load PI/2 and adjust its sign.
1807 // Return +PI/2 when ArgY = +Inf and ArgX = +/-0 or normal
1808 // Return -PI/2 when ArgY = -Inf and ArgX = +/-0 or normal
1809 //
1810 { .mmi
1811 (p7) add table_ptr1 = 32, table_ptr1 ;;
1812 (p7) ldfd Result = [table_ptr1], 8
1813 nop.i 999 ;;
1814 }
1815 { .mfi
1816 (p7) ldfd Result_lo = [table_ptr1], -8
1817 nop.f 999
1818 nop.i 999 ;;
1819 }
1820 { .mfi
1821 nop.m 999
1822 (p7) fmerge.s Result = ArgY, Result
1823 nop.i 999 ;;
1824 }
1825 { .mfi
1826 nop.m 999
1827 (p7) fmerge.s Result_lo = ArgY, Result_lo
1828 nop.i 999 ;;
1829 }
1830 { .mfb
1831 nop.m 999
1832 (p7) fadd.s0 Result = Result, Result_lo
1833 (p7) br.ret.sptk b0 ;;
1834 }
1835 //
1836 // Load PI/4 and adjust its sign.
1837 // Return +PI/4 when ArgY = +Inf and ArgX = +Inf
1838 // Return -PI/4 when ArgY = -Inf and ArgX = +Inf
1839 //
1840 { .mmi
1841 (p8) add table_ptr1 = 48, table_ptr1 ;;
1842 (p8) ldfd Result = [table_ptr1], 8
1843 nop.i 999 ;;
1844 }
1845 { .mfi
1846 (p8) ldfd Result_lo = [table_ptr1], -8
1847 nop.f 999
1848 nop.i 999 ;;
1849 }
1850 { .mfi
1851 nop.m 999
1852 (p8) fmerge.s Result = ArgY, Result
1853 nop.i 999 ;;
1854 }
1855 { .mfi
1856 nop.m 999
1857 (p8) fmerge.s Result_lo = ArgY, Result_lo
1858 nop.i 999 ;;
1859 }
1860 { .mfb
1861 nop.m 999
1862 (p8) fadd.s0 Result = Result, Result_lo
1863 (p8) br.ret.sptk b0 ;;
1864 }
1865 L(ATANL_ArgY_Not_INF):
1866 { .mfi
1867 nop.m 999
1868 //
1869 // Load PI/4 and adjust its sign.
1870 // Return +3PI/4 when ArgY = +Inf and ArgX = -Inf
1871 // Return -3PI/4 when ArgY = -Inf and ArgX = -Inf
1872 //
1873 (p0) fclass.m.unc p6, p0 = ArgX, 0x007
1874 nop.i 999
1875 }
1876 { .mfi
1877 nop.m 999
1878 (p0) fclass.m.unc p7, p0 = ArgX, 0x021
1879 nop.i 999 ;;
1880 }
1881 { .mfi
1882 nop.m 999
1883 (p0) fclass.m.unc p8, p0 = ArgX, 0x022
1884 nop.i 999 ;;
1885 }
1886 { .mmi
1887 (p6) add table_ptr1 = 16, table_ptr1 ;;
1888 (p6) ldfd Result = [table_ptr1], 8
1889 nop.i 999 ;;
1890 }
1891 { .mfi
1892 (p6) ldfd Result_lo = [table_ptr1], -8
1893 nop.f 999
1894 nop.i 999 ;;
1895 }
1896 { .mfi
1897 nop.m 999
1898 (p6) fmerge.s Result = ArgY, Result
1899 nop.i 999 ;;
1900 }
1901 { .mfi
1902 nop.m 999
1903 (p6) fmerge.s Result_lo = ArgY, Result_lo
1904 nop.i 999 ;;
1905 }
1906 { .mfb
1907 nop.m 999
1908 (p6) fadd.s0 Result = Result, Result_lo
1909 (p6) br.ret.spnt b0 ;;
1910 }
1911 { .mfi
1912 nop.m 999
1913 //
1914 // return = sign_Y * PI/2 when ArgX = 0
1915 //
1916 (p7) fmerge.s Result = ArgY, f0
1917 nop.i 999 ;;
1918 }
1919 { .mfb
1920 nop.m 999
1921 (p7) fnorm.s0 Result = Result
1922 (p7) br.ret.spnt b0 ;;
1923 }
1924 //
1925 // return = sign_Y * 0 when ArgX = Inf
1926 //
1927 { .mmi
1928 (p8) ldfd Result = [table_ptr1], 8 ;;
1929 (p8) ldfd Result_lo = [table_ptr1], -8
1930 nop.i 999 ;;
1931 }
1932 { .mfi
1933 nop.m 999
1934 (p8) fmerge.s Result = ArgY, Result
1935 nop.i 999 ;;
1936 }
1937 { .mfi
1938 nop.m 999
1939 (p8) fmerge.s Result_lo = ArgY, Result_lo
1940 nop.i 999 ;;
1941 }
1942 { .mfb
1943 nop.m 999
1944 (p8) fadd.s0 Result = Result, Result_lo
1945 (p8) br.ret.sptk b0 ;;
1946 }
1947 //
1948 // return = sign_Y * PI when ArgX = -Inf
1949 //
1950 .endp atan2l
1951 ASM_SIZE_DIRECTIVE(atan2l)
1952 ASM_SIZE_DIRECTIVE(__atan2l)
1953 ASM_SIZE_DIRECTIVE(__ieee754_atan2l)
1955 .proc __libm_error_region
1956 __libm_error_region:
1957 .prologue
1958 { .mfi
1959 add GR_Parameter_Y=-32,sp // Parameter 2 value
1960 nop.f 0
1961 .save ar.pfs,GR_SAVE_PFS
1962 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
1963 }
1964 { .mfi
1965 .fframe 64
1966 add sp=-64,sp // Create new stack
1967 nop.f 0
1968 mov GR_SAVE_GP=gp // Save gp
1969 };;
1970 { .mmi
1971 stfe [GR_Parameter_Y] = FR_Y,16 // Save Parameter 2 on stack
1972 add GR_Parameter_X = 16,sp // Parameter 1 address
1973 .save b0, GR_SAVE_B0
1974 mov GR_SAVE_B0=b0 // Save b0
1975 };;
1976 .body
1977 { .mib
1978 stfe [GR_Parameter_X] = FR_X // Store Parameter 1 on stack
1979 add GR_Parameter_RESULT = 0,GR_Parameter_Y
1980 nop.b 0 // Parameter 3 address
1981 }
1982 { .mib
1983 stfe [GR_Parameter_Y] = FR_RESULT // Store Parameter 3 on stack
1984 add GR_Parameter_Y = -16,GR_Parameter_Y
1985 br.call.sptk b0=__libm_error_support# // Call error handling function
1986 };;
1987 { .mmi
1988 nop.m 0
1989 nop.m 0
1990 add GR_Parameter_RESULT = 48,sp
1991 };;
1992 { .mmi
1993 ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack
1994 .restore sp
1995 add sp = 64,sp // Restore stack pointer
1996 mov b0 = GR_SAVE_B0 // Restore return address
1997 };;
1998 { .mib
1999 mov gp = GR_SAVE_GP // Restore gp
2000 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
2001 br.ret.sptk b0 // Return
2002 };;
2004 .endp __libm_error_region
2005 ASM_SIZE_DIRECTIVE(__libm_error_region)
2006 .type __libm_error_support#,@function
2007 .global __libm_error_support#