about summary refs log tree commit diff
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, 3115 insertions, 382 deletions
diff --git a/sysdeps/ia64/fpu/s_cos.S b/sysdeps/ia64/fpu/s_cos.S
index fc121fce19..6540aec724 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 - 2005, Intel Corporation
+// Copyright (C) 2000, 2001, Intel Corporation
 // All rights reserved.
 //
-// Contributed 2000 by the Intel Numerics Group, Intel Corporation
+// 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.
 //
 // 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,25 +35,17 @@
 //
 // Intel Corporation is the author of this code, and requests that all
 // problem reports or change requests be submitted to it directly at
-// http://www.intel.com/software/products/opensource/libraries/num.htm.
+// http://developer.intel.com/opensource.
 //
 // History
 //==============================================================
-// 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
+// 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
 // 10/18/00 Changed one table entry to ensure symmetry
-// 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
+// 1/03/01  Improved speed, fixed flag settings for small arguments.
 
 // API
 //==============================================================
@@ -71,13 +63,9 @@
 //    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)) - 
-//                        nfloat * LOWEST(pi/2^k) for increased accuracy.
+//    Do this as (x -  nfloat * HIGH(pi/2^k)) - nfloat * LOW(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
@@ -133,7 +121,7 @@
 //
 // as follows
 //
-//    S[m] = Sin(Mpi/2^k) and C[m] = Cos(Mpi/2^k)
+//    Sm = Sin(Mpi/2^k) and Cm = Cos(Mpi/2^k)
 //    rsq = r*r
 //
 //
@@ -153,31 +141,32 @@
 //
 //       P =  r + rcub * P
 //
-//    Answer = S[m] Cos(r) + [Cm] P
+//    Answer = Sm 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 + ...
 //
-//       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
+//       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
 //
 // Then,
 //
-//    Answer = Q + C[m] P
+//    Answer = Q + Cm P
 
+#include "libm_support.h"
 
 // Registers used
 //==============================================================
 // general input registers:
-// r14 -> r26
-// r32 -> r35
+// r14 -> r19
+// r32 -> r45
 
 // predicate registers used:
-// p6 -> p11
+// p6 -> p14
 
 // floating-point registers used
 // f9 -> f15
@@ -185,94 +174,99 @@
 
 // Assembly macros
 //==============================================================
-sincos_NORM_f8                 = f9
-sincos_W                       = f10
-sincos_int_Nfloat              = f11
-sincos_Nfloat                  = f12
+sind_NORM_f8                 = f9
+sind_W                       = f10
+sind_int_Nfloat              = f11
+sind_Nfloat                  = f12
 
-sincos_r                       = f13
-sincos_rsq                     = f14
-sincos_rcub                    = f15
-sincos_save_tmp                = f15
+sind_r                       = f13
+sind_rsq                     = f14
+sind_rcub                    = f15
 
-sincos_Inv_Pi_by_16            = f32
-sincos_Pi_by_16_1              = f33
-sincos_Pi_by_16_2              = f34
+sind_Inv_Pi_by_16            = f32
+sind_Pi_by_16_hi             = f33
+sind_Pi_by_16_lo             = f34
 
-sincos_Inv_Pi_by_64            = f35
+sind_Inv_Pi_by_64            = f35
+sind_Pi_by_64_hi             = f36
+sind_Pi_by_64_lo             = f37
 
-sincos_Pi_by_16_3              = f36
+sind_Sm                      = f38
+sind_Cm                      = f39
 
-sincos_r_exact                 = f37
+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_Sm                      = f38
-sincos_Cm                      = f39
+sind_P_temp1                 = f48
+sind_P_temp2                 = f49
 
-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_Q_temp1                 = f50
+sind_Q_temp2                 = f51
 
-sincos_P_temp1                 = f48
-sincos_P_temp2                 = f49
+sind_P                       = f52
+sind_Q                       = f53
 
-sincos_Q_temp1                 = f50
-sincos_Q_temp2                 = f51
+sind_srsq                    = f54
 
-sincos_P                       = f52
-sincos_Q                       = f53
+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
 
-sincos_srsq                    = f54
+fp_tmp                       = f61
 
-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
+/////////////////////////////////////////////////////////////
 
