diff options
Diffstat (limited to 'sysdeps/ia64/fpu/s_expm1l.S')
| -rw-r--r-- | sysdeps/ia64/fpu/s_expm1l.S | 1947 |
1 files changed, 879 insertions, 1068 deletions
diff --git a/sysdeps/ia64/fpu/s_expm1l.S b/sysdeps/ia64/fpu/s_expm1l.S index e53d3c8..5f135fa 100644 --- a/sysdeps/ia64/fpu/s_expm1l.S +++ b/sysdeps/ia64/fpu/s_expm1l.S @@ -1,10 +1,10 @@ -.file "exp_m1l.s" +.file "expl_m1.s" -// Copyright (C) 2000, 2001, Intel Corporation + +// Copyright (c) 2000 - 2003, 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,15 +35,23 @@ // // 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 //============================================================== -// 4/04/00 Unwind support added -// 8/15/00 Bundle added after call to __libm_error_support to properly +// 02/02/00 Initial Version +// 04/04/00 Unwind support added +// 08/15/00 Bundle added after call to __libm_error_support to properly // set [the previously overwritten] GR_Parameter_RESULT. +// 07/07/01 Improved speed of all paths +// 05/20/02 Cleaned up namespace and sf0 syntax +// 02/10/03 Reordered header: .section, .global, .proc, .align; +// used data8 for long double table values +// 03/11/03 Improved accuracy and performance, corrected missing inexact flags +// 04/17/03 Eliminated misplaced and unused data label +// 12/15/03 Eliminated call to error support on expm1l underflow // -// ********************************************************************* +//********************************************************************* // // Function: Combined expl(x) and expm1l(x), where // x @@ -51,20 +59,20 @@ // x // expm1l(x) = e - 1 for double-extended precision x values // -// ********************************************************************* +//********************************************************************* // // Resources Used: // // Floating-Point Registers: f8 (Input and Return Value) -// f9,f32-f61, f99-f102 +// f9-f15,f32-f77 // // General Purpose Registers: -// r32-r61 -// r62-r65 (Used to pass arguments to error handling routine) +// r14-r38 +// r35-r38 (Used to pass arguments to error handling routine) // // Predicate Registers: p6-p15 // -// ********************************************************************* +//********************************************************************* // // IEEE Special Conditions: // @@ -74,39 +82,37 @@ // (Error Handling Routine called for overflow and Underflow) // Inexact raised when appropriate by algorithm // -// expl(inf) = inf -// expl(-inf) = +0 -// expl(SNaN) = QNaN -// expl(QNaN) = QNaN -// expl(0) = 1 -// expl(EM_special Values) = QNaN -// expl(inf) = inf -// expm1l(-inf) = -1 -// expm1l(SNaN) = QNaN -// expm1l(QNaN) = QNaN -// expm1l(0) = 0 -// expm1l(EM_special Values) = QNaN +// exp(inf) = inf +// exp(-inf) = +0 +// exp(SNaN) = QNaN +// exp(QNaN) = QNaN +// exp(0) = 1 +// exp(EM_special Values) = QNaN +// exp(inf) = inf +// expm1(-inf) = -1 +// expm1(SNaN) = QNaN +// expm1(QNaN) = QNaN +// expm1(0) = 0 +// expm1(EM_special Values) = QNaN // -// ********************************************************************* +//********************************************************************* // // Implementation and Algorithm Notes: // // ker_exp_64( in_FR : X, -// in_GR : Flag, -// in_GR : Expo_Range // out_FR : Y_hi, // out_FR : Y_lo, // out_FR : scale, // out_PR : Safe ) // -// On input, X is in register format and -// Flag = 0 for exp, -// Flag = 1 for expm1, +// On input, X is in register format +// p6 for exp, +// p7 for expm1, // -// On output, provided X and X_cor are real numbers, then +// On output, // -// scale*(Y_hi + Y_lo) approximates expl(X) if Flag is 0 -// scale*(Y_hi + Y_lo) approximates expl(X)-1 if Flag is 1 +// scale*(Y_hi + Y_lo) approximates exp(X) if exp +// scale*(Y_hi + Y_lo) approximates exp(X)-1 if expm1 // // The accuracy is sufficient for a highly accurate 64 sig. // bit implementation. Safe is set if there is no danger of @@ -122,36 +128,36 @@ // The method consists of three cases. // // If |X| < Tiny use case exp_tiny; -// else if |X| < 2^(-6) use case exp_small; +// else if |X| < 2^(-m) use case exp_small; m=12 for exp, m=7 for expm1 // else use case exp_regular; // // Case exp_tiny: // -// 1 + X can be used to approximate expl(X) or expl(X+X_cor); -// X + X^2/2 can be used to approximate expl(X) - 1 +// 1 + X can be used to approximate exp(X) +// X + X^2/2 can be used to approximate exp(X) - 1 // // Case exp_small: // -// Here, expl(X), expl(X+X_cor), and expl(X) - 1 can all be +// Here, exp(X) and exp(X) - 1 can all be // appproximated by a relatively simple polynomial. // // This polynomial resembles the truncated Taylor series // -// expl(w) = 1 + w + w^2/2! + w^3/3! + ... + w^n/n! +// exp(w) = 1 + w + w^2/2! + w^3/3! + ... + w^n/n! // // Case exp_regular: // // Here we use a table lookup method. The basic idea is that in -// order to compute expl(X), we accurately decompose X into +// order to compute exp(X), we accurately decompose X into // // X = N * log(2)/(2^12) + r, |r| <= log(2)/2^13. // // Hence // -// expl(X) = 2^( N / 2^12 ) * expl(r). +// exp(X) = 2^( N / 2^12 ) * exp(r). // // The value 2^( N / 2^12 ) is obtained by simple combinations -// of values calculated beforehand and stored in table; expl(r) +// of values calculated beforehand and stored in table; exp(r) // is approximated by a short polynomial because |r| is small. // // We elaborate this method in 4 steps. @@ -178,13 +184,9 @@ // as a double-precision number; L_lo has 64 significant bits and // stored as a double-extended number. // -// In the case Flag = 2, we further modify r by -// -// r := r + X_cor. -// // Step 2: Approximation // -// expl(r) - 1 is approximated by a short polynomial of the form +// exp(r) - 1 is approximated by a short polynomial of the form // // r + A_1 r^2 + A_2 r^3 + A_3 r^4 . // @@ -213,19 +215,19 @@ // Define two mathematical values, delta_1 and delta_2, implicitly // such that // -// T_1 = expl( [M_1 log(2)/2^6] - delta_1 ) -// T_2 = expl( [M_2 log(2)/2^12] - delta_2 ) +// T_1 = exp( [M_1 log(2)/2^6] - delta_1 ) +// T_2 = exp( [M_2 log(2)/2^12] - delta_2 ) // // are representable as 24 significant bits. To illustrate the idea, // we show how we define delta_1: // -// T_1 := round_to_24_bits( expl( M_1 log(2)/2^6 ) ) +// T_1 := round_to_24_bits( exp( M_1 log(2)/2^6 ) ) // delta_1 = (M_1 log(2)/2^6) - log( T_1 ) // // The last equality means mathematical equality. We then tabulate // -// W_1 := expl(delta_1) - 1 -// W_2 := expl(delta_2) - 1 +// W_1 := exp(delta_1) - 1 +// W_2 := exp(delta_2) - 1 // // Both in double precision. // @@ -235,13 +237,13 @@ // T := T_1 * T_2 ...exactly // W := W_1 + (1 + W_1)*W_2 // -// W approximates expl( delta ) - 1 where delta = delta_1 + delta_2. +// W approximates exp( delta ) - 1 where delta = delta_1 + delta_2. // The mathematical product of T and (W+1) is an accurate representation // of 2^(M_1/2^6) * 2^(M_2/2^12). // // Step 4. Reconstruction // -// Finally, we can reconstruct expl(X), expl(X) - 1. +// Finally, we can reconstruct exp(X), exp(X) - 1. // Because // // X = K * log(2) + (M_1*log(2)/2^6 - delta_1) @@ -249,18 +251,18 @@ // + delta_1 + delta_2 + r ...accurately // We have // -// expl(X) ~=~ 2^K * ( T + T*[expl(delta_1+delta_2+r) - 1] ) -// ~=~ 2^K * ( T + T*[expl(delta + r) - 1] ) -// ~=~ 2^K * ( T + T*[(expl(delta)-1) -// + expl(delta)*(expl(r)-1)] ) +// exp(X) ~=~ 2^K * ( T + T*[exp(delta_1+delta_2+r) - 1] ) +// ~=~ 2^K * ( T + T*[exp(delta + r) - 1] ) +// ~=~ 2^K * ( T + T*[(exp(delta)-1) +// + exp(delta)*(exp(r)-1)] ) // ~=~ 2^K * ( T + T*( W + (1+W)*poly(r) ) ) // ~=~ 2^K * ( Y_hi + Y_lo ) // // where Y_hi = T and Y_lo = T*(W + (1+W)*poly(r)) // -// For expl(X)-1, we have +// For exp(X)-1, we have // -// expl(X)-1 ~=~ 2^K * ( Y_hi + Y_lo ) - 1 +// exp(X)-1 ~=~ 2^K * ( Y_hi + Y_lo ) - 1 // ~=~ 2^K * ( Y_hi + Y_lo - 2^(-K) ) // // and we combine Y_hi + Y_lo - 2^(-N) into the form of two @@ -278,7 +280,7 @@ // different rounding directions and a correct setting of the SAFE // flag. // -// If Flag is 1, then +// If expm1 is 1, then // SAFE := False ...possibility of underflow // Scale := 1.0 // Y_hi := X @@ -296,26 +298,25 @@ // // Let r = X // -// If Flag is not 1 ...i.e. expl( argument ) +// If exp ...i.e. exp( argument ) // // rsq := r * r; // r4 := rsq*rsq // poly_lo := P_3 + r*(P_4 + r*(P_5 + r*P_6)) // poly_hi := r + rsq*(P_1 + r*P_2) // Y_lo := poly_hi + r4 * poly_lo -// set lsb(Y_lo) to 1 // Y_hi := 1.0 // Scale := 1.0 // -// Else ...i.e. expl( argument ) - 1 +// Else ...i.e. exp( argument ) - 1 // // rsq := r * r // r4 := rsq * rsq -// r6 := rsq * r4 -// poly_lo := r6*(Q_5 + r*(Q_6 + r*Q_7)) -// poly_hi := Q_1 + r*(Q_2 + r*(Q_3 + r*Q_4)) -// Y_lo := rsq*poly_hi + poly_lo -// set lsb(Y_lo) to 1 +// poly_lo := Q_7 + r*(Q_8 + r*Q_9)) +// poly_med:= Q_3 + r*Q_4 + rsq*(Q_5 + r*Q_6) +// poly_med:= poly_med + r4*poly_lo +// poly_hi := Q_1 + r*Q_2 +// Y_lo := rsq*(poly_hi + rsq*poly_lo) // Y_hi := X // Scale := 1.0 // @@ -325,14 +326,14 @@ // // The previous description contain enough information except the // computation of poly and the final Y_hi and Y_lo in the case for -// expl(X)-1. +// exp(X)-1. // // The computation of poly for Step 2: // // rsq := r*r // poly := r + rsq*(A_1 + r*(A_2 + r*A_3)) // -// For the case expl(X) - 1, we need to incorporate 2^(-K) into +// For the case exp(X) - 1, we need to incorporate 2^(-K) into // Y_hi and Y_lo at the end of Step 4. // // If K > 10 then @@ -346,72 +347,197 @@ // End If // End If // +//======================================================= +// General Purpose Registers +// +GR_ad_Arg = r14 +GR_ad_A = r15 +GR_sig_inv_ln2 = r15 +GR_rshf_2to51 = r16 +GR_ad_PQ = r16 +GR_ad_Q = r16 +GR_signexp_x = r17 +GR_exp_x = r17 +GR_small_exp = r18 +GR_rshf = r18 +GR_exp_mask = r19 +GR_ad_W1 = r20 +GR_exp_2tom51 = r20 +GR_ad_W2 = r21 +GR_exp_underflow = r21 +GR_M2 = r22 +GR_huge_exp = r22 +GR_M1 = r23 +GR_huge_signif = r23 +GR_K = r24 +GR_one = r24 +GR_minus_one = r24 +GR_exp_bias = r25 +GR_ad_Limits = r26 +GR_N_fix = r26 +GR_exp_2_mk = r26 +GR_ad_P = r27 +GR_exp_2_k = r27 +GR_big_expo_neg = r28 +GR_very_small_exp = r29 +GR_exp_half = r29 +GR_ad_T1 = r30 +GR_ad_T2 = r31 -#include "libm_support.h" +GR_SAVE_PFS = r32 +GR_SAVE_B0 = r33 +GR_SAVE_GP = r34 +GR_Parameter_X = r35 +GR_Parameter_Y = r36 +GR_Parameter_RESULT = r37 +GR_Parameter_TAG = r38 -#ifdef _LIBC -.rodata -#else -.data -#endif +// Floating Point Registers +// +FR_norm_x = f9 +FR_RSHF_2TO51 = f10 +FR_INV_LN2_2TO63 = f11 +FR_W_2TO51_RSH = f12 +FR_2TOM51 = f13 +FR_RSHF = f14 +FR_Y_hi = f34 +FR_Y_lo = f35 +FR_scale = f36 +FR_tmp = f37 +FR_float_N = f38 +FR_N_signif = f39 +FR_L_hi = f40 +FR_L_lo = f41 +FR_r = f42 +FR_W1 = f43 +FR_T1 = f44 +FR_W2 = f45 +FR_T2 = f46 +FR_W1_p1 = f47 +FR_rsq = f48 +FR_A2 = f49 +FR_r4 = f50 +FR_A3 = f51 +FR_poly = f52 +FR_T = f53 +FR_W = f54 +FR_Wp1 = f55 +FR_p21 = f59 +FR_p210 = f59 +FR_p65 = f60 +FR_p654 = f60 +FR_p6543 = f60 +FR_2_mk = f61 +FR_P4Q7 = f61 +FR_P4 = f61 +FR_Q7 = f61 +FR_P3Q6 = f62 +FR_P3 = f62 +FR_Q6 = f62 +FR_q65 = f62 +FR_q6543 = f62 +FR_P2Q5 = f63 +FR_P2 = f63 +FR_Q5 = f63 +FR_P1Q4 = f64 +FR_P1 = f64 +FR_Q4 = f64 +FR_q43 = f64 +FR_Q3 = f65 +FR_Q2 = f66 +FR_q21 = f66 +FR_Q1 = f67 +FR_A1 = f68 +FR_P6Q9 = f68 +FR_P6 = f68 +FR_Q9 = f68 +FR_P5Q8 = f69 +FR_P5 = f69 +FR_Q8 = f69 +FR_q987 = f69 +FR_q98 = f69 +FR_q9876543 = f69 +FR_min_oflow_x = f70 +FR_huge_exp = f70 +FR_zero_uflow_x = f71 +FR_huge_signif = f71 +FR_huge = f72 +FR_small = f72 +FR_half = f73 +FR_T_scale = f74 +FR_result_lo = f75 +FR_W_T_scale = f76 +FR_Wp1_T_scale = f77 +FR_ftz = f77 +FR_half_x = f77 +// -.align 64 -Constants_exp_64_Arg: -ASM_TYPE_DIRECTIVE(Constants_exp_64_Arg,@object) -data4 0x5C17F0BC,0xB8AA3B29,0x0000400B,0x00000000 -data4 0x00000000,0xB17217F4,0x00003FF2,0x00000000 -data4 0xF278ECE6,0xF473DE6A,0x00003FD4,0x00000000 -// /* Inv_L, L_hi, L_lo */ -ASM_SIZE_DIRECTIVE(Constants_exp_64_Arg) +FR_X = f9 +FR_Y = f0 +FR_RESULT = f15 -.align 64 -Constants_exp_64_Exponents: -ASM_TYPE_DIRECTIVE(Constants_exp_64_Exponents,@object) -data4 0x0000007E,0x00000000,0xFFFFFF83,0xFFFFFFFF -data4 0x000003FE,0x00000000,0xFFFFFC03,0xFFFFFFFF -data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF -data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF -data4 0xFFFFFFE2,0xFFFFFFFF,0xFFFFFFC4,0xFFFFFFFF -data4 0xFFFFFFBA,0xFFFFFFFF,0xFFFFFFBA,0xFFFFFFFF -ASM_SIZE_DIRECTIVE(Constants_exp_64_Exponents) +// ************* DO NOT CHANGE ORDER OF THESE TABLES ******************** -.align 64 -Constants_exp_64_A: -ASM_TYPE_DIRECTIVE(Constants_exp_64_A,@object) -data4 0xB1B736A0,0xAAAAAAAB,0x00003FFA,0x00000000 -data4 0x90CD6327,0xAAAAAAAB,0x00003FFC,0x00000000 -data4 0xFFFFFFFF,0xFFFFFFFF,0x00003FFD,0x00000000 -// /* Reversed */ -ASM_SIZE_DIRECTIVE(Constants_exp_64_A) +// double-extended 1/ln(2) +// 3fff b8aa 3b29 5c17 f0bb be87fed0691d3e88 +// 3fff b8aa 3b29 5c17 f0bc +// For speed the significand will be loaded directly with a movl and setf.sig +// and the exponent will be bias+63 instead of bias+0. Thus subsequent +// computations need to scale appropriately. +// The constant 2^12/ln(2) is needed for the computation of N. This is also +// obtained by scaling the computations. +// +// Two shifting constants are loaded directly with movl and setf.d. +// 1. RSHF_2TO51 = 1.1000..00 * 2^(63-12) +// This constant is added to x*1/ln2 to shift the integer part of +// x*2^12/ln2 into the rightmost bits of the significand. +// The result of this fma is N_signif. +// 2. RSHF = 1.1000..00 * 2^(63) +// This constant is subtracted from N_signif * 2^(-51) to give +// the integer part of N, N_fix, as a floating-point number. +// The result of this fms is float_N. +RODATA .align 64 -Constants_exp_64_P: -ASM_TYPE_DIRECTIVE(Constants_exp_64_P,@object) -data4 0x43914A8A,0xD00D6C81,0x00003FF2,0x00000000 -data4 0x30304B30,0xB60BC4AC,0x00003FF5,0x00000000 -data4 0x7474C518,0x88888888,0x00003FF8,0x00000000 -data4 0x8DAE729D,0xAAAAAAAA,0x00003FFA,0x00000000 -data4 0xAAAAAF61,0xAAAAAAAA,0x00003FFC,0x00000000 -data4 0x000004C7,0x80000000,0x00003FFE,0x00000000 -// /* Reversed */ -ASM_SIZE_DIRECTIVE(Constants_exp_64_P) +LOCAL_OBJECT_START(Constants_exp_64_Arg) +//data8 0xB8AA3B295C17F0BC,0x0000400B // Inv_L = 2^12/log(2) +data8 0xB17217F400000000,0x00003FF2 // L_hi = hi part log(2)/2^12 +data8 0xF473DE6AF278ECE6,0x00003FD4 // L_lo = lo part log(2)/2^12 +LOCAL_OBJECT_END(Constants_exp_64_Arg) -.align 64 -Constants_exp_64_Q: -ASM_TYPE_DIRECTIVE(Constants_exp_64_Q,@object) -data4 0xA49EF6CA,0xD00D56F7,0x00003FEF,0x00000000 -data4 0x1C63493D,0xD00D59AB,0x00003FF2,0x00000000 -data4 0xFB50CDD2,0xB60B60B5,0x00003FF5,0x00000000 -data4 0x7BA68DC8,0x88888888,0x00003FF8,0x00000000 -data4 0xAAAAAC8D,0xAAAAAAAA,0x00003FFA,0x00000000 -data4 0xAAAAACCA,0xAAAAAAAA,0x00003FFC,0x00000000 -data4 0x00000000,0x80000000,0x00003FFE,0x00000000 -// /* Reversed */ -ASM_SIZE_DIRECTIVE(Constants_exp_64_Q) +LOCAL_OBJECT_START(Constants_exp_64_Limits) +data8 0xb17217f7d1cf79ac,0x0000400c // Smallest long dbl oflow x +data8 0xb220000000000000,0x0000c00c // Small long dbl uflow zero x +LOCAL_OBJECT_END(Constants_exp_64_Limits) -.align 64 -Constants_exp_64_T1: -ASM_TYPE_DIRECTIVE(Constants_exp_64_T1,@object) +LOCAL_OBJECT_START(Constants_exp_64_A) +data8 0xAAAAAAABB1B736A0,0x00003FFA // A3 +data8 0xAAAAAAAB90CD6327,0x00003FFC // A2 +data8 0xFFFFFFFFFFFFFFFF,0x00003FFD // A1 +LOCAL_OBJECT_END(Constants_exp_64_A) + +LOCAL_OBJECT_START(Constants_exp_64_P) +data8 0xD00D6C8143914A8A,0x00003FF2 // P6 +data8 0xB60BC4AC30304B30,0x00003FF5 // P5 +data8 0x888888887474C518,0x00003FF8 // P4 +data8 0xAAAAAAAA8DAE729D,0x00003FFA // P3 +data8 0xAAAAAAAAAAAAAF61,0x00003FFC // P2 +data8 0x80000000000004C7,0x00003FFE // P1 +LOCAL_OBJECT_END(Constants_exp_64_P) + +LOCAL_OBJECT_START(Constants_exp_64_Q) +data8 0x93F2AC5F7471F32E, 0x00003FE9 // Q9 +data8 0xB8DA0F3550B3E764, 0x00003FEC // Q8 +data8 0xD00D00D0028E89C4, 0x00003FEF // Q7 +data8 0xD00D00DAEB8C4E91, 0x00003FF2 // Q6 +data8 0xB60B60B60B60B6F5, 0x00003FF5 // Q5 +data8 0x888888888886CC23, 0x00003FF8 // Q4 +data8 0xAAAAAAAAAAAAAAAB, 0x00003FFA // Q3 +data8 0xAAAAAAAAAAAAAAAB, 0x00003FFC // Q2 +data8 0x8000000000000000, 0x00003FFE // Q1 +LOCAL_OBJECT_END(Constants_exp_64_Q) + +LOCAL_OBJECT_START(Constants_exp_64_T1) data4 0x3F800000,0x3F8164D2,0x3F82CD87,0x3F843A29 data4 0x3F85AAC3,0x3F871F62,0x3F88980F,0x3F8A14D5 data4 0x3F8B95C2,0x3F8D1ADF,0x3F8EA43A,0x3F9031DC @@ -428,11 +554,9 @@ data4 0x3FD744FD,0x3FD99D16,0x3FDBFBB8,0x3FDE60F5 data4 0x3FE0CCDF,0x3FE33F89,0x3FE5B907,0x3FE8396A data4 0x3FEAC0C7,0x3FED4F30,0x3FEFE4BA,0x3FF28177 data4 0x3FF5257D,0x3FF7D0DF,0x3FFA83B3,0x3FFD3E0C -ASM_SIZE_DIRECTIVE(Constants_exp_64_T1) +LOCAL_OBJECT_END(Constants_exp_64_T1) -.align 64 -Constants_exp_64_T2: -ASM_TYPE_DIRECTIVE(Constants_exp_64_T2,@object) +LOCAL_OBJECT_START(Constants_exp_64_T2) data4 0x3F800000,0x3F80058C,0x3F800B18,0x3F8010A4 data4 0x3F801630,0x3F801BBD,0x3F80214A,0x3F8026D7 data4 0x3F802C64,0x3F8031F2,0x3F803780,0x3F803D0E @@ -449,1124 +573,812 @@ data4 0x3F810B41,0x3F8110D8,0x3F81166F,0x3F811C07 data4 0x3F81219F,0x3F812737,0x3F812CD0,0x3F813269 data4 0x3F813802,0x3F813D9B,0x3F814334,0x3F8148CE data4 0x3F814E68,0x3F815402,0x3F81599C,0x3F815F37 -ASM_SIZE_DIRECTIVE(Constants_exp_64_T2) +LOCAL_OBJECT_END(Constants_exp_64_T2) -.align 64 -Constants_exp_64_W1: -ASM_TYPE_DIRECTIVE(Constants_exp_64_W1,@object) -data4 0x00000000,0x00000000,0x171EC4B4,0xBE384454 -data4 0x4AA72766,0xBE694741,0xD42518F8,0xBE5D32B6 -data4 0x3A319149,0x3E68D96D,0x62415F36,0xBE68F4DA -data4 0xC9C86A3B,0xBE6DDA2F,0xF49228FE,0x3E6B2E50 -data4 0x1188B886,0xBE49C0C2,0x1A4C2F1F,0x3E64BFC2 -data4 0x2CB98B54,0xBE6A2FBB,0x9A55D329,0x3E5DC5DE -data4 0x39A7AACE,0x3E696490,0x5C66DBA5,0x3E54728B -data4 0xBA1C7D7D,0xBE62B0DB,0x09F1AF5F,0x3E576E04 -data4 0x1A0DD6A1,0x3E612500,0x795FBDEF,0xBE66A419 -data4 0xE1BD41FC,0xBE5CDE8C,0xEA54964F,0xBE621376 -data4 0x476E76EE,0x3E6370BE,0x3427EB92,0x3E390D1A -data4 0x2BF82BF8,0x3E1336DE,0xD0F7BD9E,0xBE5FF1CB -data4 0x0CEB09DD,0xBE60A355,0x0980F30D,0xBE5CA37E -data4 0x4C082D25,0xBE5C541B,0x3B467D29,0xBE5BBECA -data4 0xB9D946C5,0xBE400D8A,0x07ED374A,0xBE5E2A08 -data4 0x365C8B0A,0xBE66CB28,0xD3403BCA,0x3E3AAD5B -data4 0xC7EA21E0,0x3E526055,0xE72880D6,0xBE442C75 -data4 0x85222A43,0x3E58B2BB,0x522C42BF,0xBE5AAB79 -data4 0x469DC2BC,0xBE605CB4,0xA48C40DC,0xBE589FA7 -data4 0x1AA42614,0xBE51C214,0xC37293F4,0xBE48D087 -data4 0xA2D673E0,0x3E367A1C,0x114F7A38,0xBE51BEBB -data4 0x661A4B48,0xBE6348E5,0x1D3B9962,0xBDF52643 -data4 0x35A78A53,0x3E3A3B5E,0x1CECD788,0xBE46C46C -data4 0x7857D689,0xBE60B7EC,0xD14F1AD7,0xBE594D3D -data4 0x4C9A8F60,0xBE4F9C30,0x02DFF9D2,0xBE521873 -data4 0x55E6D68F,0xBE5E4C88,0x667F3DC4,0xBE62140F -data4 0x3BF88747,0xBE36961B,0xC96EC6AA,0x3E602861 -data4 0xD57FD718,0xBE3B5151,0xFC4A627B,0x3E561CD0 -data4 0xCA913FEA,0xBE3A5217,0x9A5D193A,0x3E40A3CC -data4 0x10A9C312,0xBE5AB713,0xC5F57719,0x3E4FDADB -data4 0xDBDF59D5,0x3E361428,0x61B4180D,0x3E5DB5DB -data4 0x7408D856,0xBE42AD5F,0x31B2B707,0x3E2A3148 -ASM_SIZE_DIRECTIVE(Constants_exp_64_W1) +LOCAL_OBJECT_START(Constants_exp_64_W1) +data8 0x0000000000000000, 0xBE384454171EC4B4 +data8 0xBE6947414AA72766, 0xBE5D32B6D42518F8 +data8 0x3E68D96D3A319149, 0xBE68F4DA62415F36 +data8 0xBE6DDA2FC9C86A3B, 0x3E6B2E50F49228FE +data8 0xBE49C0C21188B886, 0x3E64BFC21A4C2F1F +data8 0xBE6A2FBB2CB98B54, 0x3E5DC5DE9A55D329 +data8 0x3E69649039A7AACE, 0x3E54728B5C66DBA5 +data8 0xBE62B0DBBA1C7D7D, 0x3E576E0409F1AF5F +data8 0x3E6125001A0DD6A1, 0xBE66A419795FBDEF +data8 0xBE5CDE8CE1BD41FC, 0xBE621376EA54964F +data8 0x3E6370BE476E76EE, 0x3E390D1A3427EB92 +data8 0x3E1336DE2BF82BF8, 0xBE5FF1CBD0F7BD9E +data8 0xBE60A3550CEB09DD, 0xBE5CA37E0980F30D +data8 0xBE5C541B4C082D25, 0xBE5BBECA3B467D29 +data8 0xBE400D8AB9D946C5, 0xBE5E2A0807ED374A +data8 0xBE66CB28365C8B0A, 0x3E3AAD5BD3403BCA +data8 0x3E526055C7EA21E0, 0xBE442C75E72880D6 +data8 0x3E58B2BB85222A43, 0xBE5AAB79522C42BF +data8 0xBE605CB4469DC2BC, 0xBE589FA7A48C40DC +data8 0xBE51C2141AA42614, 0xBE48D087C37293F4 +data8 0x3E367A1CA2D673E0, 0xBE51BEBB114F7A38 +data8 0xBE6348E5661A4B48, 0xBDF526431D3B9962 +data8 0x3E3A3B5E35A78A53, 0xBE46C46C1CECD788 +data8 0xBE60B7EC7857D689, 0xBE594D3DD14F1AD7 +data8 0xBE4F9C304C9A8F60, 0xBE52187302DFF9D2 +data8 0xBE5E4C8855E6D68F, 0xBE62140F667F3DC4 +data8 0xBE36961B3BF88747, 0x3E602861C96EC6AA +data8 0xBE3B5151D57FD718, 0x3E561CD0FC4A627B +data8 0xBE3A5217CA913FEA, 0x3E40A3CC9A5D193A +data8 0xBE5AB71310A9C312, 0x3E4FDADBC5F57719 +data8 0x3E361428DBDF59D5, 0x3E5DB5DB61B4180D +data8 0xBE42AD5F7408D856, 0x3E2A314831B2B707 +LOCAL_OBJECT_END(Constants_exp_64_W1) -.align 64 -Constants_exp_64_W2: -ASM_TYPE_DIRECTIVE(Constants_exp_64_W2,@object) -data4 0x00000000,0x00000000,0x37A3D7A2,0xBE641F25 -data4 0xAD028C40,0xBE68DD57,0xF212B1B6,0xBE5C77D8 -data4 0x1BA5B070,0x3E57878F,0x2ECAE6FE,0xBE55A36A -data4 0x569DFA3B,0xBE620608,0xA6D300A3,0xBE53B50E -data4 0x223F8F2C,0x3E5B5EF2,0xD6DE0DF4,0xBE56A0D9 -data4 0xEAE28F51,0xBE64EEF3,0x367EA80B,0xBE5E5AE2 -data4 0x5FCBC02D,0x3E47CB1A,0x9BDAFEB7,0xBE656BA0 -data4 0x805AFEE7,0x3E6E70C6,0xA3415EBA,0xBE6E0509 -data4 0x49BFF529,0xBE56856B,0x00508651,0x3E66DD33 -data4 0xC114BC13,0x3E51165F,0xC453290F,0x3E53333D -data4 0x05539FDA,0x3E6A072B,0x7C0A7696,0xBE47CD87 -data4 0xEB05C6D9,0xBE668BF4,0x6AE86C93,0xBE67C3E3 -data4 0xD0B3E84B,0xBE533904,0x556B53CE,0x3E63E8D9 -data4 0x63A98DC8,0x3E212C89,0x032A7A22,0xBE33138F -data4 0xBC584008,0x3E530FA9,0xCCB93C97,0xBE6ADF82 -data4 0x8370EA39,0x3E5F9113,0xFB6A05D8,0x3E5443A4 -data4 0x181FEE7A,0x3E63DACD,0xF0F67DEC,0xBE62B29D -data4 0x3DDE6307,0x3E65C483,0xD40A24C1,0x3E5BF030 -data4 0x14E437BE,0x3E658B8F,0xED98B6C7,0xBE631C29 -data4 0x04CF7C71,0x3E6335D2,0xE954A79D,0x3E529EED -data4 0xF64A2FB8,0x3E5D9257,0x854ED06C,0xBE6BED1B -data4 0xD71405CB,0x3E5096F6,0xACB9FDF5,0xBE3D4893 -data4 0x01B68349,0xBDFEB158,0xC6A463B9,0x3E628D35 -data4 0xADE45917,0xBE559725,0x042FC476,0xBE68C29C -data4 0x01E511FA,0xBE67593B,0x398801ED,0xBE4A4313 -data4 0xDA7C3300,0x3E699571,0x08062A9E,0x3E5349BE -data4 0x755BB28E,0x3E5229C4,0x77A1F80D,0x3E67E426 -data4 0x6B69C352,0xBE52B33F,0x084DA57F,0xBE6B3550 -data4 0xD1D09A20,0xBE6DB03F,0x2161B2C1,0xBE60CBC4 -data4 0x78A2B771,0x3E56ED9C,0x9D0FA795,0xBE508E31 -data4 0xFD1A54E9,0xBE59482A,0xB07FD23E,0xBE2A17CE -data4 0x17365712,0x3E68BF5C,0xB3785569,0x3E3956F9 -ASM_SIZE_DIRECTIVE(Constants_exp_64_W2) - -GR_SAVE_PFS = r59 -GR_SAVE_B0 = r60 -GR_SAVE_GP = r61 -GR_Parameter_X = r62 -GR_Parameter_Y = r63 -GR_Parameter_RESULT = r64 -GR_Parameter_TAG = r65 +LOCAL_OBJECT_START(Constants_exp_64_W2) +data8 0x0000000000000000, 0xBE641F2537A3D7A2 +data8 0xBE68DD57AD028C40, 0xBE5C77D8F212B1B6 +data8 0x3E57878F1BA5B070, 0xBE55A36A2ECAE6FE +data8 0xBE620608569DFA3B, 0xBE53B50EA6D300A3 +data8 0x3E5B5EF2223F8F2C, 0xBE56A0D9D6DE0DF4 +data8 0xBE64EEF3EAE28F51, 0xBE5E5AE2367EA80B +data8 0x3E47CB1A5FCBC02D, 0xBE656BA09BDAFEB7 +data8 0x3E6E70C6805AFEE7, 0xBE6E0509A3415EBA +data8 0xBE56856B49BFF529, 0x3E66DD3300508651 +data8 0x3E51165FC114BC13, 0x3E53333DC453290F +data8 0x3E6A072B05539FDA, 0xBE47CD877C0A7696 +data8 0xBE668BF4EB05C6D9, 0xBE67C3E36AE86C93 +data8 0xBE533904D0B3E84B, 0x3E63E8D9556B53CE +data8 0x3E212C8963A98DC8, 0xBE33138F032A7A22 +data8 0x3E530FA9BC584008, 0xBE6ADF82CCB93C97 +data8 0x3E5F91138370EA39, 0x3E5443A4FB6A05D8 +data8 0x3E63DACD181FEE7A, 0xBE62B29DF0F67DEC +data8 0x3E65C4833DDE6307, 0x3E5BF030D40A24C1 +data8 0x3E658B8F14E437BE, 0xBE631C29ED98B6C7 +data8 0x3E6335D204CF7C71, 0x3E529EEDE954A79D +data8 0x3E5D9257F64A2FB8, 0xBE6BED1B854ED06C +data8 0x3E5096F6D71405CB, 0xBE3D4893ACB9FDF5 +data8 0xBDFEB15801B68349, 0x3E628D35C6A463B9 +data8 0xBE559725ADE45917, 0xBE68C29C042FC476 +data8 0xBE67593B01E511FA, 0xBE4A4313398801ED +data8 0x3E699571DA7C3300, 0x3E5349BE08062A9E +data8 0x3E5229C4755BB28E, 0x3E67E42677A1F80D +data8 0xBE52B33F6B69C352, 0xBE6B3550084DA57F +data8 0xBE6DB03FD1D09A20, 0xBE60CBC42161B2C1 +data8 0x3E56ED9C78A2B771, 0xBE508E319D0FA795 +data8 0xBE59482AFD1A54E9, 0xBE2A17CEB07FD23E +data8 0x3E68BF5C17365712, 0x3E3956F9B3785569 +LOCAL_OBJECT_END(Constants_exp_64_W2) -FR_X = f9 -FR_Y = f9 -FR_RESULT = f99 .section .text -.proc expm1l# -.global expm1l# -.align 64 -expm1l: -#ifdef _LIBC -.global __expm1l# -__expm1l: -#endif -{ .mii -alloc r32 = ar.pfs,0,30,4,0 -(p0) add r33 = 1, r0 -(p0) cmp.eq.unc p7, p0 = r0, r0 -} -{ .mbb - nop.m 999 -(p0) br.cond.sptk exp_continue - nop.b 999 ;; -} + +GLOBAL_IEEE754_ENTRY(expm1l) // -// Set p7 true for expm1 -// Set Flag = r33 = 1 for expm1 +// Set p7 true for expm1, p6 false // -.endp expm1l -ASM_SIZE_DIRECTIVE(expm1l) - -#ifdef _LIBC -libm_hidden_def (__expm1l) -#endif - -.section .text -.proc expl# -.global expl# -.align 64 -expl: -#ifdef _LIBC -.global __ieee754_expl# -__ieee754_expl: -#endif -{ .mii -alloc r32 = ar.pfs,0,30,4,0 -(p0) add r33 = r0, r0 -(p0) cmp.eq.unc p0, p7 = r0, r0 ;; +{ .mlx + getf.exp GR_signexp_x = f8 // Get sign and exponent of x, redo if unorm + movl GR_sig_inv_ln2 = 0xb8aa3b295c17f0bc // significand of 1/ln2 } -exp_continue: -{ .mfi -(p0) add r32 = 2,r0 -(p0) fnorm.s1 f9 = f8 - nop.i 0 +{ .mlx + addl GR_ad_Arg = @ltoff(Constants_exp_64_Arg#),gp + movl GR_rshf_2to51 = 0x4718000000000000 // 1.10000 2^(63+51) } +;; + { .mfi -(p0) nop.m 0 + ld8 GR_ad_Arg = [GR_ad_Arg] // Point to Arg table + fclass.m p8, p0 = f8, 0x1E7 // Test x for natval, nan, inf, zero + cmp.eq p7, p6 = r0, r0 +} +{ .mfb + mov GR_exp_half = 0x0FFFE // Exponent of 0.5, for very small path + fnorm.s1 FR_norm_x = f8 // Normalize x + br.cond.sptk exp_continue +} +;; + +GLOBAL_IEEE754_END(expm1l) + + +GLOBAL_IEEE754_ENTRY(expl) // -// Set p7 false for exp -// Set Flag = r33 = 0 for exp +// Set p7 false for exp, p6 true // -(p0) fclass.m.unc p6, p8 = f8, 0x1E7 - nop.i 0;; +{ .mlx + getf.exp GR_signexp_x = f8 // Get sign and exponent of x, redo if unorm + movl GR_sig_inv_ln2 = 0xb8aa3b295c17f0bc // significand of 1/ln2 +} +{ .mlx + addl GR_ad_Arg = @ltoff(Constants_exp_64_Arg#),gp + movl GR_rshf_2to51 = 0x4718000000000000 // 1.10000 2^(63+51) } +;; + { .mfi - nop.m 999 -(p0) fclass.nm.unc p9, p0 = f8, 0x1FF - nop.i 0 + ld8 GR_ad_Arg = [GR_ad_Arg] // Point to Arg table + fclass.m p8, p0 = f8, 0x1E7 // Test x for natval, nan, inf, zero + cmp.eq p6, p7 = r0, r0 } { .mfi - nop.m 999 -(p0) mov f36 = f1 - nop.i 999 ;; + mov GR_exp_half = 0x0FFFE // Exponent of 0.5, for very small path + fnorm.s1 FR_norm_x = f8 // Normalize x + nop.i 999 } -{ .mfb - nop.m 999 -// -// Identify NatVals, NaNs, Infs, and Zeros. -// Identify EM unsupporteds. -// Save special input registers -(p0) mov f32 = f0 -// -// Create FR_X_cor = 0.0 -// GR_Flag = 0 -// GR_Expo_Range = 2 (r32) for double-extended precision -// FR_Scale = 1.0 -// -(p6) br.cond.spnt EXPL_64_SPECIAL ;; +;; + +exp_continue: +// Form two constants we need +// 1/ln2 * 2^63 to compute w = x * 1/ln2 * 128 +// 1.1000..000 * 2^(63+63-12) to right shift int(N) into the significand + +{ .mfi + setf.sig FR_INV_LN2_2TO63 = GR_sig_inv_ln2 // form 1/ln2 * 2^63 + fclass.nm.unc p9, p0 = f8, 0x1FF // Test x for unsupported + mov GR_exp_2tom51 = 0xffff-51 +} +{ .mlx + setf.d FR_RSHF_2TO51 = GR_rshf_2to51 // Form const 1.1000 * 2^(63+51) + movl GR_rshf = 0x43e8000000000000 // 1.10000 2^63 for right shift +} +;; + +{ .mfi + setf.exp FR_half = GR_exp_half // Form 0.5 for very small path + fma.s1 FR_scale = f1,f1,f0 // Scale = 1.0 + mov GR_exp_bias = 0x0FFFF // Set exponent bias } { .mib - nop.m 999 - nop.i 999 -(p9) br.cond.spnt EXPL_64_UNSUPPORTED ;; + add GR_ad_Limits = 0x20, GR_ad_Arg // Point to Limits table + mov GR_exp_mask = 0x1FFFF // Form exponent mask +(p8) br.cond.spnt EXP_64_SPECIAL // Branch if natval, nan, inf, zero } +;; + { .mfi -(p0) cmp.ne.unc p12, p13 = 0x01, r33 -// -// Branch out for special input values -// -(p0) fcmp.lt.unc.s0 p9,p0 = f8, f0 -(p0) cmp.eq.unc p15, p0 = r0, r0 + setf.exp FR_2TOM51 = GR_exp_2tom51 // Form 2^-51 for scaling float_N + nop.f 999 + add GR_ad_A = 0x40, GR_ad_Arg // Point to A table } -{ .mmi - nop.m 999 -// -// Raise possible denormal operand exception -// Normalize x -// -// This function computes expl( x + x_cor) -// Input FR 1: FR_X -// Input FR 2: FR_X_cor -// Input GR 1: GR_Flag -// Input GR 2: GR_Expo_Range -// Output FR 3: FR_Y_hi -// Output FR 4: FR_Y_lo -// Output FR 5: FR_Scale -// Output PR 1: PR_Safe -(p0) addl r34 = @ltoff(Constants_exp_64_Arg#),gp -(p0) addl r40 = @ltoff(Constants_exp_64_W1#),gp -};; -// -// Prepare to load constants -// Set Safe = True -// +{ .mib + setf.d FR_RSHF = GR_rshf // Form right shift const 1.1000 * 2^63 + add GR_ad_T1 = 0x160, GR_ad_Arg // Point to T1 table +(p9) br.cond.spnt EXP_64_UNSUPPORTED // Branch if unsupported +} +;; -{ .mmi - ld8 r34 = [r34] - ld8 r40 = [r40] -(p0) addl r41 = @ltoff(Constants_exp_64_W2#),gp +.pred.rel "mutex",p6,p7 +{ .mfi + ldfe FR_L_hi = [GR_ad_Arg],16 // Get L_hi + fcmp.eq.s0 p9,p0 = f8, f0 // Dummy op to flag denormals +(p6) add GR_ad_PQ = 0x30, GR_ad_A // Point to P table for exp +} +{ .mfi + ldfe FR_min_oflow_x = [GR_ad_Limits],16 // Get min x to cause overflow + fmpy.s1 FR_rsq = f8, f8 // rsq = x * x for small path +(p7) add GR_ad_PQ = 0x90, GR_ad_A // Point to Q table for expm1 };; { .mmi -(p0) ldfe f37 = [r34],16 -(p0) ld8 r41 = [r41] ;; + ldfe FR_L_lo = [GR_ad_Arg],16 // Get L_lo + ldfe FR_zero_uflow_x = [GR_ad_Limits],16 // Get x for zero uflow result + add GR_ad_W1 = 0x200, GR_ad_T1 // Point to W1 table } +;; -// -// N = fcvt.fx(float_N) -// Set p14 if -6 > expo_X -// -// -// Bias = 0x0FFFF -// expo_X = expo_X and Mask -// - -{ .mmi -(p0) ldfe f40 = [r34],16 - nop.m 999 -// -// Load L_lo -// Set p10 if 14 < expo_X -// -(p0) addl r50 = @ltoff(Constants_exp_64_T1#),gp +{ .mfi + ldfe FR_P6Q9 = [GR_ad_PQ],16 // P6(exp) or Q9(expm1) for small path + mov FR_r = FR_norm_x // r = X for small path + mov GR_very_small_exp = -60 // Exponent of x for very small path } -{ .mmi - nop.m 999 - nop.m 999 -(p0) addl r51 = @ltoff(Constants_exp_64_T2#),gp ;; +{ .mfi + add GR_ad_W2 = 0x400, GR_ad_T1 // Point to W2 table + nop.f 999 +(p7) mov GR_small_exp = -7 // Exponent of x for small path expm1 } -// -// Load W2_ptr -// Branch to SMALL is expo_X < -6 -// +;; -{.mmi -(p0) ld8 r50 = [r50] -(p0) ld8 r51 = [r51] -};; +{ .mmi + ldfe FR_P5Q8 = [GR_ad_PQ],16 // P5(exp) or Q8(expm1) for small path + and GR_exp_x = GR_signexp_x, GR_exp_mask +(p6) mov GR_small_exp = -12 // Exponent of x for small path exp +} +;; -{ .mlx -(p0) ldfe f41 = [r34],16 -// -// float_N = X * L_Inv -// expo_X = exponent of X -// Mask = 0x1FFFF -// -(p0) movl r58 = 0x0FFFF +// N_signif = X * Inv_log2_by_2^12 +// By adding 1.10...0*2^63 we shift and get round_int(N_signif) in significand. +// We actually add 1.10...0*2^51 to X * Inv_log2 to do the same thing. +{ .mfi + ldfe FR_P4Q7 = [GR_ad_PQ],16 // P4(exp) or Q7(expm1) for small path + fma.s1 FR_N_signif = FR_norm_x, FR_INV_LN2_2TO63, FR_RSHF_2TO51 + nop.i 999 } -{ .mlx - nop.m 999 -(p0) movl r39 = 0x1FFFF ;; +{ .mfi + sub GR_exp_x = GR_exp_x, GR_exp_bias // Get exponent + fmpy.s1 FR_r4 = FR_rsq, FR_rsq // Form r4 for small path + cmp.eq.unc p15, p0 = r0, r0 // Set Safe as default } +;; + { .mmi -(p0) getf.exp r37 = f9 - nop.m 999 -(p0) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp ;; + ldfe FR_P3Q6 = [GR_ad_PQ],16 // P3(exp) or Q6(expm1) for small path + cmp.lt p14, p0 = GR_exp_x, GR_very_small_exp // Is |x| < 2^-60? + nop.i 999 } -{ .mii -(p0) ld8 r34 = [r34] - nop.i 999 -(p0) and r37 = r37, r39 ;; +;; + +{ .mfi + ldfe FR_P2Q5 = [GR_ad_PQ],16 // P2(exp) or Q5(expm1) for small path + fmpy.s1 FR_half_x = FR_half, FR_norm_x // 0.5 * x for very small path + cmp.lt p13, p0 = GR_exp_x, GR_small_exp // Is |x| < 2^-m? } -{ .mmi -(p0) sub r37 = r37, r58 ;; -(p0) cmp.gt.unc p14, p0 = -6, r37 -(p0) cmp.lt.unc p10, p0 = 14, r37 ;; +{ .mib + nop.m 999 + nop.i 999 +(p14) br.cond.spnt EXP_VERY_SMALL // Branch if |x| < 2^-60 } +;; + { .mfi -(p0) nop.m 0 -// -// Load L_inv -// Set p12 true for Flag = 0 (exp) -// Set p13 true for Flag = 1 (expm1) -// -(p0) fmpy.s1 f38 = f9, f37 - nop.i 999 ;; + ldfe FR_A3 = [GR_ad_A],16 // Get A3 for normal path + fcmp.ge.s1 p10,p0 = FR_norm_x, FR_min_oflow_x // Will result overflow? + mov GR_big_expo_neg = -16381 // -0x3ffd } { .mfb - nop.m 999 -// -// Load L_hi -// expo_X = expo_X - Bias -// get W1_ptr -// -(p0) fcvt.fx.s1 f39 = f38 -(p14) br.cond.spnt EXPL_SMALL ;; + ldfe FR_P1Q4 = [GR_ad_PQ],16 // P1(exp) or Q4(expm1) for small path + nop.f 999 +(p13) br.cond.spnt EXP_SMALL // Branch if |x| < 2^-m + // m=12 for exp, m=7 for expm1 } -{ .mib - nop.m 999 - nop.i 999 -(p10) br.cond.spnt EXPL_HUGE ;; +;; + +// Now we are on the main path for |x| >= 2^-m, m=12 for exp, m=7 for expm1 +// +// float_N = round_int(N_signif) +// The signficand of N_signif contains the rounded integer part of X * 2^12/ln2, +// as a twos complement number in the lower bits (that is, it may be negative). +// That twos complement number (called N) is put into GR_N. + +// Since N_signif is scaled by 2^51, it must be multiplied by 2^-51 +// before the shift constant 1.10000 * 2^63 is subtracted to yield float_N. +// Thus, float_N contains the floating point version of N + + +{ .mfi + ldfe FR_A2 = [GR_ad_A],16 // Get A2 for main path + fcmp.lt.s1 p11,p0 = FR_norm_x, FR_zero_uflow_x // Certain zero, uflow? + add GR_ad_T2 = 0x100, GR_ad_T1 // Point to T2 table } -{ .mmi -(p0) shladd r34 = r32,4,r34 +{ .mfi nop.m 999 -(p0) addl r35 = @ltoff(Constants_exp_64_A#),gp ;; -} -// -// Load T_1,T_2 -// -{ .mmi - nop.m 999 - ld8 r35 =[r35] - nop.i 99 -};; -{ .mmb -(p0) ldfe f51 = [r35],16 -(p0) ld8 r45 = [r34],8 - nop.b 999 ;; + fms.s1 FR_float_N = FR_N_signif, FR_2TOM51, FR_RSHF // Form float_N + nop.i 999 } -// -// Set Safe = True if k >= big_expo_neg -// Set Safe = False if k < big_expo_neg -// -{ .mmb -(p0) ldfe f49 = [r35],16 -(p0) ld8 r48 = [r34],0 - nop.b 999 ;; +;; + +{ .mbb + getf.sig GR_N_fix = FR_N_signif // Get N from significand +(p10) br.cond.spnt EXP_OVERFLOW // Branch if result will overflow +(p11) br.cond.spnt EXP_CERTAIN_UNDERFLOW_ZERO // Branch if certain zero, uflow } +;; + { .mfi - nop.m 999 -// -// Branch to HUGE is expo_X > 14 -// -(p0) fcvt.xf f38 = f39 - nop.i 999 ;; + ldfe FR_A1 = [GR_ad_A],16 // Get A1 for main path + fnma.s1 FR_r = FR_L_hi, FR_float_N, FR_norm_x // r = -L_hi * float_N + x + extr.u GR_M1 = GR_N_fix, 6, 6 // Extract index M_1 } { .mfi -(p0) getf.sig r52 = f39 - nop.f 999 - nop.i 999 ;; -} -{ .mii - nop.m 999 -(p0) extr.u r43 = r52, 6, 6 ;; -// -// r = r - float_N * L_lo -// K = extr(N_fix,12,52) -// -(p0) shladd r40 = r43,3,r40 ;; + and GR_M2 = 0x3f, GR_N_fix // Extract index M_2 + nop.f 999 + nop.i 999 } +;; + +// N_fix is only correct up to 50 bits because of our right shift technique. +// Actually in the normal path we will have restricted K to about 14 bits. +// Somewhat arbitrarily we extract 32 bits. { .mfi -(p0) shladd r50 = r43,2,r50 -(p0) fnma.s1 f42 = f40, f38, f9 -// -// float_N = float(N) -// N_fix = signficand N -// -(p0) extr.u r42 = r52, 0, 6 + shladd GR_ad_W1 = GR_M1,3,GR_ad_W1 // Point to W1 + nop.f 999 + extr GR_K = GR_N_fix, 12, 32 // Extract limited range K } -{ .mmi -(p0) ldfd f43 = [r40],0 ;; -(p0) shladd r41 = r42,3,r41 -(p0) shladd r51 = r42,2,r51 +{ .mfi + shladd GR_ad_T1 = GR_M1,2,GR_ad_T1 // Point to T1 + nop.f 999 + shladd GR_ad_T2 = GR_M2,2,GR_ad_T2 // Point to T2 } -// -// W_1_p1 = 1 + W_1 -// +;; + { .mmi -(p0) ldfs f44 = [r50],0 ;; -(p0) ldfd f45 = [r41],0 -// -// M_2 = extr(N_fix,0,6) -// M_1 = extr(N_fix,6,6) -// r = X - float_N * L_hi -// -(p0) extr r44 = r52, 12, 52 + ldfs FR_T1 = [GR_ad_T1],0 // Get T1 + ldfd FR_W1 = [GR_ad_W1],0 // Get W1 + add GR_exp_2_k = GR_exp_bias, GR_K // Form exponent of 2^k } +;; + { .mmi -(p0) ldfs f46 = [r51],0 ;; -(p0) sub r46 = r58, r44 -(p0) cmp.gt.unc p8, p15 = r44, r45 -} -// -// W = W_1 + W_1_p1*W_2 -// Load A_2 -// Bias_m_K = Bias - K -// -{ .mii -(p0) ldfe f40 = [r35],16 -// -// load A_1 -// poly = A_2 + r*A_3 -// rsq = r * r -// neg_2_mK = exponent of Bias_m_k -// -(p0) add r47 = r58, r44 ;; -// -// Set Safe = True if k <= big_expo_pos -// Set Safe = False if k > big_expo_pos -// Load A_3 -// -(p15) cmp.lt p8,p15 = r44,r48 ;; + ldfs FR_T2 = [GR_ad_T2],0 // Get T2 + shladd GR_ad_W2 = GR_M2,3,GR_ad_W2 // Point to W2 + sub GR_exp_2_mk = GR_exp_bias, GR_K // Form exponent of 2^-k } +;; + { .mmf -(p0) setf.exp f61 = r46 -// -// Bias_p + K = Bias + K -// T = T_1 * T_2 -// -(p0) setf.exp f36 = r47 -(p0) fnma.s1 f42 = f41, f38, f42 ;; + ldfd FR_W2 = [GR_ad_W2],0 // Get W2 + setf.exp FR_scale = GR_exp_2_k // Set scale = 2^k + fnma.s1 FR_r = FR_L_lo, FR_float_N, FR_r // r = -L_lo * float_N + r } +;; + { .mfi - nop.m 999 -// -// Load W_1,W_2 -// Load big_exp_pos, load big_exp_neg -// -(p0) fadd.s1 f47 = f43, f1 - nop.i 999 ;; + setf.exp FR_2_mk = GR_exp_2_mk // Form 2^-k + fma.s1 FR_poly = FR_r, FR_A3, FR_A2 // poly = r * A3 + A2 + cmp.lt p8,p15 = GR_K,GR_big_expo_neg // Set Safe if K > big_expo_neg } { .mfi - nop.m 999 -(p0) fma.s1 f52 = f42, f51, f49 - nop.i 999 + nop.m 999 + fmpy.s1 FR_rsq = FR_r, FR_r // rsq = r * r + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fmpy.s1 f48 = f42, f42 - nop.i 999 ;; + nop.m 999 + fmpy.s1 FR_T = FR_T1, FR_T2 // T = T1 * T2 + nop.i 999 } { .mfi - nop.m 999 -(p0) fmpy.s1 f53 = f44, f46 - nop.i 999 ;; + nop.m 999 + fadd.s1 FR_W1_p1 = FR_W1, f1 // W1_p1 = W1 + 1.0 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fma.s1 f54 = f45, f47, f43 - nop.i 999 +(p7) cmp.lt.unc p8, p9 = 10, GR_K // If expm1, set p8 if K > 10 + fma.s1 FR_poly = FR_r, FR_poly, FR_A1 // poly = r * poly + A1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fneg f61 = f61 - nop.i 999 ;; +(p7) cmp.eq p15, p0 = r0, r0 // If expm1, set Safe flag + fma.s1 FR_T_scale = FR_T, FR_scale, f0 // T_scale = T * scale +(p9) cmp.gt.unc p9, p10 = -10, GR_K // If expm1, set p9 if K < -10 + // If expm1, set p10 if -10<=K<=10 } { .mfi - nop.m 999 -(p0) fma.s1 f52 = f42, f52, f40 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_W = FR_W2, FR_W1_p1, FR_W1 // W = W2 * (W1+1.0) + W1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fadd.s1 f55 = f54, f1 - nop.i 999 + nop.m 999 + mov FR_Y_hi = FR_T // Assume Y_hi = T + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// W + Wp1 * poly -// -(p0) mov f34 = f53 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_poly = FR_rsq, FR_poly, FR_r // poly = rsq * poly + r + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// A_1 + r * poly -// Scale = setf_expl(Bias_p_k) -// -(p0) fma.s1 f52 = f48, f52, f42 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_Wp1_T_scale = FR_W, FR_T_scale, FR_T_scale // (W+1)*T*scale + nop.i 999 } { .mfi - nop.m 999 -// -// poly = r + rsq(A_1 + r*poly) -// Wp1 = 1 + W -// neg_2_mK = -neg_2_mK -// -(p0) fma.s1 f35 = f55, f52, f54 - nop.i 999 ;; -} -{ .mfb - nop.m 999 -(p0) fmpy.s1 f35 = f35, f53 -// -// Y_hi = T -// Y_lo = T * (W + Wp1*poly) -// -(p12) br.cond.sptk EXPL_MAIN ;; -} -// -// Branch if expl(x) -// Continue for expl(x-1) -// -{ .mii -(p0) cmp.lt.unc p12, p13 = 10, r44 - nop.i 999 ;; -// -// Set p12 if 10 < K, Else p13 -// -(p13) cmp.gt.unc p13, p14 = -10, r44 ;; + nop.m 999 + fma.s1 FR_W_T_scale = FR_W, FR_T_scale, f0 // W*T*scale + nop.i 999 } -// -// K > 10: Y_lo = Y_lo + neg_2_mK -// K <=10: Set p13 if -10 > K, Else set p14 -// +;; + { .mfi -(p13) cmp.eq p15, p0 = r0, r0 -(p14) fadd.s1 f34 = f61, f34 - nop.i 999 ;; + nop.m 999 +(p9) fsub.s1 FR_Y_hi = f0, FR_2_mk // If expm1, if K < -10 set Y_hi + nop.i 999 } { .mfi - nop.m 999 -(p12) fadd.s1 f35 = f35, f61 - nop.i 999 ;; + nop.m 999 +(p10) fsub.s1 FR_Y_hi = FR_T, FR_2_mk // If expm1, if |K|<=10 set Y_hi + nop.i 999 } +;; + { .mfi - nop.m 999 -(p13) fadd.s1 f35 = f35, f34 - nop.i 999 -} -{ .mfb - nop.m 999 -// -// K <= 10 and K < -10, Set Safe = True -// K <= 10 and K < 10, Y_lo = Y_hi + Y_lo -// K <= 10 and K > =-10, Y_hi = Y_hi + neg_2_mk -// -(p13) mov f34 = f61 -(p0) br.cond.sptk EXPL_MAIN ;; -} -EXPL_SMALL: -{ .mmi nop.m 999 -(p0) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp -(p12) addl r35 = @ltoff(Constants_exp_64_P#),gp ;; + fma.s1 FR_result_lo = FR_Wp1_T_scale, FR_poly, FR_W_T_scale + nop.i 999 } -.pred.rel "mutex",p12,p13 -{ .mmi -(p12) ld8 r35=[r35] -nop.m 999 -(p13) addl r35 = @ltoff(Constants_exp_64_Q#),gp -};; -{ .mmi -(p13) ld8 r35=[r35] -(p0) ld8 r34=[r34] -nop.i 999 -};; +;; + +.pred.rel "mutex",p8,p9 +// If K > 10 adjust result_lo = result_lo - scale * 2^-k +// If |K| <= 10 adjust result_lo = result_lo + scale * T { .mfi -(p0) add r34 = 0x48,r34 -// -// Return -// K <= 10 and K < 10, Y_hi = neg_2_mk -// -// /*******************************************************/ -// /*********** Branch EXPL_SMALL ************************/ -// /*******************************************************/ -(p0) mov f42 = f9 - nop.i 999 ;; + nop.m 999 +(p8) fnma.s1 FR_result_lo = FR_scale, FR_2_mk, FR_result_lo // If K > 10 + nop.i 999 } -// -// Flag = 0 -// r4 = rsq * rsq -// { .mfi -(p0) ld8 r49 =[r34],0 - nop.f 999 - nop.i 999 ;; -} -{ .mii - nop.m 999 - nop.i 999 ;; -// -// Flag = 1 -// -(p0) cmp.lt.unc p14, p0 = r37, r49 ;; + nop.m 999 +(p9) fma.s1 FR_result_lo = FR_T_scale, f1, FR_result_lo // If |K| <= 10 + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// r = X -// -(p0) fmpy.s1 f48 = f42, f42 - nop.i 999 ;; + nop.m 999 + fmpy.s0 FR_tmp = FR_A1, FR_A1 // Dummy op to set inexact + nop.i 999 } { .mfb - nop.m 999 -// -// rsq = r * r -// -(p0) fmpy.s1 f50 = f48, f48 -// -// Is input very small? -// -(p14) br.cond.spnt EXPL_VERY_SMALL ;; -} -// -// Flag_not1: Y_hi = 1.0 -// Flag is 1: r6 = rsq * r4 -// -{ .mfi -(p12) ldfe f52 = [r35],16 -(p12) mov f34 = f1 -(p0) add r53 = 0x1,r0 ;; -} -{ .mfi -(p13) ldfe f51 = [r35],16 -// -// Flag_not_1: Y_lo = poly_hi + r4 * poly_lo -// -(p13) mov f34 = f9 - nop.i 999 ;; -} -{ .mmf -(p12) ldfe f53 = [r35],16 -// -// For Flag_not_1, Y_hi = X -// Scale = 1 -// Create 0x000...01 -// -(p0) setf.sig f37 = r53 -(p0) mov f36 = f1 ;; -} -{ .mmi -(p13) ldfe f52 = [r35],16 ;; -(p12) ldfe f54 = [r35],16 - nop.i 999 ;; -} -{ .mfi -(p13) ldfe f53 = [r35],16 -(p13) fmpy.s1 f58 = f48, f50 - nop.i 999 ;; -} -// -// Flag_not1: poly_lo = P_5 + r*P_6 -// Flag_1: poly_lo = Q_6 + r*Q_7 -// -{ .mmi -(p13) ldfe f54 = [r35],16 ;; -(p12) ldfe f55 = [r35],16 - nop.i 999 ;; -} -{ .mmi -(p12) ldfe f56 = [r35],16 ;; -(p13) ldfe f55 = [r35],16 - nop.i 999 ;; + nop.m 999 +(p15) fma.s0 f8 = FR_Y_hi, FR_scale, FR_result_lo // Safe result +(p15) br.ret.sptk b0 // Safe exit for normal path } -{ .mmi -(p12) ldfe f57 = [r35],0 ;; -(p13) ldfe f56 = [r35],16 - nop.i 999 ;; +;; + +// Here if unsafe, will only be here for exp with K < big_expo_neg +{ .mfb + nop.m 999 + fma.s0 FR_RESULT = FR_Y_hi, FR_scale, FR_result_lo // Prelim result + br.cond.sptk EXP_POSSIBLE_UNDERFLOW // Branch to unsafe code } +;; + + +EXP_SMALL: +// Here if 2^-60 < |x| < 2^-m, m=12 for exp, m=7 for expm1 { .mfi -(p13) ldfe f57 = [r35],0 - nop.f 999 - nop.i 999 ;; +(p7) ldfe FR_Q3 = [GR_ad_Q],16 // Get Q3 for small path, if expm1 +(p6) fma.s1 FR_p65 = FR_P6, FR_r, FR_P5 // If exp, p65 = P6 * r + P5 + nop.i 999 } { .mfi - nop.m 999 -// -// For Flag_not_1, load p5,p6,p1,p2 -// Else load p5,p6,p1,p2 -// -(p12) fma.s1 f60 = f52, f42, f53 - nop.i 999 ;; + mov GR_minus_one = -1 +(p7) fma.s1 FR_q98 = FR_Q9, FR_r, FR_Q8 // If expm1, q98 = Q9 * r + Q8 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p13) fma.s1 f60 = f51, f42, f52 - nop.i 999 ;; +(p7) ldfe FR_Q2 = [GR_ad_Q],16 // Get Q2 for small path, if expm1 +(p7) fma.s1 FR_q65 = FR_Q6, FR_r, FR_Q5 // If expm1, q65 = Q6 * r + Q5 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p12) fma.s1 f60 = f60, f42, f54 - nop.i 999 ;; + setf.sig FR_tmp = GR_minus_one // Create value to force inexact +(p6) fma.s1 FR_p21 = FR_P2, FR_r, FR_P1 // If exp, p21 = P2 * r + P1 + nop.i 999 } { .mfi - nop.m 999 -(p12) fma.s1 f59 = f56, f42, f57 - nop.i 999 ;; +(p7) ldfe FR_Q1 = [GR_ad_Q],16 // Get Q1 for small path, if expm1 +(p7) fma.s1 FR_q43 = FR_Q4, FR_r, FR_Q3 // If expm1, q43 = Q4 * r + Q3 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p13) fma.s1 f60 = f42, f60, f53 - nop.i 999 ;; + nop.m 999 +(p6) fma.s1 FR_p654 = FR_p65, FR_r, FR_P4 // If exp, p654 = p65 * r + P4 + nop.i 999 } { .mfi - nop.m 999 -(p12) fma.s1 f59 = f59, f48, f42 - nop.i 999 ;; + nop.m 999 +(p7) fma.s1 FR_q987 = FR_q98, FR_r, FR_Q7 // If expm1, q987 = q98 * r + Q7 + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Flag_1: poly_lo = Q_5 + r*(Q_6 + r*Q_7) -// Flag_not1: poly_lo = P_4 + r*(P_5 + r*P_6) -// Flag_not1: poly_hi = (P_1 + r*P_2) -// -(p13) fmpy.s1 f60 = f60, f58 - nop.i 999 ;; + nop.m 999 +(p7) fma.s1 FR_q21 = FR_Q2, FR_r, FR_Q1 // If expm1, q21 = Q2 * r + Q1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p12) fma.s1 f60 = f60, f42, f55 - nop.i 999 ;; + nop.m 999 +(p6) fma.s1 FR_p210 = FR_p21, FR_rsq, FR_r // If exp, p210 = p21 * r + P0 + nop.i 999 } { .mfi - nop.m 999 -// -// Flag_1: poly_lo = r6 *(Q_5 + ....) -// Flag_not1: poly_hi = r + rsq *(P_1 + r*P_2) -// -(p12) fma.s1 f35 = f60, f50, f59 - nop.i 999 + nop.m 999 +(p7) fma.s1 FR_q6543 = FR_q65, FR_rsq, FR_q43 // If expm1, q6543 = q65*r2+q43 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p13) fma.s1 f59 = f54, f42, f55 - nop.i 999 ;; + nop.m 999 +(p6) fma.s1 FR_p6543 = FR_p654, FR_r, FR_P3 // If exp, p6543 = p654 * r + P3 + nop.i 999 } { .mfi - nop.m 999 -// -// Flag_not1: Y_lo = rsq* poly_hi + poly_lo -// Flag_1: poly_lo = rsq* poly_hi + poly_lo -// -(p13) fma.s1 f59 = f59, f42, f56 - nop.i 999 ;; + nop.m 999 +(p7) fma.s1 FR_q9876543 = FR_q987, FR_r4, FR_q6543 // If expm1, q9876543 = ... + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Flag_not_1: (P_1 + r*P_2) -// -(p13) fma.s1 f59 = f59, f42, f57 - nop.i 999 ;; + nop.m 999 +(p6) fma.s1 FR_Y_lo = FR_p6543, FR_r4, FR_p210 // If exp, form Y_lo + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Flag_not_1: poly_hi = r + rsq * (P_1 + r*P_2) -// -(p13) fma.s1 f35 = f59, f48, f60 - nop.i 999 ;; + nop.m 999 +(p7) fma.s1 FR_Y_lo = FR_q9876543, FR_rsq, FR_q21 // If expm1, form Y_lo + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Create 0.000...01 -// -(p0) for f37 = f35, f37 - nop.i 999 ;; -} -{ .mfb - nop.m 999 -// -// Set lsb of Y_lo to 1 -// -(p0) fmerge.se f35 = f35,f37 -(p0) br.cond.sptk EXPL_MAIN ;; -} -EXPL_VERY_SMALL: -{ .mmi - nop.m 999 - nop.m 999 -(p13) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp + nop.m 999 + fmpy.s0 FR_tmp = FR_tmp, FR_tmp // Dummy op to set inexact + nop.i 999 } +;; + +.pred.rel "mutex",p6,p7 { .mfi - nop.m 999 -(p12) mov f35 = f9 - nop.i 999 ;; + nop.m 999 +(p6) fma.s0 f8 = FR_Y_lo, f1, f1 // If exp, result = 1 + Y_lo + nop.i 999 } { .mfb -(p13) ld8 r34 = [r34] -(p12) mov f34 = f1 -(p12) br.cond.sptk EXPL_MAIN ;; -} -{ .mlx -(p13) add r34 = 8,r34 -(p13) movl r39 = 0x0FFFE ;; -} -// -// Load big_exp_neg -// Create 1/2's exponent -// -{ .mii -(p13) setf.exp f56 = r39 -(p13) shladd r34 = r32,4,r34 ;; - nop.i 999 + nop.m 999 +(p7) fma.s0 f8 = FR_Y_lo, FR_rsq, FR_norm_x // If expm1, result = Y_lo*r2+x + br.ret.sptk b0 // Exit for 2^-60 <= |x| < 2^-m + // m=12 for exp, m=7 for expm1 } +;; + + +EXP_VERY_SMALL: // -// Negative exponents are stored after positive +// Here if 0 < |x| < 2^-60 +// If exp, result = 1.0 + x +// If expm1, result = x +x*x/2, but have to check for possible underflow // + { .mfi -(p13) ld8 r45 = [r34],0 -// -// Y_hi = x -// Scale = 1 -// -(p13) fmpy.s1 f35 = f9, f9 - nop.i 999 ;; +(p7) mov GR_exp_underflow = -16381 // Exponent for possible underflow +(p6) fadd.s0 f8 = f1, FR_norm_x // If exp, result = 1+x + nop.i 999 } { .mfi - nop.m 999 -// -// Reset Safe if necessary -// Create 1/2 -// -(p13) mov f34 = f9 - nop.i 999 ;; + nop.m 999 +(p7) fmpy.s1 FR_result_lo = FR_half_x, FR_norm_x // If expm1 result_lo = x*x/2 + nop.i 999 } +;; + { .mfi -(p13) cmp.lt.unc p0, p15 = r37, r45 -(p13) mov f36 = f1 - nop.i 999 ;; +(p7) cmp.lt.unc p0, p8 = GR_exp_x, GR_exp_underflow // Unsafe if expm1 x small +(p7) mov FR_Y_hi = FR_norm_x // If expm1, Y_hi = x +(p7) cmp.lt p0, p15 = GR_exp_x, GR_exp_underflow // Unsafe if expm1 x small } +;; + { .mfb - nop.m 999 -// -// Y_lo = x * x -// -(p13) fmpy.s1 f35 = f35, f56 -// -// Y_lo = x*x/2 -// -(p13) br.cond.sptk EXPL_MAIN ;; -} -EXPL_HUGE: -{ .mfi - nop.m 999 -(p0) fcmp.gt.unc.s1 p14, p0 = f9, f0 - nop.i 999 -} -{ .mlx - nop.m 999 -(p0) movl r39 = 0x15DC0 ;; -} -{ .mfi -(p14) setf.exp f34 = r39 -(p14) mov f35 = f1 -(p14) cmp.eq p0, p15 = r0, r0 ;; + nop.m 999 +(p8) fma.s0 f8 = FR_norm_x, f1, FR_result_lo // If expm1, result=x+x*x/2 +(p15) br.ret.sptk b0 // If Safe, exit } +;; + +// Here if expm1 and 0 < |x| < 2^-16381; may be possible underflow { .mfb - nop.m 999 -(p14) mov f36 = f34 -// -// If x > 0, Set Safe = False -// If x > 0, Y_hi = 2**(24,000) -// If x > 0, Y_lo = 1.0 -// If x > 0, Scale = 2**(24,000) -// -(p14) br.cond.sptk EXPL_MAIN ;; -} -{ .mlx - nop.m 999 -(p12) movl r39 = 0xA240 -} -{ .mlx - nop.m 999 -(p12) movl r38 = 0xA1DC ;; -} -{ .mmb -(p13) cmp.eq p15, p14 = r0, r0 -(p12) setf.exp f34 = r39 - nop.b 999 ;; -} -{ .mlx -(p12) setf.exp f35 = r38 -(p13) movl r39 = 0xFF9C -} -{ .mfi - nop.m 999 -(p13) fsub.s1 f34 = f0, f1 - nop.i 999 ;; + nop.m 999 + fma.s0 FR_RESULT = FR_Y_hi, FR_scale, FR_result_lo // Prelim result + br.cond.sptk EXP_POSSIBLE_UNDERFLOW // Branch to unsafe code } -{ .mfi - nop.m 999 -(p12) mov f36 = f34 -(p12) cmp.eq p0, p15 = r0, r0 ;; +;; + +EXP_CERTAIN_UNDERFLOW_ZERO: +// Here if x < zero_uflow_x +// For exp, set result to tiny+0.0 and set I, U, and branch to error handling +// For expm1, set result to tiny-1.0 and set I, and exit +{ .mmi + alloc GR_SAVE_PFS = ar.pfs,0,3,4,0 + nop.m 999 + mov GR_one = 1 } -{ .mfi -(p13) setf.exp f35 = r39 -(p13) mov f36 = f1 - nop.i 999 ;; +;; + +{ .mmi + setf.exp FR_small = GR_one // Form small value + nop.m 999 +(p6) mov GR_Parameter_TAG = 13 // Error tag for exp underflow } -EXPL_MAIN: +;; + { .mfi -(p0) cmp.ne.unc p12, p0 = 0x01, r33 -(p0) fmpy.s1 f101 = f36, f35 - nop.i 999 ;; + nop.m 999 + fmerge.s FR_X = f8,f8 // Save x for error call + nop.i 999 } +;; + +.pred.rel "mutex",p6,p7 { .mfb - nop.m 999 -(p0) fma.s0 f99 = f34, f36, f101 -(p15) br.cond.sptk EXPL_64_RETURN ;; -} -{ .mfi - nop.m 999 -(p0) fsetc.s3 0x7F,0x01 - nop.i 999 -} -{ .mlx - nop.m 999 -(p0) movl r50 = 0x00000000013FFF ;; + nop.m 999 +(p6) fma.s0 FR_RESULT = FR_small, FR_small, f0 // If exp, set I,U, tiny result +(p6) br.cond.sptk __libm_error_region // If exp, go to error handling } -// -// S0 user supplied status -// S2 user supplied status + WRE + TD (Overflows) -// S3 user supplied status + RZ + TD (Underflows) -// -// -// If (Safe) is true, then -// Compute result using user supplied status field. -// No overflow or underflow here, but perhaps inexact. -// Return -// Else -// Determine if overflow or underflow was raised. -// Fetch +/- overflow threshold for IEEE single, double, -// double extended -// -{ .mfi -(p0) setf.exp f60 = r50 -(p0) fma.s3 f102 = f34, f36, f101 - nop.i 999 +{ .mfb + nop.m 999 +(p7) fms.s0 f8 = FR_small, FR_small, f1 // If expm1, set I, result -1.0 +(p7) br.ret.sptk b0 // If expm1, exit +} +;; + + +EXP_OVERFLOW: +// Here if x >= min_oflow_x +{ .mmi + alloc GR_SAVE_PFS = ar.pfs,0,3,4,0 + mov GR_huge_exp = 0x1fffe + nop.i 999 } { .mfi - nop.m 999 -(p0) fsetc.s3 0x7F,0x40 - nop.i 999 ;; + mov GR_huge_signif = -0x1 + nop.f 999 +(p6) mov GR_Parameter_TAG = 12 // Error tag for exp overflow } -{ .mfi - nop.m 999 -// -// For Safe, no need to check for over/under. -// For expm1, handle errors like exp. -// -(p0) fsetc.s2 0x7F,0x42 - nop.i 999;; +;; + +{ .mmf + setf.exp FR_huge_exp = GR_huge_exp // Create huge value + setf.sig FR_huge_signif = GR_huge_signif // Create huge value + fmerge.s FR_X = f8,f8 // Save x for error call } +;; + { .mfi - nop.m 999 -(p0) fma.s2 f100 = f34, f36, f101 - nop.i 999 ;; + nop.m 999 + fmerge.se FR_huge = FR_huge_exp, FR_huge_signif +(p7) mov GR_Parameter_TAG = 39 // Error tag for expm1 overflow } -{ .mfi - nop.m 999 -(p0) fsetc.s2 0x7F,0x40 - nop.i 999 ;; +;; + +{ .mfb + nop.m 999 + fma.s0 FR_RESULT = FR_huge, FR_huge, FR_huge // Force I, O, and Inf + br.cond.sptk __libm_error_region // Branch to error handling } +;; + + + +EXP_POSSIBLE_UNDERFLOW: +// Here if exp and zero_uflow_x < x < about -11356 [where k < -16381] +// Here if expm1 and |x| < 2^-16381 { .mfi - nop.m 999 -(p7) fclass.m.unc p12, p0 = f102, 0x00F - nop.i 999 + alloc GR_SAVE_PFS = ar.pfs,0,3,4,0 + fsetc.s2 0x7F,0x41 // Set FTZ and disable traps + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fclass.m.unc p11, p0 = f102, 0x00F - nop.i 999 ;; + nop.m 999 + fma.s2 FR_ftz = FR_Y_hi, FR_scale, FR_result_lo // Result with FTZ + nop.i 999 } +;; + { .mfi - nop.m 999 -(p7) fcmp.ge.unc.s1 p10, p0 = f100, f60 - nop.i 999 + nop.m 999 + fsetc.s2 0x7F,0x40 // Disable traps (set s2 default) + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Create largest double exponent + 1. -// Create smallest double exponent - 1. -// -(p0) fcmp.ge.unc.s1 p8, p0 = f100, f60 - nop.i 999 ;; -} -// -// fcmp: resultS2 >= + overflow threshold -> set (a) if true -// fcmp: resultS2 <= - overflow threshold -> set (b) if true -// fclass: resultS3 is denorm/unorm/0 -> set (d) if true -// -{ .mib -(p10) mov GR_Parameter_TAG = 39 - nop.i 999 -(p10) br.cond.sptk __libm_error_region ;; -} -{ .mib -(p8) mov GR_Parameter_TAG = 12 - nop.i 999 -(p8) br.cond.sptk __libm_error_region ;; -} -// -// Report that exp overflowed -// -{ .mib -(p12) mov GR_Parameter_TAG = 40 - nop.i 999 -(p12) br.cond.sptk __libm_error_region ;; -} -{ .mib -(p11) mov GR_Parameter_TAG = 13 - nop.i 999 -(p11) br.cond.sptk __libm_error_region ;; -} -{ .mib - nop.m 999 - nop.i 999 -// -// Report that exp underflowed -// -(p0) br.cond.sptk EXPL_64_RETURN ;; + nop.m 999 +(p6) fclass.m.unc p11, p0 = FR_ftz, 0x00F // If exp, FTZ result denorm or zero? + nop.i 999 } -EXPL_64_SPECIAL: -{ .mfi - nop.m 999 -(p0) fclass.m.unc p6, p0 = f8, 0x0c3 - nop.i 999 +;; + +{ .mfb +(p11) mov GR_Parameter_TAG = 13 // exp underflow + fmerge.s FR_X = f8,f8 // Save x for error call +(p11) br.cond.spnt __libm_error_region // Branch on exp underflow } -{ .mfi - nop.m 999 -(p0) fclass.m.unc p13, p8 = f8, 0x007 - nop.i 999 ;; +;; + +{ .mfb + nop.m 999 + mov f8 = FR_RESULT // Was safe after all + br.ret.sptk b0 } +;; + + +EXP_64_SPECIAL: +// Here if x natval, nan, inf, zero +// If x natval, +inf, or if expm1 and x zero, just return x. +// The other cases must be tested for, and results set. +// These cases do not generate exceptions. { .mfi - nop.m 999 -(p7) fclass.m.unc p14, p0 = f8, 0x007 - nop.i 999 + nop.m 999 + fclass.m p8, p0 = f8, 0x0c3 // Is x nan? + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fclass.m.unc p12, p9 = f8, 0x021 - nop.i 999 ;; + nop.m 999 +(p6) fclass.m.unc p13, p0 = f8, 0x007 // If exp, is x zero? + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fclass.m.unc p11, p0 = f8, 0x022 - nop.i 999 + nop.m 999 +(p6) fclass.m.unc p11, p0 = f8, 0x022 // If exp, is x -inf? + nop.i 999 } { .mfi - nop.m 999 -(p7) fclass.m.unc p10, p0 = f8, 0x022 - nop.i 999 ;; + nop.m 999 +(p8) fadd.s0 f8 = f8, f1 // If x nan, result quietized x + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// Identify +/- 0, Inf, or -Inf -// Generate the right kind of NaN. -// -(p13) fadd.s0 f99 = f0, f1 - nop.i 999 ;; + nop.m 999 +(p7) fclass.m.unc p10, p0 = f8, 0x022 // If expm1, is x -inf? + nop.i 999 } { .mfi - nop.m 999 -(p14) mov f99 = f8 - nop.i 999 ;; -} -{ .mfb - nop.m 999 -(p6) fadd.s0 f99 = f8, f1 -// -// expl(+/-0) = 1 -// expm1l(+/-0) = +/-0 -// No exceptions raised -// -(p6) br.cond.sptk EXPL_64_RETURN ;; -} -{ .mib - nop.m 999 - nop.i 999 -(p14) br.cond.sptk EXPL_64_RETURN ;; + nop.m 999 +(p13) fadd.s0 f8 = f0, f1 // If exp and x zero, result 1.0 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p11) mov f99 = f0 - nop.i 999 ;; -} -{ .mfb - nop.m 999 -(p10) fsub.s1 f99 = f0, f1 -// -// expl(-Inf) = 0 -// expm1l(-Inf) = -1 -// No exceptions raised. -// -(p10) br.cond.sptk EXPL_64_RETURN ;; -} -{ .mfb - nop.m 999 -(p12) fmpy.s1 f99 = f8, f1 -// -// expl(+Inf) = Inf -// No exceptions raised. -// -(p0) br.cond.sptk EXPL_64_RETURN ;; + nop.m 999 +(p11) mov f8 = f0 // If exp and x -inf, result 0 + nop.i 999 } -EXPL_64_UNSUPPORTED: +;; + { .mfb - nop.m 999 -(p0) fmpy.s0 f99 = f8, f0 -(p0) br.cond.sptk EXPL_64_RETURN ;; + nop.m 999 +(p10) fsub.s1 f8 = f0, f1 // If expm1, x -inf, result -1.0 + br.ret.sptk b0 // Exit special cases } -EXPL_64_RETURN: +;; + + +EXP_64_UNSUPPORTED: +// Here if x unsupported type { .mfb nop.m 999 -(p0) mov f8 = f99 -(p0) br.ret.sptk b0 + fmpy.s0 f8 = f8, f0 // Return nan + br.ret.sptk b0 } -.endp -ASM_SIZE_DIRECTIVE(expl) +;; + +GLOBAL_IEEE754_END(expl) -.proc __libm_error_region -__libm_error_region: +LOCAL_LIBM_ENTRY(__libm_error_region) .prologue { .mfi add GR_Parameter_Y=-32,sp // Parameter 2 value @@ -1598,9 +1410,9 @@ __libm_error_region: br.call.sptk b0=__libm_error_support# // Call error handling function };; { .mmi - nop.m 0 - nop.m 0 add GR_Parameter_RESULT = 48,sp + nop.m 0 + nop.i 0 };; { .mmi ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack @@ -1613,8 +1425,7 @@ __libm_error_region: mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs br.ret.sptk b0 // Return };; -.endp __libm_error_region -ASM_SIZE_DIRECTIVE(__libm_error_region) +LOCAL_LIBM_END(__libm_error_region#) .type __libm_error_support#,@function .global __libm_error_support# |