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