-fp_tmp                         = f61
+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
 
-/////////////////////////////////////////////////////////////
+sind_r_exp                   = r39
+sind_r_17_ones               = r40
+
+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
+
+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
-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)
+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)
    data8 0x3EC71C963717C63A // P4
    data8 0x3EF9FFBA8F191AE6 // Q4
    data8 0xBF2A01A00F4E11A8 // P3
@@ -281,112 +275,125 @@ LOCAL_OBJECT_START(double_sincos_pq_k4)
    data8 0x3FA555555554DD45 // Q2
    data8 0xBFC5555555555555 // P1
    data8 0xBFDFFFFFFFFFFFFC // Q1
-LOCAL_OBJECT_END(double_sincos_pq_k4)
+ASM_SIZE_DIRECTIVE(double_sind_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
-LOCAL_OBJECT_END(double_sin_cos_beta_k4)
+ASM_SIZE_DIRECTIVE(double_sin_cos_beta_k4)
 
-.section .text
+.align 32
+.global sin#
+.global cos#
+#ifdef _LIBC
+.global __sin#
+.global __cos#
+#endif
 
 ////////////////////////////////////////////////////////
 // There are two entry points: sin and cos
@@ -395,374 +402,3100 @@ LOCAL_OBJECT_END(double_sin_cos_beta_k4)
 // If from sin, p8 is true
 // If from cos, p9 is true
 
-GLOBAL_IEEE754_ENTRY(sin)
+.section .text
+.proc  sin#
+#ifdef _LIBC
+.proc  __sin#
+#endif
+.align 32
+
+sin:
+#ifdef _LIBC
+__sin:
+#endif
 
 { .mlx
-      getf.exp      sincos_r_signexp    = f8
-      movl sincos_GR_sig_inv_pi_by_16   = 0xA2F9836E4E44152A // signd of 16/pi
+      alloc          r32=ar.pfs,1,13,0,0
+      movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
 }
 { .mlx
-      addl          sincos_AD_1         = @ltoff(double_sincos_pi), gp
-      movl sincos_GR_rshf_2to61         = 0x47b8000000000000 // 1.1 2^(63+63-2)
+      addl           sind_AD_1   = @ltoff(double_sind_pi), gp
+      movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
 }
 ;;
 
 { .mfi
-      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
+      ld8 sind_AD_1 = [sind_AD_1]
+      fnorm     sind_NORM_f8  = f8
+      cmp.eq     p8,p9         = r0, r0
 }
 { .mib
-      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
+      mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
+      mov            sind_r_sincos = 0x0
+      br.cond.sptk   L(SIND_SINCOS)
 }
 ;;
 
-GLOBAL_IEEE754_END(sin)
+.endp sin
+ASM_SIZE_DIRECTIVE(sin)
+
 
-GLOBAL_IEEE754_ENTRY(cos)
+.section .text
+.proc  cos#
+#ifdef _LIBC
+.proc  __cos#
+#endif
+.align 32
+cos:
+#ifdef _LIBC
+__cos:
+#endif
 
 { .mlx
-      getf.exp      sincos_r_signexp    = f8
-      movl sincos_GR_sig_inv_pi_by_16   = 0xA2F9836E4E44152A // signd of 16/pi
+      alloc          r32=ar.pfs,1,13,0,0
+      movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
 }
 { .mlx
-      addl          sincos_AD_1         = @ltoff(double_sincos_pi), gp
-      movl sincos_GR_rshf_2to61         = 0x47b8000000000000 // 1.1 2^(63+63-2)
+      addl           sind_AD_1   = @ltoff(double_sind_pi), gp
+      movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
 }
 ;;
 
 { .mfi
-      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
+      ld8 sind_AD_1 = [sind_AD_1]
+      fnorm.s1     sind_NORM_f8  = f8
+      cmp.eq     p9,p8         = r0, r0
 }
 { .mib
-      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
+      mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
+      mov            sind_r_sincos = 0x8
+      br.cond.sptk   L(SIND_SINCOS)
 }
 ;;
 
+
 ////////////////////////////////////////////////////////
 // All entry points end up here.
-// 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
+// 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
 
-///////////// Common sin and cos part //////////////////
-_SINCOS_COMMON:
+L(SIND_SINCOS):
 
 
 // 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      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
+      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
 }
 { .mlx
-      setf.d        sincos_RSHF_2TO61   = sincos_GR_rshf_2to61
-      movl          sincos_GR_rshf      = 0x43e8000000000000 // 1.1 2^63
-}                                                            // Right shift
+      setf.d sind_RSHF_2TO61 = sind_GR_rshf_2to61
+      movl sind_GR_rshf = 0x43e8000000000000 // 1.1000 2^63 for right shift
+}
 ;;
 
 // Form another constant
 //  2^-61 for scaling Nfloat
