summaryrefslogtreecommitdiff
path: root/sysdeps/ia64/fpu/s_cos.S
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ia64/fpu/s_cos.S')
-rw-r--r--sysdeps/ia64/fpu/s_cos.S3497
1 files changed, 382 insertions, 3115 deletions
diff --git a/sysdeps/ia64/fpu/s_cos.S b/sysdeps/ia64/fpu/s_cos.S
index 6540aec724..fc121fce19 100644
--- a/sysdeps/ia64/fpu/s_cos.S
+++ b/sysdeps/ia64/fpu/s_cos.S
@@ -1,10 +1,10 @@
.file "sincos.s"
-// Copyright (C) 2000, 2001, Intel Corporation
+
+// Copyright (c) 2000 - 2005, Intel Corporation
// All rights reserved.
//
-// Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
-// and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
+// Contributed 2000 by the Intel Numerics Group, Intel Corporation
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
@@ -20,7 +20,7 @@
// * The name of Intel Corporation may not be used to endorse or promote
// products derived from this software without specific prior written
// permission.
-//
+
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
@@ -35,17 +35,25 @@
//
// Intel Corporation is the author of this code, and requests that all
// problem reports or change requests be submitted to it directly at
-// http://developer.intel.com/opensource.
+// http://www.intel.com/software/products/opensource/libraries/num.htm.
//
// History
//==============================================================
-// 2/02/00 Initial revision
-// 4/02/00 Unwind support added.
-// 6/16/00 Updated tables to enforce symmetry
-// 8/31/00 Saved 2 cycles in main path, and 9 in other paths.
-// 9/20/00 The updated tables regressed to an old version, so reinstated them
+// 02/02/00 Initial version
+// 04/02/00 Unwind support added.
+// 06/16/00 Updated tables to enforce symmetry
+// 08/31/00 Saved 2 cycles in main path, and 9 in other paths.
+// 09/20/00 The updated tables regressed to an old version, so reinstated them
// 10/18/00 Changed one table entry to ensure symmetry
-// 1/03/01 Improved speed, fixed flag settings for small arguments.
+// 01/03/01 Improved speed, fixed flag settings for small arguments.
+// 02/18/02 Large arguments processing routine excluded
+// 05/20/02 Cleaned up namespace and sf0 syntax
+// 06/03/02 Insure inexact flag set for large arg result
+// 09/05/02 Work range is widened by reduction strengthen (3 parts of Pi/16)
+// 02/10/03 Reordered header: .section, .global, .proc, .align
+// 08/08/03 Improved performance
+// 10/28/04 Saved sincos_r_sincos to avoid clobber by dynamic loader
+// 03/31/05 Reformatted delimiters between data tables
// API
//==============================================================
@@ -63,9 +71,13 @@
// nfloat = Round result to integer (round-to-nearest)
//
// r = x - nfloat * pi/2^k
-// Do this as (x - nfloat * HIGH(pi/2^k)) - nfloat * LOW(pi/2^k) for increased accuracy.
+// Do this as ((((x - nfloat * HIGH(pi/2^k))) -
+// nfloat * LOW(pi/2^k)) -
+// nfloat * LOWEST(pi/2^k) for increased accuracy.
// pi/2^k is stored as two numbers that when added make pi/2^k.
// pi/2^k = HIGH(pi/2^k) + LOW(pi/2^k)
+// HIGH and LOW parts are rounded to zero values,
+// and LOWEST is rounded to nearest one.
//
// x = (nfloat * pi/2^k) + r
// r is small enough that we can use a polynomial approximation
@@ -121,7 +133,7 @@
//
// as follows
//
-// Sm = Sin(Mpi/2^k) and Cm = Cos(Mpi/2^k)
+// S[m] = Sin(Mpi/2^k) and C[m] = Cos(Mpi/2^k)
// rsq = r*r
//
//
@@ -141,32 +153,31 @@
//
// P = r + rcub * P
//
-// Answer = Sm Cos(r) + Cm P
+// Answer = S[m] Cos(r) + [Cm] P
//
// Cos(r) = 1 + rsq Q
// Cos(r) = 1 + r^2 Q
// Cos(r) = 1 + r^2 (q1 + r^2q2 + r^4q3 + r^6q4)
// Cos(r) = 1 + r^2q1 + r^4q2 + r^6q3 + r^8q4 + ...
//
-// Sm Cos(r) = Sm(1 + rsq Q)
-// Sm Cos(r) = Sm + Sm rsq Q
-// Sm Cos(r) = Sm + s_rsq Q
-// Q = Sm + s_rsq Q
+// S[m] Cos(r) = S[m](1 + rsq Q)
+// S[m] Cos(r) = S[m] + Sm rsq Q
+// S[m] Cos(r) = S[m] + s_rsq Q
+// Q = S[m] + s_rsq Q
//
// Then,
//
-// Answer = Q + Cm P
+// Answer = Q + C[m] P
-#include "libm_support.h"
// Registers used
//==============================================================
// general input registers:
-// r14 -> r19
-// r32 -> r45
+// r14 -> r26
+// r32 -> r35
// predicate registers used:
-// p6 -> p14
+// p6 -> p11
// floating-point registers used
// f9 -> f15
@@ -174,99 +185,94 @@
// Assembly macros
//==============================================================
-sind_NORM_f8 = f9
-sind_W = f10
-sind_int_Nfloat = f11
-sind_Nfloat = f12
+sincos_NORM_f8 = f9
+sincos_W = f10
+sincos_int_Nfloat = f11
+sincos_Nfloat = f12
-sind_r = f13
-sind_rsq = f14
-sind_rcub = f15
+sincos_r = f13
+sincos_rsq = f14
+sincos_rcub = f15
+sincos_save_tmp = f15
-sind_Inv_Pi_by_16 = f32
-sind_Pi_by_16_hi = f33
-sind_Pi_by_16_lo = f34
+sincos_Inv_Pi_by_16 = f32
+sincos_Pi_by_16_1 = f33
+sincos_Pi_by_16_2 = f34
-sind_Inv_Pi_by_64 = f35
-sind_Pi_by_64_hi = f36
-sind_Pi_by_64_lo = f37
+sincos_Inv_Pi_by_64 = f35
-sind_Sm = f38
-sind_Cm = f39
+sincos_Pi_by_16_3 = f36
-sind_P1 = f40
-sind_Q1 = f41
-sind_P2 = f42
-sind_Q2 = f43
-sind_P3 = f44
-sind_Q3 = f45
-sind_P4 = f46
-sind_Q4 = f47
+sincos_r_exact = f37
-sind_P_temp1 = f48
-sind_P_temp2 = f49
+sincos_Sm = f38
+sincos_Cm = f39
-sind_Q_temp1 = f50
-sind_Q_temp2 = f51
+sincos_P1 = f40
+sincos_Q1 = f41
+sincos_P2 = f42
+sincos_Q2 = f43
+sincos_P3 = f44
+sincos_Q3 = f45
+sincos_P4 = f46
+sincos_Q4 = f47
-sind_P = f52
-sind_Q = f53
+sincos_P_temp1 = f48
+sincos_P_temp2 = f49
-sind_srsq = f54
+sincos_Q_temp1 = f50
+sincos_Q_temp2 = f51
-sind_SIG_INV_PI_BY_16_2TO61 = f55
-sind_RSHF_2TO61 = f56
-sind_RSHF = f57
-sind_2TOM61 = f58
-sind_NFLOAT = f59
-sind_W_2TO61_RSH = f60
-
-fp_tmp = f61
-
-/////////////////////////////////////////////////////////////
+sincos_P = f52
+sincos_Q = f53
-sind_AD_1 = r33
-sind_AD_2 = r34
-sind_exp_limit = r35
-sind_r_signexp = r36
-sind_AD_beta_table = r37
-sind_r_sincos = r38
+sincos_srsq = f54
-sind_r_exp = r39
-sind_r_17_ones = r40
+sincos_SIG_INV_PI_BY_16_2TO61 = f55
+sincos_RSHF_2TO61 = f56
+sincos_RSHF = f57
+sincos_2TOM61 = f58
+sincos_NFLOAT = f59
+sincos_W_2TO61_RSH = f60
-sind_GR_sig_inv_pi_by_16 = r14
-sind_GR_rshf_2to61 = r15
-sind_GR_rshf = r16
-sind_GR_exp_2tom61 = r17
-sind_GR_n = r18
-sind_GR_m = r19
-sind_GR_32m = r19
+fp_tmp = f61
-gr_tmp = r41
-GR_SAVE_PFS = r41
-GR_SAVE_B0 = r42
-GR_SAVE_GP = r43
-
-
-#ifdef _LIBC
-.rodata
-#else
-.data
-#endif
+/////////////////////////////////////////////////////////////
+sincos_GR_sig_inv_pi_by_16 = r14
+sincos_GR_rshf_2to61 = r15
+sincos_GR_rshf = r16
+sincos_GR_exp_2tom61 = r17
+sincos_GR_n = r18
+sincos_GR_m = r19
+sincos_GR_32m = r19
+sincos_GR_all_ones = r19
+sincos_AD_1 = r20
+sincos_AD_2 = r21
+sincos_exp_limit = r22
+sincos_r_signexp = r23
+sincos_r_17_ones = r24
+sincos_r_sincos = r25
+sincos_r_exp = r26
+
+GR_SAVE_PFS = r33
+GR_SAVE_B0 = r34
+GR_SAVE_GP = r35
+GR_SAVE_r_sincos = r36
+
+
+RODATA
+
+// Pi/16 parts
.align 16
-double_sind_pi:
-ASM_TYPE_DIRECTIVE(double_sind_pi,@object)
-// data8 0xA2F9836E4E44152A, 0x00004001 // 16/pi (significand loaded w/ setf)
-// c90fdaa22168c234
- data8 0xC90FDAA22168C234, 0x00003FFC // pi/16 hi
-// c4c6628b80dc1cd1 29024e088a
- data8 0xC4C6628B80DC1CD1, 0x00003FBC // pi/16 lo
-ASM_SIZE_DIRECTIVE(double_sind_pi)
-
-double_sind_pq_k4:
-ASM_TYPE_DIRECTIVE(double_sind_pq_k4,@object)
+LOCAL_OBJECT_START(double_sincos_pi)
+ data8 0xC90FDAA22168C234, 0x00003FFC // pi/16 1st part
+ data8 0xC4C6628B80DC1CD1, 0x00003FBC // pi/16 2nd part
+ data8 0xA4093822299F31D0, 0x00003F7A // pi/16 3rd part
+LOCAL_OBJECT_END(double_sincos_pi)
+
+// Coefficients for polynomials
+LOCAL_OBJECT_START(double_sincos_pq_k4)
data8 0x3EC71C963717C63A // P4
data8 0x3EF9FFBA8F191AE6 // Q4
data8 0xBF2A01A00F4E11A8 // P3
@@ -275,125 +281,112 @@ ASM_TYPE_DIRECTIVE(double_sind_pq_k4,@object)
data8 0x3FA555555554DD45 // Q2
data8 0xBFC5555555555555 // P1
data8 0xBFDFFFFFFFFFFFFC // Q1
-ASM_SIZE_DIRECTIVE(double_sind_pq_k4)
+LOCAL_OBJECT_END(double_sincos_pq_k4)
+// Sincos table (S[m], C[m])
+LOCAL_OBJECT_START(double_sin_cos_beta_k4)
-double_sin_cos_beta_k4:
-ASM_TYPE_DIRECTIVE(double_sin_cos_beta_k4,@object)
data8 0x0000000000000000 , 0x00000000 // sin( 0 pi/16) S0
data8 0x8000000000000000 , 0x00003fff // cos( 0 pi/16) C0
-
+//
data8 0xc7c5c1e34d3055b3 , 0x00003ffc // sin( 1 pi/16) S1
data8 0xfb14be7fbae58157 , 0x00003ffe // cos( 1 pi/16) C1
-
+//
data8 0xc3ef1535754b168e , 0x00003ffd // sin( 2 pi/16) S2
data8 0xec835e79946a3146 , 0x00003ffe // cos( 2 pi/16) C2
-
+//
data8 0x8e39d9cd73464364 , 0x00003ffe // sin( 3 pi/16) S3
data8 0xd4db3148750d181a , 0x00003ffe // cos( 3 pi/16) C3
-
+//
data8 0xb504f333f9de6484 , 0x00003ffe // sin( 4 pi/16) S4
data8 0xb504f333f9de6484 , 0x00003ffe // cos( 4 pi/16) C4
-
-
+//
data8 0xd4db3148750d181a , 0x00003ffe // sin( 5 pi/16) C3
data8 0x8e39d9cd73464364 , 0x00003ffe // cos( 5 pi/16) S3
-
+//
data8 0xec835e79946a3146 , 0x00003ffe // sin( 6 pi/16) C2
data8 0xc3ef1535754b168e , 0x00003ffd // cos( 6 pi/16) S2
-
+//
data8 0xfb14be7fbae58157 , 0x00003ffe // sin( 7 pi/16) C1
data8 0xc7c5c1e34d3055b3 , 0x00003ffc // cos( 7 pi/16) S1
-
+//
data8 0x8000000000000000 , 0x00003fff // sin( 8 pi/16) C0
data8 0x0000000000000000 , 0x00000000 // cos( 8 pi/16) S0
-
-
+//
data8 0xfb14be7fbae58157 , 0x00003ffe // sin( 9 pi/16) C1
data8 0xc7c5c1e34d3055b3 , 0x0000bffc // cos( 9 pi/16) -S1
-
+//
data8 0xec835e79946a3146 , 0x00003ffe // sin(10 pi/16) C2
data8 0xc3ef1535754b168e , 0x0000bffd // cos(10 pi/16) -S2
-
+//
data8 0xd4db3148750d181a , 0x00003ffe // sin(11 pi/16) C3
data8 0x8e39d9cd73464364 , 0x0000bffe // cos(11 pi/16) -S3
-
+//
data8 0xb504f333f9de6484 , 0x00003ffe // sin(12 pi/16) S4
data8 0xb504f333f9de6484 , 0x0000bffe // cos(12 pi/16) -S4
-
-
+//
data8 0x8e39d9cd73464364 , 0x00003ffe // sin(13 pi/16) S3
data8 0xd4db3148750d181a , 0x0000bffe // cos(13 pi/16) -C3
-
+//
data8 0xc3ef1535754b168e , 0x00003ffd // sin(14 pi/16) S2
data8 0xec835e79946a3146 , 0x0000bffe // cos(14 pi/16) -C2
-
+//
data8 0xc7c5c1e34d3055b3 , 0x00003ffc // sin(15 pi/16) S1
data8 0xfb14be7fbae58157 , 0x0000bffe // cos(15 pi/16) -C1
-
+//
data8 0x0000000000000000 , 0x00000000 // sin(16 pi/16) S0
data8 0x8000000000000000 , 0x0000bfff // cos(16 pi/16) -C0
-
-
+//
data8 0xc7c5c1e34d3055b3 , 0x0000bffc // sin(17 pi/16) -S1
data8 0xfb14be7fbae58157 , 0x0000bffe // cos(17 pi/16) -C1
-
+//
data8 0xc3ef1535754b168e , 0x0000bffd // sin(18 pi/16) -S2
data8 0xec835e79946a3146 , 0x0000bffe // cos(18 pi/16) -C2
-
+//
data8 0x8e39d9cd73464364 , 0x0000bffe // sin(19 pi/16) -S3
data8 0xd4db3148750d181a , 0x0000bffe // cos(19 pi/16) -C3
-
+//
data8 0xb504f333f9de6484 , 0x0000bffe // sin(20 pi/16) -S4
data8 0xb504f333f9de6484 , 0x0000bffe // cos(20 pi/16) -S4
-
-
+//
data8 0xd4db3148750d181a , 0x0000bffe // sin(21 pi/16) -C3
data8 0x8e39d9cd73464364 , 0x0000bffe // cos(21 pi/16) -S3
-
+//
data8 0xec835e79946a3146 , 0x0000bffe // sin(22 pi/16) -C2
data8 0xc3ef1535754b168e , 0x0000bffd // cos(22 pi/16) -S2
-
+//
data8 0xfb14be7fbae58157 , 0x0000bffe // sin(23 pi/16) -C1
data8 0xc7c5c1e34d3055b3 , 0x0000bffc // cos(23 pi/16) -S1
-
+//
data8 0x8000000000000000 , 0x0000bfff // sin(24 pi/16) -C0
data8 0x0000000000000000 , 0x00000000 // cos(24 pi/16) S0
-
-
+//
data8 0xfb14be7fbae58157 , 0x0000bffe // sin(25 pi/16) -C1
data8 0xc7c5c1e34d3055b3 , 0x00003ffc // cos(25 pi/16) S1
-
+//
data8 0xec835e79946a3146 , 0x0000bffe // sin(26 pi/16) -C2
data8 0xc3ef1535754b168e , 0x00003ffd // cos(26 pi/16) S2
-
+//
data8 0xd4db3148750d181a , 0x0000bffe // sin(27 pi/16) -C3
data8 0x8e39d9cd73464364 , 0x00003ffe // cos(27 pi/16) S3
-
+//
data8 0xb504f333f9de6484 , 0x0000bffe // sin(28 pi/16) -S4
data8 0xb504f333f9de6484 , 0x00003ffe // cos(28 pi/16) S4
-
-
+//
data8 0x8e39d9cd73464364 , 0x0000bffe // sin(29 pi/16) -S3
data8 0xd4db3148750d181a , 0x00003ffe // cos(29 pi/16) C3
-
+//
data8 0xc3ef1535754b168e , 0x0000bffd // sin(30 pi/16) -S2
data8 0xec835e79946a3146 , 0x00003ffe // cos(30 pi/16) C2
-
+//
data8 0xc7c5c1e34d3055b3 , 0x0000bffc // sin(31 pi/16) -S1
data8 0xfb14be7fbae58157 , 0x00003ffe // cos(31 pi/16) C1
-
+//
data8 0x0000000000000000 , 0x00000000 // sin(32 pi/16) S0
data8 0x8000000000000000 , 0x00003fff // cos(32 pi/16) C0
-ASM_SIZE_DIRECTIVE(double_sin_cos_beta_k4)
+LOCAL_OBJECT_END(double_sin_cos_beta_k4)
-.align 32
-.global sin#
-.global cos#
-#ifdef _LIBC
-.global __sin#
-.global __cos#
-#endif
+.section .text
////////////////////////////////////////////////////////
// There are two entry points: sin and cos
@@ -402,3100 +395,374 @@ ASM_SIZE_DIRECTIVE(double_sin_cos_beta_k4)
// If from sin, p8 is true
// If from cos, p9 is true
-.section .text
-.proc sin#
-#ifdef _LIBC
-.proc __sin#
-#endif
-.align 32
-
-sin:
-#ifdef _LIBC
-__sin:
-#endif
+GLOBAL_IEEE754_ENTRY(sin)
{ .mlx
- alloc r32=ar.pfs,1,13,0,0
- movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
+ getf.exp sincos_r_signexp = f8
+ movl sincos_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // signd of 16/pi
}
{ .mlx
- addl sind_AD_1 = @ltoff(double_sind_pi), gp
- movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
+ addl sincos_AD_1 = @ltoff(double_sincos_pi), gp
+ movl sincos_GR_rshf_2to61 = 0x47b8000000000000 // 1.1 2^(63+63-2)
}
;;
{ .mfi
- ld8 sind_AD_1 = [sind_AD_1]
- fnorm sind_NORM_f8 = f8
- cmp.eq p8,p9 = r0, r0
+ ld8 sincos_AD_1 = [sincos_AD_1]
+ fnorm.s0 sincos_NORM_f8 = f8 // Normalize argument
+ cmp.eq p8,p9 = r0, r0 // set p8 (clear p9) for sin
}
{ .mib
- mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
- mov sind_r_sincos = 0x0
- br.cond.sptk L(SIND_SINCOS)
+ mov sincos_GR_exp_2tom61 = 0xffff-61 // exponent of scale 2^-61
+ mov sincos_r_sincos = 0x0 // sincos_r_sincos = 0 for sin
+ br.cond.sptk _SINCOS_COMMON // go to common part
}
;;
-.endp sin
-ASM_SIZE_DIRECTIVE(sin)
-
+GLOBAL_IEEE754_END(sin)
-.section .text
-.proc cos#
-#ifdef _LIBC
-.proc __cos#
-#endif
-.align 32
-cos:
-#ifdef _LIBC
-__cos:
-#endif
+GLOBAL_IEEE754_ENTRY(cos)
{ .mlx
- alloc r32=ar.pfs,1,13,0,0
- movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
+ getf.exp sincos_r_signexp = f8
+ movl sincos_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // signd of 16/pi
}
{ .mlx
- addl sind_AD_1 = @ltoff(double_sind_pi), gp
- movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
+ addl sincos_AD_1 = @ltoff(double_sincos_pi), gp
+ movl sincos_GR_rshf_2to61 = 0x47b8000000000000 // 1.1 2^(63+63-2)
}
;;
{ .mfi
- ld8 sind_AD_1 = [sind_AD_1]
- fnorm.s1 sind_NORM_f8 = f8
- cmp.eq p9,p8 = r0, r0
+ ld8 sincos_AD_1 = [sincos_AD_1]
+ fnorm.s1 sincos_NORM_f8 = f8 // Normalize argument
+ cmp.eq p9,p8 = r0, r0 // set p9 (clear p8) for cos
}
{ .mib
- mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
- mov sind_r_sincos = 0x8
- br.cond.sptk L(SIND_SINCOS)
+ mov sincos_GR_exp_2tom61 = 0xffff-61 // exp of scale 2^-61
+ mov sincos_r_sincos = 0x8 // sincos_r_sincos = 8 for cos
+ nop.b 999
}
;;
-
////////////////////////////////////////////////////////
// All entry points end up here.
-// If from sin, sind_r_sincos is 0 and p8 is true
-// If from cos, sind_r_sincos is 8 = 2^(k-1) and p9 is true
-// We add sind_r_sincos to N
+// If from sin, sincos_r_sincos is 0 and p8 is true
+// If from cos, sincos_r_sincos is 8 = 2^(k-1) and p9 is true
+// We add sincos_r_sincos to N
-L(SIND_SINCOS):
+///////////// Common sin and cos part //////////////////
+_SINCOS_COMMON:
// Form two constants we need
// 16/pi * 2^-2 * 2^63, scaled by 2^61 since we just loaded the significand
// 1.1000...000 * 2^(63+63-2) to right shift int(W) into the low significand
-// fcmp used to set denormal, and invalid on snans
{ .mfi
- setf.sig sind_SIG_INV_PI_BY_16_2TO61 = sind_GR_sig_inv_pi_by_16
- fcmp.eq.s0 p12,p0=f8,f0
- mov sind_r_17_ones = 0x1ffff
+ setf.sig sincos_SIG_INV_PI_BY_16_2TO61 = sincos_GR_sig_inv_pi_by_16
+ fclass.m p6,p0 = f8, 0xe7 // if x = 0,inf,nan
+ mov sincos_exp_limit = 0x1001a
}
{ .mlx
- setf.d sind_RSHF_2TO61 = sind_GR_rshf_2to61
- movl sind_GR_rshf = 0x43e8000000000000 // 1.1000 2^63 for right shift
-}
+ setf.d sincos_RSHF_2TO61 = sincos_GR_rshf_2to61
+ movl sincos_GR_rshf = 0x43e8000000000000 // 1.1 2^63
+} // Right shift
;;
// Form another constant
// 2^-61 for scaling Nfloat
-// 0x10009 is register_bias + 10.
-// So if f8 > 2^10 = Gamma, go to DBX
+// 0x1001a is register_bias + 27.
+// So if f8 >= 2^27, go to large argument routines
{ .mfi
- setf.exp sind_2TOM61 = sind_GR_exp_2tom61
- fclass.m p13,p0 = f8, 0x23 // Test for x inf
- mov sind_exp_limit = 0x10009
+ alloc r32 = ar.pfs, 1, 4, 0, 0
+ fclass.m p11,p0 = f8, 0x0b // Test for x=unorm
+ mov sincos_GR_all_ones = -1 // For "inexect" constant create
+}
+{ .mib
+ setf.exp sincos_2TOM61 = sincos_GR_exp_2tom61
+ nop.i 999
+(p6) br.cond.spnt _SINCOS_SPECIAL_ARGS
}
;;
// Load the two pieces of pi/16
// Form another constant
// 1.1000...000 * 2^63, the right shift constant
-{ .mmf
- ldfe sind_Pi_by_16_hi = [sind_AD_1],16
- setf.d sind_RSHF = sind_GR_rshf
- fclass.m p14,p0 = f8, 0xc3 // Test for x nan
+{ .mmb
+ ldfe sincos_Pi_by_16_1 = [sincos_AD_1],16
+ setf.d sincos_RSHF = sincos_GR_rshf
+(p11) br.cond.spnt _SINCOS_UNORM // Branch if x=unorm
}
;;
-{ .mfi
- ldfe sind_Pi_by_16_lo = [sind_AD_1],16
-(p13) frcpa.s0 f8,p12=f0,f0 // force qnan indef for x=inf
- addl gr_tmp = -1,r0
-}
-{ .mfb
- addl sind_AD_beta_table = @ltoff(double_sin_cos_beta_k4), gp
- nop.f 999
-(p13) br.ret.spnt b0 ;; // Exit for x=inf
-}
-
-// Start loading P, Q coefficients
-// SIN(0)
-{ .mfi
- ldfpd sind_P4,sind_Q4 = [sind_AD_1],16
-(p8) fclass.m.unc p6,p0 = f8, 0x07 // Test for sin(0)
- nop.i 999
-}
-{ .mfb
- addl sind_AD_beta_table = @ltoff(double_sin_cos_beta_k4), gp
-(p14) fma.d f8=f8,f1,f0 // qnan for x=nan
-(p14) br.ret.spnt b0 ;; // Exit for x=nan
-}
-
+_SINCOS_COMMON2:
+// Return here if x=unorm
+// Create constant used to set inexact
+{ .mmi
+ ldfe sincos_Pi_by_16_2 = [sincos_AD_1],16
+ setf.sig fp_tmp = sincos_GR_all_ones
+ nop.i 999
+};;
-// COS(0)
+// Select exponent (17 lsb)
{ .mfi
- getf.exp sind_r_signexp = f8
-(p9) fclass.m.unc p7,p0 = f8, 0x07 // Test for sin(0)
- nop.i 999
-}
-{ .mfi
- ld8 sind_AD_beta_table = [sind_AD_beta_table]
- nop.f 999
- nop.i 999 ;;
-}
-
-{ .mmb
- ldfpd sind_P3,sind_Q3 = [sind_AD_1],16
- setf.sig fp_tmp = gr_tmp // Create constant such that fmpy sets inexact
-(p6) br.ret.spnt b0 ;;
-}
-
-{ .mfb
- and sind_r_exp = sind_r_17_ones, sind_r_signexp
-(p7) fmerge.s f8 = f1,f1
-(p7) br.ret.spnt b0 ;;
-}
+ ldfe sincos_Pi_by_16_3 = [sincos_AD_1],16
+ nop.f 999
+ dep.z sincos_r_exp = sincos_r_signexp, 0, 17
+};;
+// Polynomial coefficients (Q4, P4, Q3, P3, Q2, Q1, P2, P1) loading
// p10 is true if we must call routines to handle larger arguments
-// p10 is true if f8 exp is > 0x10009
-
-{ .mfi
- ldfpd sind_P2,sind_Q2 = [sind_AD_1],16
- nop.f 999
- cmp.ge p10,p0 = sind_r_exp,sind_exp_limit
-}
-;;
+// p10 is true if f8 exp is >= 0x1001a (2^27)
+{ .mmb
+ ldfpd sincos_P4,sincos_Q4 = [sincos_AD_1],16
+ cmp.ge p10,p0 = sincos_r_exp,sincos_exp_limit
+(p10) br.cond.spnt _SINCOS_LARGE_ARGS // Go to "large args" routine
+};;
-// sind_W = x * sind_Inv_Pi_by_16
+// sincos_W = x * sincos_Inv_Pi_by_16
// Multiply x by scaled 16/pi and add large const to shift integer part of W to
// rightmost bits of significand
{ .mfi
- ldfpd sind_P1,sind_Q1 = [sind_AD_1]
- fma.s1 sind_W_2TO61_RSH = sind_NORM_f8,sind_SIG_INV_PI_BY_16_2TO61,sind_RSHF_2TO61
- nop.i 999
-}
-{ .mbb
-(p10) cmp.ne.unc p11,p12=sind_r_sincos,r0 // p11 call __libm_cos_double_dbx
- // p12 call __libm_sin_double_dbx
-(p11) br.cond.spnt L(COSD_DBX)
-(p12) br.cond.spnt L(SIND_DBX)
-}
-;;
-
+ ldfpd sincos_P3,sincos_Q3 = [sincos_AD_1],16
+ fma.s1 sincos_W_2TO61_RSH = sincos_NORM_f8,sincos_SIG_INV_PI_BY_16_2TO61,sincos_RSHF_2TO61
+ nop.i 999
+};;
-// sind_NFLOAT = Round_Int_Nearest(sind_W)
+// get N = (int)sincos_int_Nfloat
+// sincos_NFLOAT = Round_Int_Nearest(sincos_W)
// This is done by scaling back by 2^-61 and subtracting the shift constant
-{ .mfi
- nop.m 999
- fms.s1 sind_NFLOAT = sind_W_2TO61_RSH,sind_2TOM61,sind_RSHF
- nop.i 999 ;;
-}
-
-
-// get N = (int)sind_int_Nfloat
-{ .mfi
- getf.sig sind_GR_n = sind_W_2TO61_RSH
- nop.f 999
- nop.i 999 ;;
-}
-
-// Add 2^(k-1) (which is in sind_r_sincos) to N
-// sind_r = -sind_Nfloat * sind_Pi_by_16_hi + x
-// sind_r = sind_r -sind_Nfloat * sind_Pi_by_16_lo
-{ .mfi
- add sind_GR_n = sind_GR_n, sind_r_sincos
- fnma.s1 sind_r = sind_NFLOAT, sind_Pi_by_16_hi, sind_NORM_f8
- nop.i 999 ;;
-}
-
-
-// Get M (least k+1 bits of N)
-{ .mmi
- and sind_GR_m = 0x1f,sind_GR_n ;;
- nop.m 999
- shl sind_GR_32m = sind_GR_m,5 ;;
-}
-
-// Add 32*M to address of sin_cos_beta table
-{ .mmi
- add sind_AD_2 = sind_GR_32m, sind_AD_beta_table
- nop.m 999
- nop.i 999 ;;
-}
-
-{ .mfi
- ldfe sind_Sm = [sind_AD_2],16
-(p8) fclass.m.unc p10,p0=f8,0x0b // If sin, note denormal input to set uflow
- nop.i 999 ;;
-}
-
-{ .mfi
- ldfe sind_Cm = [sind_AD_2]
- fnma.s1 sind_r = sind_NFLOAT, sind_Pi_by_16_lo, sind_r
- nop.i 999 ;;
-}
-
-// get rsq
-{ .mfi
- nop.m 999
- fma.s1 sind_rsq = sind_r, sind_r, f0
- nop.i 999
-}
-{ .mfi
- nop.m 999
- fmpy.s0 fp_tmp = fp_tmp,fp_tmp // fmpy forces inexact flag
- nop.i 999 ;;
-}
-
-// form P and Q series
-{ .mfi
- nop.m 999
- fma.s1 sind_P_temp1 = sind_rsq, sind_P4, sind_P3
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_Q_temp1 = sind_rsq, sind_Q4, sind_Q3
- nop.i 999 ;;
-}
-
-// get rcube and sm*rsq
-{ .mfi
- nop.m 999
- fmpy.s1 sind_srsq = sind_Sm,sind_rsq
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fmpy.s1 sind_rcub = sind_r, sind_rsq
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_Q_temp2 = sind_rsq, sind_Q_temp1, sind_Q2
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_P_temp2 = sind_rsq, sind_P_temp1, sind_P2
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_Q = sind_rsq, sind_Q_temp2, sind_Q1
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_P = sind_rsq, sind_P_temp2, sind_P1
- nop.i 999 ;;
-}
-
-// Get final P and Q
-{ .mfi
- nop.m 999
- fma.s1 sind_Q = sind_srsq,sind_Q, sind_Sm
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fma.s1 sind_P = sind_rcub,sind_P, sind_r
- nop.i 999 ;;
-}
-
-// If sin(denormal), force inexact to be set
-{ .mfi
- nop.m 999
-(p10) fmpy.d.s0 fp_tmp = f8,f8
- nop.i 999 ;;
-}
-
-// Final calculation
-{ .mfb
- nop.m 999
- fma.d f8 = sind_Cm, sind_P, sind_Q
- br.ret.sptk b0 ;;
-}
-.endp cos#
-ASM_SIZE_DIRECTIVE(cos#)
-
-
-
-.proc __libm_callout_1s
-__libm_callout_1s:
-L(SIND_DBX):
-.prologue
-{ .mfi
- nop.m 0
- nop.f 0
-.save ar.pfs,GR_SAVE_PFS
- mov GR_SAVE_PFS=ar.pfs
-}
-;;
-
-{ .mfi
- mov GR_SAVE_GP=gp
- nop.f 0
-.save b0, GR_SAVE_B0
- mov GR_SAVE_B0=b0
-}
-
-.body
-{ .mib
- nop.m 999
- nop.i 999
- br.call.sptk.many b0=__libm_sin_double_dbx# ;;
-}
-;;
-
-
-{ .mfi
- mov gp = GR_SAVE_GP
- nop.f 999
- mov b0 = GR_SAVE_B0
-}
-;;
-
-{ .mib
- nop.m 999
- mov ar.pfs = GR_SAVE_PFS
- br.ret.sptk b0 ;;
-}
-.endp __libm_callout_1s
-ASM_SIZE_DIRECTIVE(__libm_callout_1s)
-
-
-.proc __libm_callout_1c
-__libm_callout_1c:
-L(COSD_DBX):
-.prologue
-{ .mfi
- nop.m 0
- nop.f 0
-.save ar.pfs,GR_SAVE_PFS
- mov GR_SAVE_PFS=ar.pfs
-}
-;;
-
-{ .mfi
- mov GR_SAVE_GP=gp
- nop.f 0
-.save b0, GR_SAVE_B0
- mov GR_SAVE_B0=b0
-}
-
-.body
-{ .mib
- nop.m 999
- nop.i 999
- br.call.sptk.many b0=__libm_cos_double_dbx# ;;
-}
-;;
-
-
-{ .mfi
- mov gp = GR_SAVE_GP
- nop.f 999
- mov b0 = GR_SAVE_B0
-}
-;;
-
-{ .mib
- nop.m 999
- mov ar.pfs = GR_SAVE_PFS
- br.ret.sptk b0 ;;
-}
-.endp __libm_callout_1c
-ASM_SIZE_DIRECTIVE(__libm_callout_1c)
-
-
-// ====================================================================
-// ====================================================================
-
-// These functions calculate the sin and cos for inputs
-// greater than 2^10
-// __libm_sin_double_dbx# and __libm_cos_double_dbx#
-
-// *********************************************************************
-// *********************************************************************
-//
-// Function: Combined sin(x) and cos(x), where
-//
-// sin(x) = sine(x), for double precision x values
-// cos(x) = cosine(x), for double precision x values
-//
-// *********************************************************************
-//
-// Accuracy: Within .7 ulps for 80-bit floating point values
-// Very accurate for double precision values
-//
-// *********************************************************************
-//
-// Resources Used:
-//
-// Floating-Point Registers: f8 (Input and Return Value)
-// f32-f99
-//
-// General Purpose Registers:
-// r32-r43
-// r44-r45 (Used to pass arguments to pi_by_2 reduce routine)
-//
-// Predicate Registers: p6-p13
-//
-// *********************************************************************
-//
-// IEEE Special Conditions:
-//
-// Denormal fault raised on denormal inputs
-// Overflow exceptions do not occur
-// Underflow exceptions raised when appropriate for sin
-// (No specialized error handling for this routine)
-// Inexact raised when appropriate by algorithm
-//
-// sin(SNaN) = QNaN
-// sin(QNaN) = QNaN
-// sin(inf) = QNaN
-// sin(+/-0) = +/-0
-// cos(inf) = QNaN
-// cos(SNaN) = QNaN
-// cos(QNaN) = QNaN
-// cos(0) = 1
-//
-// *********************************************************************
-//
-// Mathematical Description
-// ========================
-//
-// The computation of FSIN and FCOS is best handled in one piece of
-// code. The main reason is that given any argument Arg, computation
-// of trigonometric functions first calculate N and an approximation
-// to alpha where
-//
-// Arg = N pi/2 + alpha, |alpha| <= pi/4.
-//
-// Since
-//
-// cos( Arg ) = sin( (N+1) pi/2 + alpha ),
-//
-// therefore, the code for computing sine will produce cosine as long
-// as 1 is added to N immediately after the argument reduction
-// process.
-//
-// Let M = N if sine
-// N+1 if cosine.
-//
-// Now, given
-//
-// Arg = M pi/2 + alpha, |alpha| <= pi/4,
-//
-// let I = M mod 4, or I be the two lsb of M when M is represented
-// as 2's complement. I = [i_0 i_1]. Then
-//
-// sin( Arg ) = (-1)^i_0 sin( alpha ) if i_1 = 0,
-// = (-1)^i_0 cos( alpha ) if i_1 = 1.
-//
-// For example:
-// if M = -1, I = 11
-// sin ((-pi/2 + alpha) = (-1) cos (alpha)
-// if M = 0, I = 00
-// sin (alpha) = sin (alpha)
-// if M = 1, I = 01
-// sin (pi/2 + alpha) = cos (alpha)
-// if M = 2, I = 10
-// sin (pi + alpha) = (-1) sin (alpha)
-// if M = 3, I = 11
-// sin ((3/2)pi + alpha) = (-1) cos (alpha)
-//
-// The value of alpha is obtained by argument reduction and
-// represented by two working precision numbers r and c where
-//
-// alpha = r + c accurately.
-//
-// The reduction method is described in a previous write up.
-// The argument reduction scheme identifies 4 cases. For Cases 2
-// and 4, because |alpha| is small, sin(r+c) and cos(r+c) can be
-// computed very easily by 2 or 3 terms of the Taylor series
-// expansion as follows:
-//
-// Case 2:
-// -------
-//
-// sin(r + c) = r + c - r^3/6 accurately
-// cos(r + c) = 1 - 2^(-67) accurately
-//
-// Case 4:
-// -------
-//
-// sin(r + c) = r + c - r^3/6 + r^5/120 accurately
-// cos(r + c) = 1 - r^2/2 + r^4/24 accurately
-//
-// The only cases left are Cases 1 and 3 of the argument reduction
-// procedure. These two cases will be merged since after the
-// argument is reduced in either cases, we have the reduced argument
-// represented as r + c and that the magnitude |r + c| is not small
-// enough to allow the usage of a very short approximation.
-//
-// The required calculation is either
-//
-// sin(r + c) = sin(r) + correction, or
-// cos(r + c) = cos(r) + correction.
-//
-// Specifically,
-//
-// sin(r + c) = sin(r) + c sin'(r) + O(c^2)
-// = sin(r) + c cos (r) + O(c^2)
-// = sin(r) + c(1 - r^2/2) accurately.
-// Similarly,
-//
-// cos(r + c) = cos(r) - c sin(r) + O(c^2)
-// = cos(r) - c(r - r^3/6) accurately.
-//
-// We therefore concentrate on accurately calculating sin(r) and
-// cos(r) for a working-precision number r, |r| <= pi/4 to within
-// 0.1% or so.
-//
-// The greatest challenge of this task is that the second terms of
-// the Taylor series
-//
-// r - r^3/3! + r^r/5! - ...
-//
-// and
-//
-// 1 - r^2/2! + r^4/4! - ...
-//
-// are not very small when |r| is close to pi/4 and the rounding
-// errors will be a concern if simple polynomial accumulation is
-// used. When |r| < 2^-3, however, the second terms will be small
-// enough (6 bits or so of right shift) that a normal Horner
-// recurrence suffices. Hence there are two cases that we consider
-// in the accurate computation of sin(r) and cos(r), |r| <= pi/4.
-//
-// Case small_r: |r| < 2^(-3)
-// --------------------------
-//
-// Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1],
-// we have
-//
-// sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
-// = (-1)^i_0 * cos(r + c) if i_1 = 1
-//
-// can be accurately approximated by
-//
-// sin(Arg) = (-1)^i_0 * [sin(r) + c] if i_1 = 0
-// = (-1)^i_0 * [cos(r) - c*r] if i_1 = 1
-//
-// because |r| is small and thus the second terms in the correction
-// are unneccessary.
-//
-// Finally, sin(r) and cos(r) are approximated by polynomials of
-// moderate lengths.
-//
-// sin(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11
-// cos(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10
-//
-// We can make use of predicates to selectively calculate
-// sin(r) or cos(r) based on i_1.
-//
-// Case normal_r: 2^(-3) <= |r| <= pi/4
-// ------------------------------------
-//
-// This case is more likely than the previous one if one considers
-// r to be uniformly distributed in [-pi/4 pi/4]. Again,
-//
-// sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
-// = (-1)^i_0 * cos(r + c) if i_1 = 1.
-//
-// Because |r| is now larger, we need one extra term in the
-// correction. sin(Arg) can be accurately approximated by
-//
-// sin(Arg) = (-1)^i_0 * [sin(r) + c(1-r^2/2)] if i_1 = 0
-// = (-1)^i_0 * [cos(r) - c*r*(1 - r^2/6)] i_1 = 1.
-//
-// Finally, sin(r) and cos(r) are approximated by polynomials of
-// moderate lengths.
-//
-// sin(r) = r + PP_1_hi r^3 + PP_1_lo r^3 +
-// PP_2 r^5 + ... + PP_8 r^17
-//
-// cos(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16
-//
-// where PP_1_hi is only about 16 bits long and QQ_1 is -1/2.
-// The crux in accurate computation is to calculate
-//
-// r + PP_1_hi r^3 or 1 + QQ_1 r^2
-//
-// accurately as two pieces: U_hi and U_lo. The way to achieve this
-// is to obtain r_hi as a 10 sig. bit number that approximates r to
-// roughly 8 bits or so of accuracy. (One convenient way is
-//
-// r_hi := frcpa( frcpa( r ) ).)
-//
-// This way,
-//
-// r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 +
-// PP_1_hi (r^3 - r_hi^3)
-// = [r + PP_1_hi r_hi^3] +
-// [PP_1_hi (r - r_hi)
-// (r^2 + r_hi r + r_hi^2) ]
-// = U_hi + U_lo
-//
-// Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long,
-// PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed
-// exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign
-// and that there is no more than 8 bit shift off between r and
-// PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus
-// calculated without any error. Finally, the fact that
-//
-// |U_lo| <= 2^(-8) |U_hi|
-//
-// says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly
-// 8 extra bits of accuracy.
-//
-// Similarly,
-//
-// 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] +
-// [QQ_1 (r - r_hi)(r + r_hi)]
-// = U_hi + U_lo.
-//
-// Summarizing, we calculate r_hi = frcpa( frcpa( r ) ).
-//
-// If i_1 = 0, then
-//
-// U_hi := r + PP_1_hi * r_hi^3
-// U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2)
-// poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17
-// correction := c * ( 1 + C_1 r^2 )
-//
-// Else ...i_1 = 1
-//
-// U_hi := 1 + QQ_1 * r_hi * r_hi
-// U_lo := QQ_1 * (r - r_hi) * (r + r_hi)
-// poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16
-// correction := -c * r * (1 + S_1 * r^2)
-//
-// End
-//
-// Finally,
-//
-// V := poly + ( U_lo + correction )
-//
-// / U_hi + V if i_0 = 0
-// result := |
-// \ (-U_hi) - V if i_0 = 1
-//
-// It is important that in the last step, negation of U_hi is
-// performed prior to the subtraction which is to be performed in
-// the user-set rounding mode.
-//
-//
-// Algorithmic Description
-// =======================
-//
-// The argument reduction algorithm is tightly integrated into FSIN
-// and FCOS which share the same code. The following is complete and
-// self-contained. The argument reduction description given
-// previously is repeated below.
-//
-//
-// Step 0. Initialization.
-//
-// If FSIN is invoked, set N_inc := 0; else if FCOS is invoked,
-// set N_inc := 1.
-//
-// Step 1. Check for exceptional and special cases.
-//
-// * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special
-// handling.
-// * If |Arg| < 2^24, go to Step 2 for reduction of moderate
-// arguments. This is the most likely case.
-// * If |Arg| < 2^63, go to Step 8 for pre-reduction of large
-// arguments.
-// * If |Arg| >= 2^63, go to Step 10 for special handling.
-//
-// Step 2. Reduction of moderate arguments.
-//
-// If |Arg| < pi/4 ...quick branch
-// N_fix := N_inc (integer)
-// r := Arg
-// c := 0.0
-// Branch to Step 4, Case_1_complete
-// Else ...cf. argument reduction
-// N := Arg * two_by_PI (fp)
-// N_fix := fcvt.fx( N ) (int)
-// N := fcvt.xf( N_fix )
-// N_fix := N_fix + N_inc
-// s := Arg - N * P_1 (first piece of pi/2)
-// w := -N * P_2 (second piece of pi/2)
-//
-// If |s| >= 2^(-33)
-// go to Step 3, Case_1_reduce
-// Else
-// go to Step 7, Case_2_reduce
-// Endif
-// Endif
-//
-// Step 3. Case_1_reduce.
-//
-// r := s + w
-// c := (s - r) + w ...observe order
-//
-// Step 4. Case_1_complete
-//
-// ...At this point, the reduced argument alpha is
-// ...accurately represented as r + c.
-// If |r| < 2^(-3), go to Step 6, small_r.
-//
-// Step 5. Normal_r.
-//
-// Let [i_0 i_1] by the 2 lsb of N_fix.
-// FR_rsq := r * r
-// r_hi := frcpa( frcpa( r ) )
-// r_lo := r - r_hi
-//
-// If i_1 = 0, then
-// poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8))
-// U_hi := r + PP_1_hi*r_hi*r_hi*r_hi ...any order
-// U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi)
-// correction := c + c*C_1*FR_rsq ...any order
-// Else
-// poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8))
-// U_hi := 1 + QQ_1 * r_hi * r_hi ...any order
-// U_lo := QQ_1 * r_lo * (r + r_hi)
-// correction := -c*(r + S_1*FR_rsq*r) ...any order
-// Endif
-//
-// V := poly + (U_lo + correction) ...observe order
-//
-// result := (i_0 == 0? 1.0 : -1.0)
-//
-// Last instruction in user-set rounding mode
-//
-// result := (i_0 == 0? result*U_hi + V :
-// result*U_hi - V)
-//
-// Return
-//
-// Step 6. Small_r.
-//
-// ...Use flush to zero mode without causing exception
-// Let [i_0 i_1] be the two lsb of N_fix.
-//
-// FR_rsq := r * r
-//
-// If i_1 = 0 then
-// z := FR_rsq*FR_rsq; z := FR_rsq*z *r
-// poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5)
-// poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2)
-// correction := c
-// result := r
-// Else
-// z := FR_rsq*FR_rsq; z := FR_rsq*z
-// poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5)
-// poly_hi := FR_rsq*(C_1 + FR_rsq*C_2)
-// correction := -c*r
-// result := 1
-// Endif
-//
-// poly := poly_hi + (z * poly_lo + correction)
-//
-// If i_0 = 1, result := -result
-//
-// Last operation. Perform in user-set rounding mode
-//
-// result := (i_0 == 0? result + poly :
-// result - poly )
-// Return
-//
-// Step 7. Case_2_reduce.
-//
-// ...Refer to the write up for argument reduction for
-// ...rationale. The reduction algorithm below is taken from
-// ...argument reduction description and integrated this.
-//
-// w := N*P_3
-// U_1 := N*P_2 + w ...FMA
-// U_2 := (N*P_2 - U_1) + w ...2 FMA
-// ...U_1 + U_2 is N*(P_2+P_3) accurately
-//
-// r := s - U_1
-// c := ( (s - r) - U_1 ) - U_2
-//
-// ...The mathematical sum r + c approximates the reduced
-// ...argument accurately. Note that although compared to
-// ...Case 1, this case requires much more work to reduce
-// ...the argument, the subsequent calculation needed for
-// ...any of the trigonometric function is very little because
-// ...|alpha| < 1.01*2^(-33) and thus two terms of the
-// ...Taylor series expansion suffices.
-//
-// If i_1 = 0 then
-// poly := c + S_1 * r * r * r ...any order
-// result := r
-// Else
-// poly := -2^(-67)
-// result := 1.0
-// Endif
-//
-// If i_0 = 1, result := -result
-//
-// Last operation. Perform in user-set rounding mode
-//
-// result := (i_0 == 0? result + poly :
-// result - poly )
-//
-// Return
-//
-//
-// Step 8. Pre-reduction of large arguments.
-//
-// ...Again, the following reduction procedure was described
-// ...in the separate write up for argument reduction, which
-// ...is tightly integrated here.
-
-// N_0 := Arg * Inv_P_0
-// N_0_fix := fcvt.fx( N_0 )
-// N_0 := fcvt.xf( N_0_fix)
-
-// Arg' := Arg - N_0 * P_0
-// w := N_0 * d_1
-// N := Arg' * two_by_PI
-// N_fix := fcvt.fx( N )
-// N := fcvt.xf( N_fix )
-// N_fix := N_fix + N_inc
-//
-// s := Arg' - N * P_1
-// w := w - N * P_2
-//
-// If |s| >= 2^(-14)
-// go to Step 3
-// Else
-// go to Step 9
-// Endif
-//
-// Step 9. Case_4_reduce.
-//
-// ...first obtain N_0*d_1 and -N*P_2 accurately
-// U_hi := N_0 * d_1 V_hi := -N*P_2
-// U_lo := N_0 * d_1 - U_hi V_lo := -N*P_2 - U_hi ...FMAs
-//
-// ...compute the contribution from N_0*d_1 and -N*P_3
-// w := -N*P_3
-// w := w + N_0*d_2
-// t := U_lo + V_lo + w ...any order
-//
-// ...at this point, the mathematical value
-// ...s + U_hi + V_hi + t approximates the true reduced argument
-// ...accurately. Just need to compute this accurately.
-//
-// ...Calculate U_hi + V_hi accurately:
-// A := U_hi + V_hi
-// if |U_hi| >= |V_hi| then
-// a := (U_hi - A) + V_hi
-// else
-// a := (V_hi - A) + U_hi
-// endif
-// ...order in computing "a" must be observed. This branch is
-// ...best implemented by predicates.
-// ...A + a is U_hi + V_hi accurately. Moreover, "a" is
-// ...much smaller than A: |a| <= (1/2)ulp(A).
-//
-// ...Just need to calculate s + A + a + t
-// C_hi := s + A t := t + a
-// C_lo := (s - C_hi) + A
-// C_lo := C_lo + t
-//
-// ...Final steps for reduction
-// r := C_hi + C_lo
-// c := (C_hi - r) + C_lo
-//
-// ...At this point, we have r and c
-// ...And all we need is a couple of terms of the corresponding
-// ...Taylor series.
-//
-// If i_1 = 0
-// poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2)
-// result := r
-// Else
-// poly := FR_rsq*(C_1 + FR_rsq*C_2)
-// result := 1
-// Endif
-//
-// If i_0 = 1, result := -result
-//
-// Last operation. Perform in user-set rounding mode
-//
-// result := (i_0 == 0? result + poly :
-// result - poly )
-// Return
-//
-// Large Arguments: For arguments above 2**63, a Payne-Hanek
-// style argument reduction is used and pi_by_2 reduce is called.
-//
-
-
-#ifdef _LIBC
-.rodata
-#else
-.data
-#endif
-.align 64
-
-FSINCOS_CONSTANTS:
-ASM_TYPE_DIRECTIVE(FSINCOS_CONSTANTS,@object)
-data4 0x4B800000, 0xCB800000, 0x00000000,0x00000000 // two**24, -two**24
-data4 0x4E44152A, 0xA2F9836E, 0x00003FFE,0x00000000 // Inv_pi_by_2
-data4 0xCE81B9F1, 0xC84D32B0, 0x00004016,0x00000000 // P_0
-data4 0x2168C235, 0xC90FDAA2, 0x00003FFF,0x00000000 // P_1
-data4 0xFC8F8CBB, 0xECE675D1, 0x0000BFBD,0x00000000 // P_2
-data4 0xACC19C60, 0xB7ED8FBB, 0x0000BF7C,0x00000000 // P_3
-data4 0x5F000000, 0xDF000000, 0x00000000,0x00000000 // two_to_63, -two_to_63
-data4 0x6EC6B45A, 0xA397E504, 0x00003FE7,0x00000000 // Inv_P_0
-data4 0xDBD171A1, 0x8D848E89, 0x0000BFBF,0x00000000 // d_1
-data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C,0x00000000 // d_2
-data4 0x2168C234, 0xC90FDAA2, 0x00003FFE,0x00000000 // pi_by_4
-data4 0x2168C234, 0xC90FDAA2, 0x0000BFFE,0x00000000 // neg_pi_by_4
-data4 0x3E000000, 0xBE000000, 0x00000000,0x00000000 // two**-3, -two**-3
-data4 0x2F000000, 0xAF000000, 0x9E000000,0x00000000 // two**-33, -two**-33, -two**-67
-data4 0xA21C0BC9, 0xCC8ABEBC, 0x00003FCE,0x00000000 // PP_8
-data4 0x720221DA, 0xD7468A05, 0x0000BFD6,0x00000000 // PP_7
-data4 0x640AD517, 0xB092382F, 0x00003FDE,0x00000000 // PP_6
-data4 0xD1EB75A4, 0xD7322B47, 0x0000BFE5,0x00000000 // PP_5
-data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1
-data4 0x00000000, 0xAAAA0000, 0x0000BFFC,0x00000000 // PP_1_hi
-data4 0xBAF69EEA, 0xB8EF1D2A, 0x00003FEC,0x00000000 // PP_4
-data4 0x0D03BB69, 0xD00D00D0, 0x0000BFF2,0x00000000 // PP_3
-data4 0x88888962, 0x88888888, 0x00003FF8,0x00000000 // PP_2
-data4 0xAAAB0000, 0xAAAAAAAA, 0x0000BFEC,0x00000000 // PP_1_lo
-data4 0xC2B0FE52, 0xD56232EF, 0x00003FD2,0x00000000 // QQ_8
-data4 0x2B48DCA6, 0xC9C99ABA, 0x0000BFDA,0x00000000 // QQ_7
-data4 0x9C716658, 0x8F76C650, 0x00003FE2,0x00000000 // QQ_6
-data4 0xFDA8D0FC, 0x93F27DBA, 0x0000BFE9,0x00000000 // QQ_5
-data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1
-data4 0x00000000, 0x80000000, 0x0000BFFE,0x00000000 // QQ_1
-data4 0x0C6E5041, 0xD00D00D0, 0x00003FEF,0x00000000 // QQ_4
-data4 0x0B607F60, 0xB60B60B6, 0x0000BFF5,0x00000000 // QQ_3
-data4 0xAAAAAA9B, 0xAAAAAAAA, 0x00003FFA,0x00000000 // QQ_2
-data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1
-data4 0xAAAA719F, 0xAAAAAAAA, 0x00003FFA,0x00000000 // C_2
-data4 0x0356F994, 0xB60B60B6, 0x0000BFF5,0x00000000 // C_3
-data4 0xB2385EA9, 0xD00CFFD5, 0x00003FEF,0x00000000 // C_4
-data4 0x292A14CD, 0x93E4BD18, 0x0000BFE9,0x00000000 // C_5
-data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1
-data4 0x888868DB, 0x88888888, 0x00003FF8,0x00000000 // S_2
-data4 0x055EFD4B, 0xD00D00D0, 0x0000BFF2,0x00000000 // S_3
-data4 0x839730B9, 0xB8EF1C5D, 0x00003FEC,0x00000000 // S_4
-data4 0xE5B3F492, 0xD71EA3A4, 0x0000BFE5,0x00000000 // S_5
-data4 0x38800000, 0xB8800000, 0x00000000 // two**-14, -two**-14
-ASM_SIZE_DIRECTIVE(FSINCOS_CONSTANTS)
-
-FR_Input_X = f8
-FR_Neg_Two_to_M3 = f32
-FR_Two_to_63 = f32
-FR_Two_to_24 = f33
-FR_Pi_by_4 = f33
-FR_Two_to_M14 = f34
-FR_Two_to_M33 = f35
-FR_Neg_Two_to_24 = f36
-FR_Neg_Pi_by_4 = f36
-FR_Neg_Two_to_M14 = f37
-FR_Neg_Two_to_M33 = f38
-FR_Neg_Two_to_M67 = f39
-FR_Inv_pi_by_2 = f40
-FR_N_float = f41
-FR_N_fix = f42
-FR_P_1 = f43
-FR_P_2 = f44
-FR_P_3 = f45
-FR_s = f46
-FR_w = f47
-FR_c = f48
-FR_r = f49
-FR_Z = f50
-FR_A = f51
-FR_a = f52
-FR_t = f53
-FR_U_1 = f54
-FR_U_2 = f55
-FR_C_1 = f56
-FR_C_2 = f57
-FR_C_3 = f58
-FR_C_4 = f59
-FR_C_5 = f60
-FR_S_1 = f61
-FR_S_2 = f62
-FR_S_3 = f63
-FR_S_4 = f64
-FR_S_5 = f65
-FR_poly_hi = f66
-FR_poly_lo = f67
-FR_r_hi = f68
-FR_r_lo = f69
-FR_rsq = f70
-FR_r_cubed = f71
-FR_C_hi = f72
-FR_N_0 = f73
-FR_d_1 = f74
-FR_V = f75
-FR_V_hi = f75
-FR_V_lo = f76
-FR_U_hi = f77
-FR_U_lo = f78
-FR_U_hiabs = f79
-FR_V_hiabs = f80
-FR_PP_8 = f81
-FR_QQ_8 = f81
-FR_PP_7 = f82
-FR_QQ_7 = f82
-FR_PP_6 = f83
-FR_QQ_6 = f83
-FR_PP_5 = f84
-FR_QQ_5 = f84
-FR_PP_4 = f85
-FR_QQ_4 = f85
-FR_PP_3 = f86
-FR_QQ_3 = f86
-FR_PP_2 = f87
-FR_QQ_2 = f87
-FR_QQ_1 = f88
-FR_N_0_fix = f89
-FR_Inv_P_0 = f90
-FR_corr = f91
-FR_poly = f92
-FR_d_2 = f93
-FR_Two_to_M3 = f94
-FR_Neg_Two_to_63 = f94
-FR_P_0 = f95
-FR_C_lo = f96
-FR_PP_1 = f97
-FR_PP_1_lo = f98
-FR_ArgPrime = f99
-
-GR_Table_Base = r32
-GR_Table_Base1 = r33
-GR_i_0 = r34
-GR_i_1 = r35
-GR_N_Inc = r36
-GR_Sin_or_Cos = r37
-
-GR_SAVE_B0 = r39
-GR_SAVE_GP = r40
-GR_SAVE_PFS = r41
-
-.section .text
-.proc __libm_sin_double_dbx#
-.align 64
-__libm_sin_double_dbx:
-
-{ .mlx
-alloc GR_Table_Base = ar.pfs,0,12,2,0
- movl GR_Sin_or_Cos = 0x0 ;;
-}
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
-
-{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
-;;
-
-
-{ .mib
- nop.m 999
- nop.i 999
- br.cond.sptk L(SINCOS_CONTINUE) ;;
-}
-
-.endp __libm_sin_double_dbx#
-ASM_SIZE_DIRECTIVE(__libm_sin_double_dbx)
-
-.section .text
-.proc __libm_cos_double_dbx#
-__libm_cos_double_dbx:
-
-{ .mlx
-alloc GR_Table_Base= ar.pfs,0,12,2,0
- movl GR_Sin_or_Cos = 0x1 ;;
-}
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
-
-{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
-;;
-
-//
-// Load Table Address
-//
-L(SINCOS_CONTINUE):
-
-{ .mmi
- add GR_Table_Base1 = 96, GR_Table_Base
- ldfs FR_Two_to_24 = [GR_Table_Base], 4
- nop.i 999
-}
-;;
-
-{ .mmi
- nop.m 999
-//
-// Load 2**24, load 2**63.
-//
- ldfs FR_Neg_Two_to_24 = [GR_Table_Base], 12
- mov r41 = ar.pfs ;;
-}
-
-{ .mfi
- ldfs FR_Two_to_63 = [GR_Table_Base1], 4
-//
-// Check for unnormals - unsupported operands. We do not want
-// to generate denormal exception
-// Check for NatVals, QNaNs, SNaNs, +/-Infs
-// Check for EM unsupporteds
-// Check for Zero
-//
- fclass.m.unc p6, p8 = FR_Input_X, 0x1E3
- mov r40 = gp ;;
-}
-
-{ .mfi
- nop.m 999
- fclass.nm.unc p8, p0 = FR_Input_X, 0x1FF
-// GR_Sin_or_Cos denotes
- mov r39 = b0
-}
-
-{ .mfb
- ldfs FR_Neg_Two_to_63 = [GR_Table_Base1], 12
- fclass.m.unc p10, p0 = FR_Input_X, 0x007
-(p6) br.cond.spnt L(SINCOS_SPECIAL) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-(p8) br.cond.spnt L(SINCOS_SPECIAL) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-//
-// Branch if +/- NaN, Inf.
-// Load -2**24, load -2**63.
-//
-(p10) br.cond.spnt L(SINCOS_ZERO) ;;
-}
-
-{ .mmb
- ldfe FR_Inv_pi_by_2 = [GR_Table_Base], 16
- ldfe FR_Inv_P_0 = [GR_Table_Base1], 16
- nop.b 999 ;;
-}
-
-{ .mmb
- nop.m 999
- ldfe FR_d_1 = [GR_Table_Base1], 16
- nop.b 999 ;;
-}
-//
-// Raise possible denormal operand flag with useful fcmp
-// Is x <= -2**63
-// Load Inv_P_0 for pre-reduction
-// Load Inv_pi_by_2
-//
-
-{ .mmb
- ldfe FR_P_0 = [GR_Table_Base], 16
- ldfe FR_d_2 = [GR_Table_Base1], 16
- nop.b 999 ;;
-}
-//
-// Load P_0
-// Load d_1
-// Is x >= 2**63
-// Is x <= -2**24?
-//
-
-{ .mmi
- ldfe FR_P_1 = [GR_Table_Base], 16 ;;
-//
-// Load P_1
-// Load d_2
-// Is x >= 2**24?
-//
- ldfe FR_P_2 = [GR_Table_Base], 16
- nop.i 999 ;;
-}
-
-{ .mmf
- nop.m 999
- ldfe FR_P_3 = [GR_Table_Base], 16
- fcmp.le.unc.s1 p7, p8 = FR_Input_X, FR_Neg_Two_to_24
-}
-
-{ .mfi
- nop.m 999
-//
-// Branch if +/- zero.
-// Decide about the paths to take:
-// If -2**24 < FR_Input_X < 2**24 - CASE 1 OR 2
-// OTHERWISE - CASE 3 OR 4
-//
- fcmp.le.unc.s0 p10, p11 = FR_Input_X, FR_Neg_Two_to_63
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p8) fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24
- nop.i 999
-}
-
-{ .mfi
- ldfe FR_Pi_by_4 = [GR_Table_Base1], 16
-(p11) fcmp.ge.s1 p10, p0 = FR_Input_X, FR_Two_to_63
- nop.i 999 ;;
-}
-
-{ .mmi
- ldfe FR_Neg_Pi_by_4 = [GR_Table_Base1], 16 ;;
- ldfs FR_Two_to_M3 = [GR_Table_Base1], 4
- nop.i 999 ;;
-}
-
-{ .mib
- ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1], 12
- nop.i 999
-//
-// Load P_2
-// Load P_3
-// Load pi_by_4
-// Load neg_pi_by_4
-// Load 2**(-3)
-// Load -2**(-3).
-//
-(p10) br.cond.spnt L(SINCOS_ARG_TOO_LARGE) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-//
-// Branch out if x >= 2**63. Use Payne-Hanek Reduction
-//
-(p7) br.cond.spnt L(SINCOS_LARGER_ARG) ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Branch if Arg <= -2**24 or Arg >= 2**24 and use pre-reduction.
-//
- fma.s1 FR_N_float = FR_Input_X, FR_Inv_pi_by_2, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
- fcmp.lt.unc.s1 p6, p7 = FR_Input_X, FR_Pi_by_4
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Select the case when |Arg| < pi/4
-// Else Select the case when |Arg| >= pi/4
-//
- fcvt.fx.s1 FR_N_fix = FR_N_float
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// N = Arg * 2/pi
-// Check if Arg < pi/4
-//
-(p6) fcmp.gt.s1 p6, p7 = FR_Input_X, FR_Neg_Pi_by_4
- nop.i 999 ;;
-}
-//
-// Case 2: Convert integer N_fix back to normalized floating-point value.
-// Case 1: p8 is only affected when p6 is set
-//
-
-{ .mfi
-(p7) ldfs FR_Two_to_M33 = [GR_Table_Base1], 4
-//
-// Grab the integer part of N and call it N_fix
-//
-(p6) fmerge.se FR_r = FR_Input_X, FR_Input_X
-// If |x| < pi/4, r = x and c = 0
-// lf |x| < pi/4, is x < 2**(-3).
-// r = Arg
-// c = 0
-(p6) mov GR_N_Inc = GR_Sin_or_Cos ;;
-}
-
{ .mmf
- nop.m 999
-(p7) ldfs FR_Neg_Two_to_M33 = [GR_Table_Base1], 4
-(p6) fmerge.se FR_c = f0, f0
-}
-
-{ .mfi
- nop.m 999
-(p6) fcmp.lt.unc.s1 p8, p9 = FR_Input_X, FR_Two_to_M3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8.
-// If |x| >= pi/4,
-// Create the right N for |x| < pi/4 and otherwise
-// Case 2: Place integer part of N in GP register
-//
-(p7) fcvt.xf FR_N_float = FR_N_fix
- nop.i 999 ;;
-}
-
-{ .mmf
- nop.m 999
-(p7) getf.sig GR_N_Inc = FR_N_fix
-(p8) fcmp.gt.s1 p8, p0 = FR_Input_X, FR_Neg_Two_to_M3 ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-//
-// Load 2**(-33), -2**(-33)
-//
-(p8) br.cond.spnt L(SINCOS_SMALL_R) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-(p6) br.cond.sptk L(SINCOS_NORMAL_R) ;;
-}
-//
-// if |x| < pi/4, branch based on |x| < 2**(-3) or otherwise.
-//
-//
-// In this branch, |x| >= pi/4.
-//
-
-{ .mfi
- ldfs FR_Neg_Two_to_M67 = [GR_Table_Base1], 8
-//
-// Load -2**(-67)
-//
- fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X
-//
-// w = N * P_2
-// s = -N * P_1 + Arg
-//
- add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos
-}
-
-{ .mfi
- nop.m 999
- fma.s1 FR_w = FR_N_float, FR_P_2, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Adjust N_fix by N_inc to determine whether sine or
-// cosine is being calculated
-//
- fcmp.lt.unc.s1 p7, p6 = FR_s, FR_Two_to_M33
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-// Remember x >= pi/4.
-// Is s <= -2**(-33) or s >= 2**(-33) (p6)
-// or -2**(-33) < s < 2**(-33) (p7)
-(p6) fms.s1 FR_r = FR_s, f1, FR_w
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p6) fms.s1 FR_c = FR_s, f1, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// For big s: r = s - w: No futher reduction is necessary
-// For small s: w = N * P_3 (change sign) More reduction
-//
-(p6) fcmp.lt.unc.s1 p8, p9 = FR_r, FR_Two_to_M3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p8) fcmp.gt.s1 p8, p9 = FR_r, FR_Neg_Two_to_M3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fms.s1 FR_r = FR_s, f1, FR_U_1
- nop.i 999
-}
-
-{ .mfb
- nop.m 999
-//
-// For big s: Is |r| < 2**(-3)?
-// For big s: c = S - r
-// For small s: U_1 = N * P_2 + w
-//
-// If p8 is set, prepare to branch to Small_R.
-// If p9 is set, prepare to branch to Normal_R.
-// For big s, r is complete here.
-//
-(p6) fms.s1 FR_c = FR_c, f1, FR_w
-//
-// For big s: c = c + w (w has not been negated.)
-// For small s: r = S - U_1
-//
-(p8) br.cond.spnt L(SINCOS_SMALL_R) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-(p9) br.cond.sptk L(SINCOS_NORMAL_R) ;;
-}
-
-{ .mfi
-(p7) add GR_Table_Base1 = 224, GR_Table_Base1
-//
-// Branch to SINCOS_SMALL_R or SINCOS_NORMAL_R
-//
-(p7) fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
-//
-// c = S - U_1
-// r = S_1 * r
-//
-//
-(p7) extr.u GR_i_1 = GR_N_Inc, 0, 1
-}
-
-{ .mmi
- nop.m 999 ;;
-//
-// Get [i_0,i_1] - two lsb of N_fix_gr.
-// Do dummy fmpy so inexact is always set.
-//
-(p7) cmp.eq.unc p9, p10 = 0x0, GR_i_1
-(p7) extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
-}
-//
-// For small s: U_2 = N * P_2 - U_1
-// S_1 stored constant - grab the one stored with the
-// coefficients.
-//
-
-{ .mfi
-(p7) ldfe FR_S_1 = [GR_Table_Base1], 16
-//
-// Check if i_1 and i_0 != 0
-//
-(p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67
-(p7) cmp.eq.unc p11, p12 = 0x0, GR_i_0 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fms.s1 FR_s = FR_s, f1, FR_r
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// S = S - r
-// U_2 = U_2 + w
-// load S_1
-//
-(p7) fma.s1 FR_rsq = FR_r, FR_r, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fma.s1 FR_U_2 = FR_U_2, f1, FR_w
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p7) fmerge.se FR_Input_X = FR_r, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_Input_X = f0, f1, f1
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// FR_rsq = r * r
-// Save r as the result.
-//
-(p7) fms.s1 FR_c = FR_s, f1, FR_U_1
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if ( i_1 ==0) poly = c + S_1*r*r*r
-// else Result = 1
-//
-(p12) fnma.s1 FR_Input_X = FR_Input_X, f1, f0
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p7) fma.s1 FR_r = FR_S_1, FR_r, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fma.d.s0 FR_S_1 = FR_S_1, FR_S_1, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// If i_1 != 0, poly = 2**(-67)
-//
-(p7) fms.s1 FR_c = FR_c, f1, FR_U_2
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// c = c - U_2
-//
-(p9) fma.s1 FR_poly = FR_r, FR_rsq, FR_c
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// i_0 != 0, so Result = -Result
-//
-(p11) fma.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
- nop.i 999 ;;
-}
-
-{ .mfb
- nop.m 999
-(p12) fms.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
-//
-// if (i_0 == 0), Result = Result + poly
-// else Result = Result - poly
-//
- br.ret.sptk b0 ;;
-}
-L(SINCOS_LARGER_ARG):
+ getf.sig sincos_GR_n = sincos_W_2TO61_RSH
+ ldfpd sincos_P2,sincos_Q2 = [sincos_AD_1],16
+ fms.s1 sincos_NFLOAT = sincos_W_2TO61_RSH,sincos_2TOM61,sincos_RSHF
+};;
+// sincos_r = -sincos_Nfloat * sincos_Pi_by_16_1 + x
{ .mfi
- nop.m 999
- fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0
- nop.i 999
-}
-;;
-
-// This path for argument > 2*24
-// Adjust table_ptr1 to beginning of table.
-//
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
+ ldfpd sincos_P1,sincos_Q1 = [sincos_AD_1],16
+ fnma.s1 sincos_r = sincos_NFLOAT, sincos_Pi_by_16_1, sincos_NORM_f8
+ nop.i 999
+};;
+// Add 2^(k-1) (which is in sincos_r_sincos) to N
{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
+ add sincos_GR_n = sincos_GR_n, sincos_r_sincos
;;
+// Get M (least k+1 bits of N)
+ and sincos_GR_m = 0x1f,sincos_GR_n
+ nop.i 999
+};;
-
-//
-// Point to 2*-14
-// N_0 = Arg * Inv_P_0
-//
-
-{ .mmi
- add GR_Table_Base = 688, GR_Table_Base ;;
- ldfs FR_Two_to_M14 = [GR_Table_Base], 4
- nop.i 999 ;;
-}
-
-{ .mfi
- ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0
- nop.f 999
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Load values 2**(-14) and -2**(-14)
-//
- fcvt.fx.s1 FR_N_0_fix = FR_N_0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// N_0_fix = integer part of N_0
-//
- fcvt.xf FR_N_0 = FR_N_0_fix
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Make N_0 the integer part
-//
- fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
- fma.s1 FR_w = FR_N_0, FR_d_1, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Arg' = -N_0 * P_0 + Arg
-// w = N_0 * d_1
-//
- fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// N = A' * 2/pi
-//
- fcvt.fx.s1 FR_N_fix = FR_N_float
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// N_fix is the integer part
-//
- fcvt.xf FR_N_float = FR_N_fix
- nop.i 999 ;;
-}
-
+// sincos_r = sincos_r -sincos_Nfloat * sincos_Pi_by_16_2
{ .mfi
- getf.sig GR_N_Inc = FR_N_fix
- nop.f 999
- nop.i 999 ;;
-}
-
-{ .mii
- nop.m 999
- nop.i 999 ;;
- add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;;
-}
+ nop.m 999
+ fnma.s1 sincos_r = sincos_NFLOAT, sincos_Pi_by_16_2, sincos_r
+ shl sincos_GR_32m = sincos_GR_m,5
+};;
+// Add 32*M to address of sin_cos_beta table
+// For sin denorm. - set uflow
{ .mfi
- nop.m 999
-//
-// N is the integer part of the reduced-reduced argument.
-// Put the integer in a GP register
-//
- fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
- nop.i 999
-}
+ add sincos_AD_2 = sincos_GR_32m, sincos_AD_1
+(p8) fclass.m.unc p10,p0 = f8,0x0b
+ nop.i 999
+};;
+// Load Sin and Cos table value using obtained index m (sincosf_AD_2)
{ .mfi
- nop.m 999
- fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
- nop.i 999 ;;
-}
+ ldfe sincos_Sm = [sincos_AD_2],16
+ nop.f 999
+ nop.i 999
+};;
+// get rsq = r*r
{ .mfi
- nop.m 999
-//
-// s = -N*P_1 + Arg'
-// w = -N*P_2 + w
-// N_fix_gr = N_fix_gr + N_inc
-//
- fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14
- nop.i 999 ;;
+ ldfe sincos_Cm = [sincos_AD_2]
+ fma.s1 sincos_rsq = sincos_r, sincos_r, f0 // r^2 = r*r
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-(p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14
- nop.i 999 ;;
-}
+ nop.m 999
+ fmpy.s0 fp_tmp = fp_tmp,fp_tmp // forces inexact flag
+ nop.i 999
+};;
+// sincos_r_exact = sincos_r -sincos_Nfloat * sincos_Pi_by_16_3
{ .mfi
- nop.m 999
-//
-// For |s| > 2**(-14) r = S + w (r complete)
-// Else U_hi = N_0 * d_1
-//
-(p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0
- nop.i 999
-}
+ nop.m 999
+ fnma.s1 sincos_r_exact = sincos_NFLOAT, sincos_Pi_by_16_3, sincos_r
+ nop.i 999
+};;
+// Polynomials calculation
+// P_1 = P4*r^2 + P3
+// Q_2 = Q4*r^2 + Q3
{ .mfi
- nop.m 999
-(p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
- nop.i 999 ;;
+ nop.m 999
+ fma.s1 sincos_P_temp1 = sincos_rsq, sincos_P4, sincos_P3
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-//
-// Either S <= -2**(-14) or S >= 2**(-14)
-// or -2**(-14) < s < 2**(-14)
-//
-(p8) fma.s1 FR_r = FR_s, f1, FR_w
- nop.i 999
-}
+ nop.m 999
+ fma.s1 sincos_Q_temp1 = sincos_rsq, sincos_Q4, sincos_Q3
+ nop.i 999
+};;
+// get rcube = r^3 and S[m]*r^2
{ .mfi
- nop.m 999
-(p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0
- nop.i 999 ;;
+ nop.m 999
+ fmpy.s1 sincos_srsq = sincos_Sm,sincos_rsq
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-//
-// We need abs of both U_hi and V_hi - don't
-// worry about switched sign of V_hi.
-//
-(p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi
- nop.i 999
-}
+ nop.m 999
+ fmpy.s1 sincos_rcub = sincos_r_exact, sincos_rsq
+ nop.i 999
+};;
+// Polynomials calculation
+// Q_2 = Q_1*r^2 + Q2
+// P_1 = P_1*r^2 + P2
{ .mfi
- nop.m 999
-//
-// Big s: finish up c = (S - r) + w (c complete)
-// Case 4: A = U_hi + V_hi
-// Note: Worry about switched sign of V_hi, so subtract instead of add.
-//
-(p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi
- nop.i 999 ;;
+ nop.m 999
+ fma.s1 sincos_Q_temp2 = sincos_rsq, sincos_Q_temp1, sincos_Q2
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-(p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
- nop.i 999 ;;
-}
+ nop.m 999
+ fma.s1 sincos_P_temp2 = sincos_rsq, sincos_P_temp1, sincos_P2
+ nop.i 999
+};;
+// Polynomials calculation
+// Q = Q_2*r^2 + Q1
+// P = P_2*r^2 + P1
{ .mfi
- nop.m 999
-(p9) fmerge.s FR_V_hiabs = f0, FR_V_hi
- nop.i 999
+ nop.m 999
+ fma.s1 sincos_Q = sincos_rsq, sincos_Q_temp2, sincos_Q1
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-// For big s: c = S - r
-// For small s do more work: U_lo = N_0 * d_1 - U_hi
-//
-(p9) fmerge.s FR_U_hiabs = f0, FR_U_hi
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// For big s: Is |r| < 2**(-3)
-// For big s: if p12 set, prepare to branch to Small_R.
-// For big s: If p13 set, prepare to branch to Normal_R.
-//
-(p8) fms.s1 FR_c = FR_s, f1, FR_r
- nop.i 999
-}
+ nop.m 999
+ fma.s1 sincos_P = sincos_rsq, sincos_P_temp2, sincos_P1
+ nop.i 999
+};;
+// Get final P and Q
+// Q = Q*S[m]*r^2 + S[m]
+// P = P*r^3 + r
{ .mfi
- nop.m 999
-//
-// For small S: V_hi = N * P_2
-// w = N * P_3
-// Note the product does not include the (-) as in the writeup
-// so (-) missing for V_hi and w.
-//
-(p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3
- nop.i 999 ;;
+ nop.m 999
+ fma.s1 sincos_Q = sincos_srsq,sincos_Q, sincos_Sm
+ nop.i 999
}
-
{ .mfi
- nop.m 999
-(p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
- nop.i 999 ;;
-}
+ nop.m 999
+ fma.s1 sincos_P = sincos_rcub,sincos_P, sincos_r_exact
+ nop.i 999
+};;
+// If sin(denormal), force underflow to be set
{ .mfi
- nop.m 999
-(p8) fma.s1 FR_c = FR_c, f1, FR_w
- nop.i 999
-}
+ nop.m 999
+(p10) fmpy.d.s0 fp_tmp = sincos_NORM_f8,sincos_NORM_f8
+ nop.i 999
+};;
+// Final calculation
+// result = C[m]*P + Q
{ .mfb
- nop.m 999
-(p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w
-(p12) br.cond.spnt L(SINCOS_SMALL_R) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-(p13) br.cond.sptk L(SINCOS_NORMAL_R) ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true.
-// The remaining stuff is for Case 4.
-// Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup)
-// Note: the (-) is still missing for V_lo.
-// Small s: w = w + N_0 * d_2
-// Note: the (-) is now incorporated in w.
-//
-(p9) fcmp.ge.unc.s1 p10, p11 = FR_U_hiabs, FR_V_hiabs
- extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// C_hi = S + A
-//
-(p9) fma.s1 FR_t = FR_U_lo, f1, FR_V_lo
- extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// t = U_lo + V_lo
-//
-//
-(p10) fms.s1 FR_a = FR_U_hi, f1, FR_A
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p11) fma.s1 FR_a = FR_V_hi, f1, FR_A
- nop.i 999
-}
-;;
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
-
-{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
-;;
-
-
-{ .mfi
- add GR_Table_Base = 528, GR_Table_Base
-//
-// Is U_hiabs >= V_hiabs?
-//
-(p9) fma.s1 FR_C_hi = FR_s, f1, FR_A
- nop.i 999 ;;
-}
-
-{ .mmi
- ldfe FR_C_1 = [GR_Table_Base], 16 ;;
- ldfe FR_C_2 = [GR_Table_Base], 64
- nop.i 999 ;;
-}
-
-{ .mmf
- nop.m 999
-//
-// c = c + C_lo finished.
-// Load C_2
-//
- ldfe FR_S_1 = [GR_Table_Base], 16
-//
-// C_lo = S - C_hi
-//
- fma.s1 FR_t = FR_t, f1, FR_w ;;
-}
-//
-// r and c have been computed.
-// Make sure ftz mode is set - should be automatic when using wre
-// |r| < 2**(-3)
-// Get [i_0,i_1] - two lsb of N_fix.
-// Load S_1
-//
-
-{ .mfi
- ldfe FR_S_2 = [GR_Table_Base], 64
-//
-// t = t + w
-//
-(p10) fms.s1 FR_a = FR_a, f1, FR_V_hi
- cmp.eq.unc p9, p10 = 0x0, GR_i_0
-}
-
-{ .mfi
- nop.m 999
-//
-// For larger u than v: a = U_hi - A
-// Else a = V_hi - A (do an add to account for missing (-) on V_hi
-//
- fms.s1 FR_C_lo = FR_s, f1, FR_C_hi
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p11) fms.s1 FR_a = FR_U_hi, f1, FR_a
- cmp.eq.unc p11, p12 = 0x0, GR_i_1
-}
-
-{ .mfi
- nop.m 999
-//
-// If u > v: a = (U_hi - A) + V_hi
-// Else a = (V_hi - A) + U_hi
-// In each case account for negative missing from V_hi.
-//
- fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// C_lo = (S - C_hi) + A
-//
- fma.s1 FR_t = FR_t, f1, FR_a
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// t = t + a
-//
- fma.s1 FR_C_lo = FR_C_lo, f1, FR_t
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// C_lo = C_lo + t
-// Adjust Table_Base to beginning of table
-//
- fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Load S_2
-//
- fma.s1 FR_rsq = FR_r, FR_r, f0
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// Table_Base points to C_1
-// r = C_hi + C_lo
-//
- fms.s1 FR_c = FR_C_hi, f1, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if i_1 ==0: poly = S_2 * FR_rsq + S_1
-// else poly = C_2 * FR_rsq + C_1
-//
-(p11) fma.s1 FR_Input_X = f0, f1, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fma.s1 FR_Input_X = f0, f1, f1
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Compute r_cube = FR_rsq * r
-//
-(p11) fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// Compute FR_rsq = r * r
-// Is i_1 == 0 ?
-//
- fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// c = C_hi - r
-// Load C_1
-//
- fma.s1 FR_c = FR_c, f1, FR_C_lo
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// if i_1 ==0: poly = r_cube * poly + c
-// else poly = FR_rsq * poly
-//
-(p10) fms.s1 FR_Input_X = f0, f1, FR_Input_X
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if i_1 ==0: Result = r
-// else Result = 1.0
-//
-(p11) fma.s1 FR_poly = FR_r_cubed, FR_poly, FR_c
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fma.s1 FR_poly = FR_rsq, FR_poly, f0
- nop.i 999 ;;
-}
+ nop.m 999
+ fma.d.s0 f8 = sincos_Cm, sincos_P, sincos_Q
+ br.ret.sptk b0 // Exit for common path
+};;
+////////// x = 0/Inf/NaN path //////////////////
+_SINCOS_SPECIAL_ARGS:
+.pred.rel "mutex",p8,p9
+// sin(+/-0) = +/-0
+// sin(Inf) = NaN
+// sin(NaN) = NaN
{ .mfi
- nop.m 999
-//
-// if i_0 !=0: Result = -Result
-//
-(p9) fma.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
- nop.i 999 ;;
+ nop.m 999
+(p8) fma.d.s0 f8 = f8, f0, f0 // sin(+/-0,NaN,Inf)
+ nop.i 999
}
-
+// cos(+/-0) = 1.0
+// cos(Inf) = NaN
+// cos(NaN) = NaN
{ .mfb
- nop.m 999
-(p10) fms.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
-//
-// if i_0 == 0: Result = Result + poly
-// else Result = Result - poly
-//
- br.ret.sptk b0 ;;
-}
-L(SINCOS_SMALL_R):
-
-{ .mii
- nop.m 999
- extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
-//
-//
-// Compare both i_1 and i_0 with 0.
-// if i_1 == 0, set p9.
-// if i_0 == 0, set p11.
-//
- cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
-}
-
-{ .mfi
- nop.m 999
- fma.s1 FR_rsq = FR_r, FR_r, f0
- extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Z = Z * FR_rsq
-//
-(p10) fnma.s1 FR_c = FR_c, FR_r, f0
- cmp.eq.unc p11, p12 = 0x0, GR_i_0
-}
-;;
-
-// ******************************************************************
-// ******************************************************************
-// ******************************************************************
-// r and c have been computed.
-// We know whether this is the sine or cosine routine.
-// Make sure ftz mode is set - should be automatic when using wre
-// |r| < 2**(-3)
-//
-// Set table_ptr1 to beginning of constant table.
-// Get [i_0,i_1] - two lsb of N_fix_gr.
-//
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
-
-{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
-;;
-
-
-//
-// Set table_ptr1 to point to S_5.
-// Set table_ptr1 to point to C_5.
-// Compute FR_rsq = r * r
-//
-
-{ .mfi
-(p9) add GR_Table_Base = 672, GR_Table_Base
-(p10) fmerge.s FR_r = f1, f1
-(p10) add GR_Table_Base = 592, GR_Table_Base ;;
-}
-//
-// Set table_ptr1 to point to S_5.
-// Set table_ptr1 to point to C_5.
-//
-
-{ .mmi
-(p9) ldfe FR_S_5 = [GR_Table_Base], -16 ;;
-//
-// if (i_1 == 0) load S_5
-// if (i_1 != 0) load C_5
-//
-(p9) ldfe FR_S_4 = [GR_Table_Base], -16
- nop.i 999 ;;
-}
-
-{ .mmf
-(p10) ldfe FR_C_5 = [GR_Table_Base], -16
-//
-// Z = FR_rsq * FR_rsq
-//
-(p9) ldfe FR_S_3 = [GR_Table_Base], -16
-//
-// Compute FR_rsq = r * r
-// if (i_1 == 0) load S_4
-// if (i_1 != 0) load C_4
-//
- fma.s1 FR_Z = FR_rsq, FR_rsq, f0 ;;
-}
-//
-// if (i_1 == 0) load S_3
-// if (i_1 != 0) load C_3
-//
-
-{ .mmi
-(p9) ldfe FR_S_2 = [GR_Table_Base], -16 ;;
-//
-// if (i_1 == 0) load S_2
-// if (i_1 != 0) load C_2
-//
-(p9) ldfe FR_S_1 = [GR_Table_Base], -16
- nop.i 999
-}
-
-{ .mmi
-(p10) ldfe FR_C_4 = [GR_Table_Base], -16 ;;
-(p10) ldfe FR_C_3 = [GR_Table_Base], -16
- nop.i 999 ;;
-}
-
-{ .mmi
-(p10) ldfe FR_C_2 = [GR_Table_Base], -16 ;;
-(p10) ldfe FR_C_1 = [GR_Table_Base], -16
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 != 0):
-// poly_lo = FR_rsq * C_5 + C_4
-// poly_hi = FR_rsq * C_2 + C_1
-//
-(p9) fma.s1 FR_Z = FR_Z, FR_r, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 == 0) load S_1
-// if (i_1 != 0) load C_1
-//
-(p9) fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// c = -c * r
-// dummy fmpy's to flag inexact.
-//
-(p9) fma.d.s0 FR_S_4 = FR_S_4, FR_S_4, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// poly_lo = FR_rsq * poly_lo + C_3
-// poly_hi = FR_rsq * poly_hi
-//
- fma.s1 FR_Z = FR_Z, FR_rsq, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 == 0):
-// poly_lo = FR_rsq * S_5 + S_4
-// poly_hi = FR_rsq * S_2 + S_1
-//
-(p10) fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 == 0):
-// Z = Z * r for only one of the small r cases - not there
-// in original implementation notes.
-//
-(p9) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.d.s0 FR_C_1 = FR_C_1, FR_C_1, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// poly_lo = FR_rsq * poly_lo + S_3
-// poly_hi = FR_rsq * poly_hi
-//
-(p10) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 == 0): dummy fmpy's to flag inexact
-// r = 1
-//
-(p9) fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// poly_hi = r * poly_hi
-//
- fma.s1 FR_poly = FR_Z, FR_poly_lo, FR_c
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fms.s1 FR_r = f0, f1, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// poly_hi = Z * poly_lo + c
-// if i_0 == 1: r = -r
-//
- fma.s1 FR_poly = FR_poly, f1, FR_poly_hi
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fms.d.s0 FR_Input_X = FR_r, f1, FR_poly
- nop.i 999
-}
+ nop.m 999
+(p9) fma.d.s0 f8 = f8, f0, f1 // cos(+/-0,NaN,Inf)
+ br.ret.sptk b0 // Exit for x = 0/Inf/NaN path
+};;
+_SINCOS_UNORM:
+// Here if x=unorm
{ .mfb
- nop.m 999
-//
-// poly = poly + poly_hi
-//
-(p11) fma.d.s0 FR_Input_X = FR_r, f1, FR_poly
-//
-// if (i_0 == 0) Result = r + poly
-// if (i_0 != 0) Result = r - poly
-//
- br.ret.sptk b0 ;;
-}
-L(SINCOS_NORMAL_R):
-
-{ .mii
- nop.m 999
- extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
-//
-// Set table_ptr1 and table_ptr2 to base address of
-// constant table.
- cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
-}
+ getf.exp sincos_r_signexp = sincos_NORM_f8 // Get signexp of x
+ fcmp.eq.s0 p11,p0 = f8, f0 // Dummy op to set denorm flag
+ br.cond.sptk _SINCOS_COMMON2 // Return to main path
+};;
-{ .mfi
- nop.m 999
- fma.s1 FR_rsq = FR_r, FR_r, f0
- extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
-}
+GLOBAL_IEEE754_END(cos)
+//////////// x >= 2^27 - large arguments routine call ////////////
+LOCAL_LIBM_ENTRY(__libm_callout_sincos)
+_SINCOS_LARGE_ARGS:
+.prologue
{ .mfi
- nop.m 999
- frcpa.s1 FR_r_hi, p6 = f1, FR_r
- cmp.eq.unc p11, p12 = 0x0, GR_i_0
+ mov GR_SAVE_r_sincos = sincos_r_sincos // Save sin or cos
+ nop.f 999
+.save ar.pfs,GR_SAVE_PFS
+ mov GR_SAVE_PFS = ar.pfs
}
;;
-// ******************************************************************
-// ******************************************************************
-// ******************************************************************
-//
-// r and c have been computed.
-// We known whether this is the sine or cosine routine.
-// Make sure ftz mode is set - should be automatic when using wre
-// Get [i_0,i_1] - two lsb of N_fix_gr alone.
-//
-
-{ .mmi
- nop.m 999
- addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
- nop.i 999
-}
-;;
-
-{ .mmi
- ld8 GR_Table_Base = [GR_Table_Base]
- nop.m 999
- nop.i 999
-}
-;;
-
-
-{ .mfi
-(p10) add GR_Table_Base = 384, GR_Table_Base
-(p12) fms.s1 FR_Input_X = f0, f1, f1
-(p9) add GR_Table_Base = 224, GR_Table_Base ;;
-}
-
-{ .mmf
- nop.m 999
-(p10) ldfe FR_QQ_8 = [GR_Table_Base], 16
-//
-// if (i_1==0) poly = poly * FR_rsq + PP_1_lo
-// else poly = FR_rsq * poly
-//
-(p11) fma.s1 FR_Input_X = f0, f1, f1 ;;
-}
-
-{ .mmf
-(p10) ldfe FR_QQ_7 = [GR_Table_Base], 16
-//
-// Adjust table pointers based on i_0
-// Compute rsq = r * r
-//
-(p9) ldfe FR_PP_8 = [GR_Table_Base], 16
- fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 ;;
-}
-
-{ .mmf
-(p9) ldfe FR_PP_7 = [GR_Table_Base], 16
-(p10) ldfe FR_QQ_6 = [GR_Table_Base], 16
-//
-// Load PP_8 and QQ_8; PP_7 and QQ_7
-//
- frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi ;;
-}
-//
-// if (i_1==0) poly = PP_7 + FR_rsq * PP_8.
-// else poly = QQ_7 + FR_rsq * QQ_8.
-//
-
-{ .mmb
-(p9) ldfe FR_PP_6 = [GR_Table_Base], 16
-(p10) ldfe FR_QQ_5 = [GR_Table_Base], 16
- nop.b 999 ;;
-}
-
-{ .mmb
-(p9) ldfe FR_PP_5 = [GR_Table_Base], 16
-(p10) ldfe FR_S_1 = [GR_Table_Base], 16
- nop.b 999 ;;
-}
-
-{ .mmb
-(p10) ldfe FR_QQ_1 = [GR_Table_Base], 16
-(p9) ldfe FR_C_1 = [GR_Table_Base], 16
- nop.b 999 ;;
-}
-
-{ .mmi
-(p10) ldfe FR_QQ_4 = [GR_Table_Base], 16 ;;
-(p9) ldfe FR_PP_1 = [GR_Table_Base], 16
- nop.i 999 ;;
-}
-
-{ .mmf
-(p10) ldfe FR_QQ_3 = [GR_Table_Base], 16
-//
-// if (i_1=0) corr = corr + c*c
-// else corr = corr * c
-//
-(p9) ldfe FR_PP_4 = [GR_Table_Base], 16
-(p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 ;;
-}
-//
-// if (i_1=0) poly = rsq * poly + PP_5
-// else poly = rsq * poly + QQ_5
-// Load PP_4 or QQ_4
-//
-
-{ .mmf
-(p9) ldfe FR_PP_3 = [GR_Table_Base], 16
-(p10) ldfe FR_QQ_2 = [GR_Table_Base], 16
-//
-// r_hi = frcpa(frcpa(r)).
-// r_cube = r * FR_rsq.
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 ;;
-}
-//
-// Do dummy multiplies so inexact is always set.
-//
-
-{ .mfi
-(p9) ldfe FR_PP_2 = [GR_Table_Base], 16
-//
-// r_lo = r - r_hi
-//
-(p9) fma.s1 FR_U_lo = FR_r_hi, FR_r_hi, f0
- nop.i 999 ;;
-}
-
-{ .mmf
- nop.m 999
-(p9) ldfe FR_PP_1_lo = [GR_Table_Base], 16
-(p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1=0) U_lo = r_hi * r_hi
-// else U_lo = r_hi + r
-//
-(p9) fma.s1 FR_corr = FR_C_1, FR_rsq, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1=0) corr = C_1 * rsq
-// else corr = S_1 * r_cubed + r
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1=0) U_hi = r_hi + U_hi
-// else U_hi = QQ_1 * U_hi + 1
-//
-(p9) fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_U_lo
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-//
-// U_hi = r_hi * r_hi
-//
- fms.s1 FR_r_lo = FR_r, f1, FR_r_hi
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Load PP_1, PP_6, PP_5, and C_1
-// Load QQ_1, QQ_6, QQ_5, and S_1
-//
- fma.s1 FR_U_hi = FR_r_hi, FR_r_hi, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fnma.s1 FR_corr = FR_corr, FR_c, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1=0) U_lo = r * r_hi + U_lo
-// else U_lo = r_lo * U_lo
-//
-(p9) fma.s1 FR_corr = FR_corr, FR_c, FR_c
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5
- nop.i 999
-}
-
{ .mfi
- nop.m 999
-//
-// if (i_1 =0) U_hi = r + U_hi
-// if (i_1 =0) U_lo = r_lo * U_lo
-//
-//
-(p9) fma.d.s0 FR_PP_5 = FR_PP_5, FR_PP_4, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1=0) poly = poly * rsq + PP_6
-// else poly = poly * rsq + QQ_6
-//
-(p9) fma.s1 FR_U_hi = FR_r_hi, FR_U_hi, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_U_hi = FR_QQ_1, FR_U_hi, f1
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.d.s0 FR_QQ_5 = FR_QQ_5, FR_QQ_5, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1!=0) U_hi = PP_1 * U_hi
-// if (i_1!=0) U_lo = r * r + U_lo
-// Load PP_3 or QQ_3
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
- nop.i 999
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// Load PP_2, QQ_2
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3
- nop.i 999 ;;
+ mov GR_SAVE_GP = gp
+ nop.f 999
+.save b0, GR_SAVE_B0
+ mov GR_SAVE_B0 = b0
}
-{ .mfi
- nop.m 999
-//
-// if (i_1==0) poly = FR_rsq * poly + PP_3
-// else poly = FR_rsq * poly + QQ_3
-// Load PP_1_lo
-//
-(p9) fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1 =0) poly = poly * rsq + pp_r4
-// else poly = poly * rsq + qq_r4
-//
-(p9) fma.s1 FR_U_hi = FR_r, f1, FR_U_hi
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1==0) U_lo = PP_1_hi * U_lo
-// else U_lo = QQ_1 * U_lo
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_0==0) Result = 1
-// else Result = -1
-//
- fma.s1 FR_V = FR_U_lo, f1, FR_corr
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1==0) poly = FR_rsq * poly + PP_2
-// else poly = FR_rsq * poly + QQ_2
-//
-(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// V = U_lo + corr
-//
-(p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, f0
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-//
-// if (i_1==0) poly = r_cube * poly
-// else poly = FR_rsq * poly
-//
- fma.s1 FR_V = FR_poly, f1, FR_V
- nop.i 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p12) fms.d.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
- nop.i 999
-}
-
-{ .mfb
- nop.m 999
-//
-// V = V + poly
-//
-(p11) fma.d.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
-//
-// if (i_0==0) Result = Result * U_hi + V
-// else Result = Result * U_hi - V
-//
- br.ret.sptk b0 ;;
-}
-
-//
-// If cosine, FR_Input_X = 1
-// If sine, FR_Input_X = +/-Zero (Input FR_Input_X)
-// Results are exact, no exceptions
-//
-L(SINCOS_ZERO):
-
-{ .mmb
- cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
- nop.m 999
- nop.b 999 ;;
-}
-
-{ .mfi
- nop.m 999
-(p7) fmerge.s FR_Input_X = FR_Input_X, FR_Input_X
- nop.i 999
-}
-
-{ .mfb
- nop.m 999
-(p6) fmerge.s FR_Input_X = f1, f1
- br.ret.sptk b0 ;;
-}
-
-L(SINCOS_SPECIAL):
-
-//
-// Path for Arg = +/- QNaN, SNaN, Inf
-// Invalid can be raised. SNaNs
-// become QNaNs
-//
-
-{ .mfb
- nop.m 999
- fmpy.d.s0 FR_Input_X = FR_Input_X, f0
- br.ret.sptk b0 ;;
-}
-.endp __libm_cos_double_dbx#
-ASM_SIZE_DIRECTIVE(__libm_cos_double_dbx#)
-
-
-
-//
-// Call int pi_by_2_reduce(double* x, double *y)
-// for |arguments| >= 2**63
-// Address to save r and c as double
-//
-//
-// psp sp+64
-// sp+48 -> f0 c
-// r45 sp+32 -> f0 r
-// r44 -> sp+16 -> InputX
-// sp sp -> scratch provided to callee
-
+.body
+{ .mbb
+ setf.sig sincos_save_tmp = sincos_GR_all_ones// inexact set
+ nop.b 999
+(p8) br.call.sptk.many b0 = __libm_sin_large# // sin(large_X)
+};;
-.proc __libm_callout_2
-__libm_callout_2:
-L(SINCOS_ARG_TOO_LARGE):
+{ .mbb
+ cmp.ne p9,p0 = GR_SAVE_r_sincos, r0 // set p9 if cos
+ nop.b 999
+(p9) br.call.sptk.many b0 = __libm_cos_large# // cos(large_X)
+};;
-.prologue
{ .mfi
- add r45=-32,sp // Parameter: r address
- nop.f 0
-.save ar.pfs,GR_SAVE_PFS
- mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
+ mov gp = GR_SAVE_GP
+ fma.d.s0 f8 = f8, f1, f0 // Round result to double
+ mov b0 = GR_SAVE_B0
}
+// Force inexact set
{ .mfi
-.fframe 64
- add sp=-64,sp // Create new stack
- nop.f 0
- mov GR_SAVE_GP=gp // Save gp
-};;
-{ .mmi
- stfe [r45] = f0,16 // Clear Parameter r on stack
- add r44 = 16,sp // Parameter x address
-.save b0, GR_SAVE_B0
- mov GR_SAVE_B0=b0 // Save b0
-};;
-.body
-{ .mib
- stfe [r45] = f0,-16 // Clear Parameter c on stack
- nop.i 0
- nop.b 0
-}
-{ .mib
- stfe [r44] = FR_Input_X // Store Parameter x on stack
- nop.i 0
- br.call.sptk b0=__libm_pi_by_2_reduce# ;;
+ nop.m 999
+ fmpy.s0 sincos_save_tmp = sincos_save_tmp, sincos_save_tmp
+ nop.i 999
};;
-
-{ .mii
- ldfe FR_Input_X =[r44],16
-//
-// Get r and c off stack
-//
- adds GR_Table_Base1 = -16, GR_Table_Base1
-//
-// Get r and c off stack
-//
- add GR_N_Inc = GR_Sin_or_Cos,r8 ;;
-}
-{ .mmb
- ldfe FR_r =[r45],16
-//
-// Get X off the stack
-// Readjust Table ptr
-//
- ldfs FR_Two_to_M3 = [GR_Table_Base1],4
- nop.b 999 ;;
-}
-{ .mmb
- ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1],0
- ldfe FR_c =[r45]
- nop.b 999 ;;
-}
-
-{ .mfi
-.restore sp
- add sp = 64,sp // Restore stack pointer
- fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3
- mov b0 = GR_SAVE_B0 // Restore return address
-};;
{ .mib
- mov gp = GR_SAVE_GP // Restore gp
- mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
- nop.b 0
+ nop.m 999
+ mov ar.pfs = GR_SAVE_PFS
+ br.ret.sptk b0 // Exit for large arguments routine call
};;
+LOCAL_LIBM_END(__libm_callout_sincos)
-{ .mfi
- nop.m 999
-(p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3
- nop.i 999 ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
-(p6) br.cond.spnt L(SINCOS_SMALL_R) ;;
-}
-
-{ .mib
- nop.m 999
- nop.i 999
- br.cond.sptk L(SINCOS_NORMAL_R) ;;
-}
-
-.endp __libm_callout_2
-ASM_SIZE_DIRECTIVE(__libm_callout_2)
-
-.type __libm_pi_by_2_reduce#,@function
-.global __libm_pi_by_2_reduce#
-
+.type __libm_sin_large#,@function
+.global __libm_sin_large#
+.type __libm_cos_large#,@function
+.global __libm_cos_large#
-.type __libm_sin_double_dbx#,@function
-.global __libm_sin_double_dbx#
-.type __libm_cos_double_dbx#,@function
-.global __libm_cos_double_dbx#