| blob | af24165d8e7d2a75dfaaebde3bf35f1dfe951c0a |
1 .file "asinf.s"
4 // Copyright (c) 2000 - 2003, Intel Corporation
5 // All rights reserved.
6 //
7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
11 // met:
12 //
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
15 //
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
19 //
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
22 // permission.
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://www.intel.com/software/products/opensource/libraries/num.htm.
40 // History
41 //==============================================================
42 // 02/02/00 Initial version
43 // 06/28/00 Improved speed
44 // 06/31/00 Changed register allocation because of some duplicate macros
45 // moved nan exit bundle up to gain a cycle.
46 // 08/08/00 Improved speed by avoiding SIR flush.
47 // 08/15/00 Bundle added after call to __libm_error_support to properly
48 // set [the previously overwritten] GR_Parameter_RESULT.
49 // 08/17/00 Changed predicate register macro-usage to direct predicate
50 // names due to an assembler bug.
51 // 10/17/00 Improved speed of x=0 and x=1 paths, set D flag if x denormal.
52 // 03/13/01 Corrected sign of imm1 value in dep instruction.
53 // 05/20/02 Cleaned up namespace and sf0 syntax
54 // 02/06/03 Reordered header: .section, .global, .proc, .align
57 // Description
58 //=========================================
59 // The asinf function computes the arc sine of x in the range [-pi,+pi].
60 // A doman error occurs for arguments not in the range [-1,+1].
61 // asinf(+-0) returns +-0
62 // asinf(x) returns a Nan and raises the invalid exception for |x| >1
64 // The acosf function returns the arc cosine in the range [0, +pi] radians.
65 // A doman error occurs for arguments not in the range [-1,+1].
66 // acosf(1) returns +0
67 // acosf(x) returns a Nan and raises the invalid exception for |x| >1
70 // |x| <= sqrt(2)/2. get Ax and Bx
72 // poly_p1 = x p1
73 // poly_p3 = x2 p4 + p3
74 // poly_p1 = x2 (poly_p1) + x = x2(x p1) + x
75 // poly_p2 = x2( poly_p3) + p2 = x2(x2 p4 + p3) + p2
77 // poly_Ax = x5(x2( poly_p3) + p2) + x2(x p1) + x
78 // = x5(x2(x2 p4 + p3) + p2) + x2(x p1) + x
80 // poly_p7 = x2 p8 + p7
81 // poly_p5 = x2 p6 + p5
83 // poly_p7 = x4 p9 + (poly_p7)
84 // poly_p7 = x4 p9 + (x2 p8 + p7)
85 // poly_Bx = x4 (x4 p9 + (x2 p8 + p7)) + x2 p6 + p5
87 // answer1 = x11(x4 (x4 p9 + (x2 p8 + p7)) + x2 p6 + p5) + x5(x2(x2 p4 + p3) + p2) + x2(x p1) + x
88 // = x19 p9 + x17 p8 + x15 p7 x13 p6 + x11 p5 + x9 p4 + x7 p3 + x5 p2 + x3 p1 + x
92 // |x| > sqrt(2)/2
94 // Get z = sqrt(1-x2)
96 // Get polynomial in t = 1-x2
98 // t2 = t t
99 // t4 = t2 t2
101 // poly_p4 = t p5 + p4
102 // poly_p1 = t p1 + 1
104 // poly_p6 = t p7 + p6
105 // poly_p2 = t p3 + p2
107 // poly_p8 = t p9 + p8
109 // poly_p4 = t2 poly_p6 + poly_p4
110 // = t2 (t p7 + p6) + (t p5 + p4)
112 // poly_p2 = t2 poly_p2 + poly_p1
113 // = t2 (t p3 + p2) + (t p1 + 1)
115 // poly_p4 = t4 poly_p8 + poly_p4
116 // = t4 (t p9 + p8) + (t2 (t p7 + p6) + (t p5 + p4))
118 // P(t) = poly_p2 + t4 poly_p8
119 // = t2 (t p3 + p2) + (t p1 + 1) + t4 (t4 (t p9 + p8) + (t2 (t p7 + p6) + (t p5 + p4)))
120 // = t3 p3 + t2 p2 + t p1 + 1 + t9 p9 + t8 p8 + t7 p7 + t6 p6 + t5 p5 + t4 p4
123 // answer2 = - sign(x) z P(t) + (sign(x) pi/2)
124 //
127 // Assembly macros
128 //=========================================
130 // predicate registers
131 //asinf_pred_LEsqrt2by2 = p7
132 //asinf_pred_GTsqrt2by2 = p8
134 // integer registers
135 ASINF_Addr1 = r33
136 ASINF_Addr2 = r34
137 ASINF_GR_1by2 = r35
139 ASINF_GR_3by2 = r36
140 ASINF_GR_5by2 = r37
142 GR_SAVE_B0 = r38
143 GR_SAVE_PFS = r39
144 GR_SAVE_GP = r40
146 GR_Parameter_X = r41
147 GR_Parameter_Y = r42
148 GR_Parameter_RESULT = r43
149 GR_Parameter_TAG = r44
151 // floating point registers
153 asinf_y = f32
154 asinf_abs_x = f33
155 asinf_x2 = f34
156 asinf_sgn_x = f35
158 asinf_1by2 = f36
159 asinf_3by2 = f37
160 asinf_5by2 = f38
161 asinf_coeff_P3 = f39
162 asinf_coeff_P8 = f40
164 asinf_coeff_P1 = f41
165 asinf_coeff_P4 = f42
166 asinf_coeff_P5 = f43
167 asinf_coeff_P2 = f44
168 asinf_coeff_P7 = f45
170 asinf_coeff_P6 = f46
171 asinf_coeff_P9 = f47
172 asinf_x2 = f48
173 asinf_x3 = f49
174 asinf_x4 = f50
176 asinf_x8 = f51
177 asinf_x5 = f52
178 asinf_const_piby2 = f53
179 asinf_const_sqrt2by2 = f54
180 asinf_x11 = f55
182 asinf_poly_p1 = f56
183 asinf_poly_p3 = f57
184 asinf_sinf1 = f58
185 asinf_poly_p2 = f59
186 asinf_poly_Ax = f60
188 asinf_poly_p7 = f61
189 asinf_poly_p5 = f62
190 asinf_sgnx_t4 = f63
191 asinf_poly_Bx = f64
192 asinf_t = f65
194 asinf_yby2 = f66
195 asinf_B = f67
196 asinf_B2 = f68
197 asinf_Az = f69
198 asinf_dz = f70
200 asinf_Sz = f71
201 asinf_d2z = f72
202 asinf_Fz = f73
203 asinf_z = f74
204 asinf_sgnx_z = f75
206 asinf_t2 = f76
207 asinf_2poly_p4 = f77
208 asinf_2poly_p6 = f78
209 asinf_2poly_p1 = f79
210 asinf_2poly_p2 = f80
212 asinf_2poly_p8 = f81
213 asinf_t4 = f82
214 asinf_Pt = f83
215 asinf_sgnx_2poly_p2 = f84
216 asinf_sgn_x_piby2 = f85
218 asinf_poly_p7a = f86
219 asinf_2poly_p4a = f87
220 asinf_2poly_p4b = f88
221 asinf_2poly_p2a = f89
222 asinf_poly_p1a = f90
228 // Data tables
229 //==============================================================
231 RODATA
233 .align 16
235 LOCAL_OBJECT_START(asinf_coeff_1_table)
236 data8 0x3FC5555607DCF816 // P1
237 data8 0x3F9CF81AD9BAB2C6 // P4
238 data8 0x3FC59E0975074DF3 // P7
239 data8 0xBFA6F4CC2780AA1D // P6
240 data8 0x3FC2DD45292E93CB // P9
241 data8 0x3fe6a09e667f3bcd // sqrt(2)/2
242 LOCAL_OBJECT_END(asinf_coeff_1_table)
244 LOCAL_OBJECT_START(asinf_coeff_2_table)
245 data8 0x3FA6F108E31EFBA6 // P3
246 data8 0xBFCA31BF175D82A0 // P8
247 data8 0x3FA30C0337F6418B // P5
248 data8 0x3FB332C9266CB1F9 // P2
249 data8 0x3ff921fb54442d18 // pi_by_2
250 LOCAL_OBJECT_END(asinf_coeff_2_table)
253 .section .text
254 GLOBAL_LIBM_ENTRY(asinf)
256 // Load the addresses of the two tables.
257 // Then, load the coefficients and other constants.
259 { .mfi
260 alloc r32 = ar.pfs,1,8,4,0
261 fnma.s1 asinf_t = f8,f8,f1
262 dep.z ASINF_GR_1by2 = 0x3f,24,8 // 0x3f000000
263 }
264 { .mfi
265 addl ASINF_Addr1 = @ltoff(asinf_coeff_1_table),gp
266 fma.s1 asinf_x2 = f8,f8,f0
267 addl ASINF_Addr2 = @ltoff(asinf_coeff_2_table),gp ;;
268 }
271 { .mfi
272 ld8 ASINF_Addr1 = [ASINF_Addr1]
273 fmerge.s asinf_abs_x = f1,f8
274 dep ASINF_GR_3by2 = -1,r0,22,8 // 0x3fc00000
275 }
276 { .mlx
277 nop.m 999
278 movl ASINF_GR_5by2 = 0x40200000;;
279 }
283 { .mfi
284 setf.s asinf_1by2 = ASINF_GR_1by2
285 fmerge.s asinf_sgn_x = f8,f1
286 nop.i 999
287 }
288 { .mfi
289 ld8 ASINF_Addr2 = [ASINF_Addr2]
290 nop.f 0
291 nop.i 999;;
292 }
295 { .mfi
296 setf.s asinf_5by2 = ASINF_GR_5by2
297 fcmp.lt.s1 p11,p12 = f8,f0
298 nop.i 999;;
299 }
301 { .mmf
302 ldfpd asinf_coeff_P1,asinf_coeff_P4 = [ASINF_Addr1],16
303 setf.s asinf_3by2 = ASINF_GR_3by2
304 fclass.m.unc p8,p0 = f8, 0xc3 ;; //@qnan | @snan
305 }
308 { .mfi
309 ldfpd asinf_coeff_P7,asinf_coeff_P6 = [ASINF_Addr1],16
310 fma.s1 asinf_t2 = asinf_t,asinf_t,f0
311 nop.i 999
312 }
313 { .mfi
314 ldfpd asinf_coeff_P3,asinf_coeff_P8 = [ASINF_Addr2],16
315 fma.s1 asinf_x4 = asinf_x2,asinf_x2,f0
316 nop.i 999;;
317 }
320 { .mfi
321 ldfpd asinf_coeff_P9,asinf_const_sqrt2by2 = [ASINF_Addr1]
322 fclass.m.unc p10,p0 = f8, 0x07 //@zero
323 nop.i 999
324 }
325 { .mfi
326 ldfpd asinf_coeff_P5,asinf_coeff_P2 = [ASINF_Addr2],16
327 fma.s1 asinf_x3 = f8,asinf_x2,f0
328 nop.i 999;;
329 }
332 { .mfi
333 ldfd asinf_const_piby2 = [ASINF_Addr2]
334 frsqrta.s1 asinf_B,p0 = asinf_t
335 nop.i 999
336 }
337 { .mfb
338 nop.m 999
339 (p8) fma.s.s0 f8 = f8,f1,f0
340 (p8) br.ret.spnt b0 ;; // Exit if x=nan
341 }
344 { .mfb
345 nop.m 999
346 fcmp.eq.s1 p6,p0 = asinf_abs_x,f1
347 (p10) br.ret.spnt b0 ;; // Exit if x=0
348 }
350 { .mfi
351 nop.m 999
352 fcmp.gt.s1 p9,p0 = asinf_abs_x,f1
353 nop.i 999;;
354 }
356 { .mfi
357 nop.m 999
358 fma.s1 asinf_x8 = asinf_x4,asinf_x4,f0
359 nop.i 999
360 }
361 { .mfb
362 nop.m 999
363 fma.s1 asinf_t4 = asinf_t2,asinf_t2,f0
364 (p6) br.cond.spnt ASINF_ABS_ONE ;; // Branch if |x|=1
365 }
367 { .mfi
368 nop.m 999
369 fma.s1 asinf_x5 = asinf_x2,asinf_x3,f0
370 nop.i 999
371 }
372 { .mfb
373 (p9) mov GR_Parameter_TAG = 62
374 fma.s1 asinf_yby2 = asinf_t,asinf_1by2,f0
375 (p9) br.cond.spnt __libm_error_region ;; // Branch if |x|>1
376 }
379 { .mfi
380 nop.m 999
381 fma.s1 asinf_Az = asinf_t,asinf_B,f0
382 nop.i 999
383 }
384 { .mfi
385 nop.m 999
386 fma.s1 asinf_B2 = asinf_B,asinf_B,f0
387 nop.i 999;;
388 }
390 { .mfi
391 nop.m 999
392 fma.s1 asinf_poly_p1 = f8,asinf_coeff_P1,f0
393 nop.i 999
394 }
395 { .mfi
396 nop.m 999
397 fma.s1 asinf_2poly_p1 = asinf_coeff_P1,asinf_t,f1
398 nop.i 999;;
399 }
401 { .mfi
402 nop.m 999
403 fma.s1 asinf_poly_p3 = asinf_coeff_P4,asinf_x2,asinf_coeff_P3
404 nop.i 999
405 }
406 { .mfi
407 nop.m 999
408 fma.s1 asinf_2poly_p6 = asinf_coeff_P7,asinf_t,asinf_coeff_P6
409 nop.i 999;;
410 }
412 { .mfi
413 nop.m 999
414 fma.s1 asinf_poly_p7 = asinf_x2,asinf_coeff_P8,asinf_coeff_P7
415 nop.i 999
416 }
417 { .mfi
418 nop.m 999
419 fma.s1 asinf_2poly_p2 = asinf_coeff_P3,asinf_t,asinf_coeff_P2
420 nop.i 999;;
421 }
424 { .mfi
425 nop.m 999
426 fma.s1 asinf_poly_p5 = asinf_x2,asinf_coeff_P6,asinf_coeff_P5
427 nop.i 999
428 }
429 { .mfi
430 nop.m 999
431 fma.s1 asinf_2poly_p4 = asinf_coeff_P5,asinf_t,asinf_coeff_P4
432 nop.i 999;;
433 }
436 { .mfi
437 nop.m 999
438 fma.d.s1 asinf_x11 = asinf_x8,asinf_x3,f0
439 nop.i 999
440 }
441 { .mfi
442 nop.m 999
443 fnma.s1 asinf_dz = asinf_B2,asinf_yby2,asinf_1by2
444 nop.i 999;;
445 }
448 { .mfi
449 nop.m 999
450 fma.s1 asinf_poly_p1a = asinf_x2,asinf_poly_p1,f8
451 nop.i 999
452 }
453 { .mfi
454 nop.m 999
455 fma.s1 asinf_2poly_p8 = asinf_coeff_P9,asinf_t,asinf_coeff_P8
456 nop.i 999;;
457 }
460 // Get the absolute value of x and determine the region in which x lies
462 { .mfi
463 nop.m 999
464 fcmp.le.s1 p7,p8 = asinf_abs_x,asinf_const_sqrt2by2
465 nop.i 999
466 }
467 { .mfi
468 nop.m 999
469 fma.s1 asinf_poly_p2 = asinf_x2,asinf_poly_p3,asinf_coeff_P2
470 nop.i 999;;
471 }
474 { .mfi
475 nop.m 999
476 fma.s1 asinf_poly_p7a = asinf_x4,asinf_coeff_P9,asinf_poly_p7
477 nop.i 999
478 }
479 { .mfi
480 nop.m 999
481 fma.s1 asinf_2poly_p2a = asinf_2poly_p2,asinf_t2,asinf_2poly_p1
482 nop.i 999;;
483 }
486 { .mfi
487 nop.m 999
488 (p8) fma.s1 asinf_sgnx_t4 = asinf_sgn_x,asinf_t4,f0
489 nop.i 999
490 }
491 { .mfi
492 nop.m 999
493 (p8) fma.s1 asinf_2poly_p4a = asinf_2poly_p6,asinf_t2,asinf_2poly_p4
494 nop.i 999;;
495 }
498 { .mfi
499 nop.m 999
500 (p8) fma.s1 asinf_Sz = asinf_5by2,asinf_dz,asinf_3by2
501 nop.i 999
502 }
503 { .mfi
504 nop.m 999
505 (p8) fma.s1 asinf_d2z = asinf_dz,asinf_dz,f0
506 nop.i 999;;
507 }
510 { .mfi
511 nop.m 999
512 (p8) fma.s1 asinf_sgn_x_piby2 = asinf_sgn_x,asinf_const_piby2,f0
513 nop.i 999
514 }
515 { .mfi
516 nop.m 999
517 (p7) fma.d.s1 asinf_poly_Ax = asinf_x5,asinf_poly_p2,asinf_poly_p1a
518 nop.i 999;;
519 }
521 { .mfi
522 nop.m 999
523 (p7) fma.d.s1 asinf_poly_Bx = asinf_x4,asinf_poly_p7a,asinf_poly_p5
524 nop.i 999
525 }
526 { .mfi
527 nop.m 999
528 (p8) fma.s1 asinf_sgnx_2poly_p2 = asinf_sgn_x,asinf_2poly_p2a,f0
529 nop.i 999;;
530 }
532 { .mfi
533 nop.m 999
534 fcmp.eq.s0 p6,p0 = f8,f0 // Only purpose is to set D if x denormal
535 nop.i 999
536 }
537 { .mfi
538 nop.m 999
539 (p8) fma.s1 asinf_2poly_p4b = asinf_2poly_p8,asinf_t4,asinf_2poly_p4a
540 nop.i 999;;
541 }
544 { .mfi
545 nop.m 999
546 (p8) fma.s1 asinf_Fz = asinf_d2z,asinf_Sz,asinf_dz
547 nop.i 999;;
548 }
551 { .mfi
552 nop.m 999
553 (p8) fma.d.s1 asinf_Pt = asinf_2poly_p4b,asinf_sgnx_t4,asinf_sgnx_2poly_p2
554 nop.i 999;;
555 }
557 { .mfi
558 nop.m 999
559 (p8) fma.d.s1 asinf_z = asinf_Az,asinf_Fz,asinf_Az
560 nop.i 999;;
561 }
563 .pred.rel "mutex",p8,p7 //asinf_pred_GTsqrt2by2,asinf_pred_LEsqrt2by2
564 { .mfi
565 nop.m 999
566 (p8) fnma.s.s0 f8 = asinf_z,asinf_Pt,asinf_sgn_x_piby2
567 nop.i 999
568 }
570 { .mfb
571 nop.m 999
572 (p7) fma.s.s0 f8 = asinf_x11,asinf_poly_Bx,asinf_poly_Ax
573 br.ret.sptk b0 ;;
574 }
576 ASINF_ABS_ONE:
577 // Here for short exit if |x|=1
578 { .mfb
579 nop.m 999
580 fma.s.s0 f8 = asinf_sgn_x,asinf_const_piby2,f0
581 br.ret.sptk b0
582 }
583 ;;
585 GLOBAL_LIBM_END(asinf)
587 // Stack operations when calling error support.
588 // (1) (2)
589 // sp -> + psp -> +
590 // | |
591 // | | <- GR_Y
592 // | |
593 // | <-GR_Y Y2->|
594 // | |
595 // | | <- GR_X
596 // | |
597 // sp-64 -> + sp -> +
598 // save ar.pfs save b0
599 // save gp
602 // Stack operations when calling error support.
603 // (3) (call) (4)
604 // psp -> + sp -> +
605 // | |
606 // R3 ->| <- GR_RESULT | -> f8
607 // | |
608 // Y2 ->| <- GR_Y |
609 // | |
610 // X1 ->| |
611 // | |
612 // sp -> + +
613 // restore gp
614 // restore ar.pfs
616 LOCAL_LIBM_ENTRY(__libm_error_region)
617 .prologue
618 { .mfi
619 add GR_Parameter_Y=-32,sp // Parameter 2 value
620 nop.f 999
621 .save ar.pfs,GR_SAVE_PFS
622 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
623 }
624 { .mfi
625 .fframe 64
626 add sp=-64,sp // Create new stack
627 nop.f 0
628 mov GR_SAVE_GP=gp // Save gp
629 };;
630 { .mmi
631 stfs [GR_Parameter_Y] = f1,16 // Store Parameter 2 on stack
632 add GR_Parameter_X = 16,sp // Parameter 1 address
633 .save b0, GR_SAVE_B0
634 mov GR_SAVE_B0=b0 // Save b0
635 };;
637 .body
638 { .mfi
639 nop.m 0
640 frcpa.s0 f9,p0 = f0,f0
641 nop.i 0
642 };;
644 { .mib
645 stfs [GR_Parameter_X] = f8 // Store Parameter 1 on stack
646 add GR_Parameter_RESULT = 0,GR_Parameter_Y
647 nop.b 0 // Parameter 3 address
648 }
649 { .mib
650 stfs [GR_Parameter_Y] = f9 // Store Parameter 3 on stack
651 add GR_Parameter_Y = -16,GR_Parameter_Y
652 br.call.sptk b0=__libm_error_support# // Call error handling function
653 };;
654 { .mmi
655 nop.m 0
656 nop.m 0
657 add GR_Parameter_RESULT = 48,sp
658 };;
660 { .mmi
661 ldfs f8 = [GR_Parameter_RESULT] // Get return result off stack
662 .restore sp
663 add sp = 64,sp // Restore stack pointer
664 mov b0 = GR_SAVE_B0 // Restore return address
665 };;
666 { .mib
667 mov gp = GR_SAVE_GP // Restore gp
668 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
669 br.ret.sptk b0 // Return
670 };;
672 LOCAL_LIBM_END(__libm_error_region)
674 .type __libm_error_support#,@function
675 .global __libm_error_support#