-// 0x1001a is register_bias + 27.
-// So if f8 >= 2^27, go to large argument routines
+// 0x10009 is register_bias + 10.
+// So if f8 > 2^10 = Gamma, go to DBX
 { .mfi
-      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
+      setf.exp sind_2TOM61 = sind_GR_exp_2tom61
+      fclass.m  p13,p0 = f8, 0x23           // Test for x inf
+      mov       sind_exp_limit = 0x10009
 }
 ;;
 
 // Load the two pieces of pi/16
 // Form another constant
 //  1.1000...000 * 2^63, the right shift constant
-{ .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
+{ .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
 }
 ;;
 
-_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
-};;
+{ .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
+}
 
-// Select exponent (17 lsb)
+// Start loading P, Q coefficients
+// SIN(0)
 { .mfi
-      ldfe          sincos_Pi_by_16_3   = [sincos_AD_1],16
-      nop.f         999
-      dep.z         sincos_r_exp        = sincos_r_signexp, 0, 17 
-};;
+      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
+}
+
+
+// COS(0)
+{ .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 ;;
+}
 
-// 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 >= 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
-};;
+      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 ;;
+}
 
-// sincos_W          = x * sincos_Inv_Pi_by_16
+// 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
+}
+;;
+
+// sind_W          = x * sind_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         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
-};;
+      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)
+}
+;;
+
 
-// get N = (int)sincos_int_Nfloat
-// sincos_NFLOAT = Round_Int_Nearest(sincos_W)
+// sind_NFLOAT = Round_Int_Nearest(sind_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
-      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
-};;
+      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
+}
 
-// sincos_r          = -sincos_Nfloat * sincos_Pi_by_16_1 + x
 { .mfi
-      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 
-};;
+      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
+}
 
-// Add 2^(k-1) (which is in sincos_r_sincos) to N
 { .mmi
-      add           sincos_GR_n         = sincos_GR_n, sincos_r_sincos
+      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):
+
+{ .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
+}
+;;
+
+{ .mmi
+      ld8 GR_Table_Base = [GR_Table_Base]
+      nop.m 999
+      nop.i 999
+}
 ;;
-// Get M (least k+1 bits of N)
-      and           sincos_GR_m         = 0x1f,sincos_GR_n
-      nop.i         999 
-};;
 
-// sincos_r          = sincos_r -sincos_Nfloat * sincos_Pi_by_16_2
+
+//
+//     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
-      nop.m         999
-      fnma.s1 sincos_r = sincos_NFLOAT, sincos_Pi_by_16_2,  sincos_r
-      shl           sincos_GR_32m       = sincos_GR_m,5
-};;
+       ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0
+      nop.f 999
+      nop.i 999 ;;
+}
 
-// Add 32*M to address of sin_cos_beta table
-// For sin denorm. - set uflow
 { .mfi
-      add           sincos_AD_2         = sincos_GR_32m, sincos_AD_1
-(p8)  fclass.m.unc  p10,p0              = f8,0x0b
-      nop.i         999 
-};;
+      nop.m 999
+//
+//     Load values 2**(-14) and -2**(-14)
+//
+       fcvt.fx.s1 FR_N_0_fix = FR_N_0
+      nop.i 999 ;;
+}
 
-// Load Sin and Cos table value using obtained index m  (sincosf_AD_2)
 { .mfi
-      ldfe          sincos_Sm           = [sincos_AD_2],16
-      nop.f         999 
-      nop.i         999 
-};;
+      nop.m 999
+//
+//     N_0_fix  = integer part of N_0
+//
+       fcvt.xf FR_N_0 = FR_N_0_fix
+      nop.i 999 ;;
+}
 
