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