-// get rsq = r*r
 { .mfi
-      ldfe          sincos_Cm           = [sincos_AD_2]
-      fma.s1        sincos_rsq          = sincos_r, sincos_r,   f0 // r^2 = r*r
-      nop.i         999
+      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
-      fmpy.s0       fp_tmp              = fp_tmp,fp_tmp // forces inexact flag
-      nop.i         999 
-};;
+      nop.m 999
+       fma.s1 FR_w = FR_N_0, FR_d_1, f0
+      nop.i 999 ;;
+}
 
-// sincos_r_exact = sincos_r -sincos_Nfloat * sincos_Pi_by_16_3
 { .mfi
-      nop.m         999
-      fnma.s1 sincos_r_exact = sincos_NFLOAT, sincos_Pi_by_16_3, sincos_r
-      nop.i         999 
-};;
+      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 ;;
+}
 
-// Polynomials calculation 
-// P_1 = P4*r^2 + P3
-// Q_2 = Q4*r^2 + Q3
 { .mfi
-      nop.m         999
-      fma.s1        sincos_P_temp1      = sincos_rsq, sincos_P4, sincos_P3
-      nop.i         999
+      nop.m 999
+//
+//     N = A' * 2/pi
+//
+       fcvt.fx.s1 FR_N_fix = FR_N_float
+      nop.i 999 ;;
 }
+
 { .mfi
-      nop.m         999
-      fma.s1        sincos_Q_temp1      = sincos_rsq, sincos_Q4, sincos_Q3
-      nop.i         999 
-};;
+      nop.m 999
+//
+//     N_fix is the integer part
+//
+       fcvt.xf FR_N_float = FR_N_fix
+      nop.i 999 ;;
+}
 
-// get rcube = r^3 and S[m]*r^2
 { .mfi
-      nop.m         999
-      fmpy.s1       sincos_srsq         = sincos_Sm,sincos_rsq
-      nop.i         999
+       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 ;;
 }
+
 { .mfi
-      nop.m         999
-      fmpy.s1       sincos_rcub         = sincos_r_exact, sincos_rsq
-      nop.i         999 
-};;
+      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
+}
 
-// Polynomials calculation 
-// Q_2 = Q_1*r^2 + Q2
-// P_1 = P_1*r^2 + P2
 { .mfi
-      nop.m         999
-      fma.s1        sincos_Q_temp2      = sincos_rsq, sincos_Q_temp1, sincos_Q2
-      nop.i         999
+      nop.m 999
+       fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
+      nop.i 999 ;;
 }
+
 { .mfi
-      nop.m         999
-      fma.s1        sincos_P_temp2      = sincos_rsq, sincos_P_temp1, sincos_P2
-      nop.i         999 
-};;
+      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 ;;
+}
 
-// Polynomials calculation 
-// Q = Q_2*r^2 + Q1
-// P = P_2*r^2 + P1
 { .mfi
-      nop.m         999
-      fma.s1        sincos_Q            = sincos_rsq, sincos_Q_temp2, sincos_Q1
-      nop.i         999
+      nop.m 999
+(p9)   fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14
+      nop.i 999 ;;
 }
+
 { .mfi
-      nop.m         999
-      fma.s1        sincos_P            = sincos_rsq, sincos_P_temp2, sincos_P1
-      nop.i         999 
-};;
+      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
+}
 
-// Get final P and Q
-// Q = Q*S[m]*r^2 + S[m]
-// P = P*r^3 + r
 { .mfi
-      nop.m         999
-      fma.s1        sincos_Q            = sincos_srsq,sincos_Q, sincos_Sm
-      nop.i         999
+      nop.m 999
+(p9)   fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
+      nop.i 999 ;;
 }
+
 { .mfi
-      nop.m         999
-      fma.s1        sincos_P            = sincos_rcub,sincos_P, sincos_r_exact
-      nop.i         999 
-};;
+      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
+}
 
-// If sin(denormal), force underflow to be set
 { .mfi
-      nop.m         999
-(p10) fmpy.d.s0     fp_tmp              = sincos_NORM_f8,sincos_NORM_f8
-      nop.i         999
-};;
+      nop.m 999
+(p9)   fma.s1 FR_w = FR_N_float, FR_P_3, f0
+      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
+}
+
+{ .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 ;;
+}
+
+{ .mfi
+      nop.m 999
+(p9)   fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
+      nop.i 999 ;;
+}
+
+{ .mfi
+      nop.m 999
+(p9)   fmerge.s FR_V_hiabs = f0, FR_V_hi
+      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
+}
+
+{ .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 ;;
+}
+
+{ .mfi
+      nop.m 999
+(p12)  fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
+      nop.i 999 ;;
+}
+
+{ .mfi
+      nop.m 999
+(p8)   fma.s1 FR_c = FR_c, f1, FR_w
+      nop.i 999
+}
 
-// Final calculation
-// result = C[m]*P + Q
 { .mfb
-      nop.m         999
-      fma.d.s0      f8                  = sincos_Cm, sincos_P, sincos_Q
-      br.ret.sptk   b0  // Exit for common path
-};;
+      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
+}
 
-////////// 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
-(p8)  fma.d.s0      f8                  = f8, f0, f0 // sin(+/-0,NaN,Inf)
-      nop.i         999
+      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 ;;
 }
-// cos(+/-0) = 1.0
-// cos(Inf)  = NaN
-// cos(NaN)  = NaN
+
+{ .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 ;;
+}
+
+{ .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 ;;
+}
+
 { .mfb
-      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
-};;
+      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
+}
 
-_SINCOS_UNORM:
-// Here if x=unorm
 { .mfb
-      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
-};;
+      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):
 
-GLOBAL_IEEE754_END(cos)
+{ .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 ;;
+}
+
+{ .mfi
+      nop.m 999
+    	fma.s1	FR_rsq = FR_r, FR_r, f0
+    	extr.u	GR_i_0 = GR_N_Inc, 1, 1 ;;
+}
 
-//////////// x >= 2^27 - large arguments routine call ////////////
-LOCAL_LIBM_ENTRY(__libm_callout_sincos)
-_SINCOS_LARGE_ARGS:
-.prologue
 { .mfi
-      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
+      nop.m 999
+    	frcpa.s1 FR_r_hi, p6 = f1, FR_r
+    	cmp.eq.unc p11, p12 = 0x0, GR_i_0
 }
 ;;
 
+// ******************************************************************
+// ******************************************************************
+// ******************************************************************
+//
+//      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
-      mov           GR_SAVE_GP          = gp
-      nop.f         999
-.save b0, GR_SAVE_B0
-      mov           GR_SAVE_B0          = b0
+(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 ;;
 }
 
-.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)
+{ .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 ;;
+}
 
-{ .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)
-};;
+{ .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
-      mov           gp                  = GR_SAVE_GP
-      fma.d.s0      f8                  = f8, f1, f0 // Round result to double
-      mov           b0                  = GR_SAVE_B0
+      nop.m 999
+//
+//      Load PP_2, QQ_2
+//
+(p9)	fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3
+      nop.i 999 ;;
 }
-// Force inexact set
+
 { .mfi
-      nop.m         999
-      fmpy.s0       sincos_save_tmp     = sincos_save_tmp, sincos_save_tmp
-      nop.i         999 
+      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
+
+
+
+.proc __libm_callout_2
+__libm_callout_2:
+L(SINCOS_ARG_TOO_LARGE):
+
+.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
+}
+{ .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# ;;
 };;
 
+
+{ .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
-      nop.m         999
-      mov           ar.pfs              = GR_SAVE_PFS
-      br.ret.sptk   b0 // Exit for large arguments routine call
+        mov   gp = GR_SAVE_GP                  // Restore gp
+        mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
+        nop.b 0
 };;
 
-LOCAL_LIBM_END(__libm_callout_sincos)
 
-.type    __libm_sin_large#,@function
-.global  __libm_sin_large#
-.type    __libm_cos_large#,@function
-.global  __libm_cos_large#
+{ .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_double_dbx#,@function
+.global __libm_sin_double_dbx#
+.type __libm_cos_double_dbx#,@function
+.global __libm_cos_double_dbx#