summaryrefslogtreecommitdiff
path: root/sysdeps/ia64/fpu/libm_reduce.S
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ia64/fpu/libm_reduce.S')
-rw-r--r--sysdeps/ia64/fpu/libm_reduce.S1492
1 files changed, 765 insertions, 727 deletions
diff --git a/sysdeps/ia64/fpu/libm_reduce.S b/sysdeps/ia64/fpu/libm_reduce.S
index 1c7f4e1e88..8bdf91d6de 100644
--- a/sysdeps/ia64/fpu/libm_reduce.S
+++ b/sysdeps/ia64/fpu/libm_reduce.S
@@ -1,10 +1,10 @@
.file "libm_reduce.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,304 +20,310 @@
// * 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
+
+// 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
-// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
+// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
-// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-//
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
// 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.
+// problem reports or change requests be submitted to it directly at
+// http://www.intel.com/software/products/opensource/libraries/num.htm.
//
-// History: 02/02/00 Initial Version
+// History:
+// 02/02/00 Initial Version
+// 05/13/02 Rescheduled for speed, changed interface to pass
+// parameters in fp registers
+// 02/10/03 Reordered header: .section, .global, .proc, .align;
+// used data8 for long double data storage
//
-// *********************************************************************
-// *********************************************************************
+//*********************************************************************
+//*********************************************************************
//
// Function: __libm_pi_by_two_reduce(x) return r, c, and N where
// x = N * pi/4 + (r+c) , where |r+c| <= pi/4.
// This function is not designed to be used by the
// general user.
//
-// *********************************************************************
+//*********************************************************************
//
// Accuracy: Returns double-precision values
//
-// *********************************************************************
+//*********************************************************************
//
// Resources Used:
//
-// Floating-Point Registers: f32-f70
+// Floating-Point Registers:
+// f8 = Input x, return value r
+// f9 = return value c
+// f32-f70
//
// General Purpose Registers:
// r8 = return value N
-// r32 = Address of x
-// r33 = Address of where to place r and then c
// r34-r64
//
// Predicate Registers: p6-p14
//
-// *********************************************************************
+//*********************************************************************
//
// IEEE Special Conditions:
//
-// No condions should be raised.
+// No condions should be raised.
//
-// *********************************************************************
+//*********************************************************************
//
// I. Introduction
// ===============
//
// For the forward trigonometric functions sin, cos, sincos, and
-// tan, the original algorithms for IA 64 handle arguments up to
+// tan, the original algorithms for IA 64 handle arguments up to
// 1 ulp less than 2^63 in magnitude. For double-extended arguments x,
-// |x| >= 2^63, this routine returns CASE, N and r_hi, r_lo where
-//
+// |x| >= 2^63, this routine returns N and r_hi, r_lo where
+//
// x is accurately approximated by
// 2*K*pi + N * pi/2 + r_hi + r_lo, |r_hi+r_lo| <= pi/4.
// CASE = 1 or 2.
// CASE is 1 unless |r_hi + r_lo| < 2^(-33).
-//
+//
// The exact value of K is not determined, but that information is
// not required in trigonometric function computations.
-//
-// We first assume the argument x in question satisfies x >= 2^(63).
+//
+// We first assume the argument x in question satisfies x >= 2^(63).
// In particular, it is positive. Negative x can be handled by symmetry:
-//
+//
// -x is accurately approximated by
// -2*K*pi + (-N) * pi/2 - (r_hi + r_lo), |r_hi+r_lo| <= pi/4.
-//
+//
// The idea of the reduction is that
-//
-// x * 2/pi = N_big + N + f, |f| <= 1/2
-//
+//
+// x * 2/pi = N_big + N + f, |f| <= 1/2
+//
// Moreover, for double extended x, |f| >= 2^(-75). (This is an
// non-obvious fact found by enumeration using a special algorithm
-// involving continued fraction.) The algorithm described below
+// involving continued fraction.) The algorithm described below
// calculates N and an accurate approximation of f.
-//
-// Roughly speaking, an appropriate 256-bit (4 X 64) portion of
+//
+// Roughly speaking, an appropriate 256-bit (4 X 64) portion of
// 2/pi is multiplied with x to give the desired information.
-//
+//
// II. Representation of 2/PI
// ==========================
-//
+//
// The value of 2/pi in binary fixed-point is
-//
+//
// .101000101111100110......
-//
+//
// We store 2/pi in a table, starting at the position corresponding
-// to bit position 63
-//
+// to bit position 63
+//
// bit position 63 62 ... 0 -1 -2 -3 -4 -5 -6 -7 .... -16576
-//
-// 0 0 ... 0 . 1 0 1 0 1 0 1 .... X
-//
+//
+// 0 0 ... 0 . 1 0 1 0 1 0 1 .... X
+//
// ^
-// |__ implied binary pt
-//
+// |__ implied binary pt
+//
// III. Algorithm
// ==============
-//
+//
// This describes the algorithm in the most natural way using
-// unsigned interger multiplication. The implementation section
+// unsigned interger multiplication. The implementation section
// describes how the integer arithmetic is simulated.
-//
+//
// STEP 0. Initialization
// ----------------------
-//
-// Let the input argument x be
-//
+//
+// Let the input argument x be
+//
// x = 2^m * ( 1. b_1 b_2 b_3 ... b_63 ), 63 <= m <= 16383.
-//
-// The first crucial step is to fetch four 64-bit portions of 2/pi.
+//
+// The first crucial step is to fetch four 64-bit portions of 2/pi.
// To fulfill this goal, we calculate the bit position L of the
// beginning of these 256-bit quantity by
-//
+//
// L := 62 - m.
-//
-// Note that -16321 <= L <= -1 because 63 <= m <= 16383; and that
+//
+// Note that -16321 <= L <= -1 because 63 <= m <= 16383; and that
// the storage of 2/pi is adequate.
-//
+//
// Fetch P_1, P_2, P_3, P_4 beginning at bit position L thus:
-//
+//
// bit position L L-1 L-2 ... L-63
-//
+//
// P_1 = b b b ... b
-//
+//
// each b can be 0 or 1. Also, let P_0 be the two bits correspoding to
// bit positions L+2 and L+1. So, when each of the P_j is interpreted
// with appropriate scaling, we have
//
// 2/pi = P_big + P_0 + (P_1 + P_2 + P_3 + P_4) + P_small
-//
+//
// Note that P_big and P_small can be ignored. The reasons are as follow.
// First, consider P_big. If P_big = 0, we can certainly ignore it.
-// Otherwise, P_big >= 2^(L+3). Now,
-//
+// Otherwise, P_big >= 2^(L+3). Now,
+//
// P_big * ulp(x) >= 2^(L+3) * 2^(m-63)
-// >= 2^(65-m + m-63 )
-// >= 2^2
-//
+// >= 2^(65-m + m-63 )
+// >= 2^2
+//
// Thus, P_big * x is an integer of the form 4*K. So
-//
-// x = 4*K * (pi/2) + x*(P_0 + P_1 + P_2 + P_3 + P_4)*(pi/2)
+//
+// x = 4*K * (pi/2) + x*(P_0 + P_1 + P_2 + P_3 + P_4)*(pi/2)
// + x*P_small*(pi/2).
-//
+//
// Hence, P_big*x corresponds to information that can be ignored for
// trigonometic function evaluation.
-//
+//
// Next, we must estimate the effect of ignoring P_small. The absolute
// error made by ignoring P_small is bounded by
-//
+//
// |P_small * x| <= ulp(P_4) * x
-// <= 2^(L-255) * 2^(m+1)
-// <= 2^(62-m-255 + m + 1)
-// <= 2^(-192)
-//
-// Since for double-extended precision, x * 2/pi = integer + f,
+// <= 2^(L-255) * 2^(m+1)
+// <= 2^(62-m-255 + m + 1)
+// <= 2^(-192)
+//
+// Since for double-extended precision, x * 2/pi = integer + f,
// 0.5 >= |f| >= 2^(-75), the relative error introduced by ignoring
// P_small is bounded by 2^(-192+75) <= 2^(-117), which is acceptable.
-//
+//
// Further note that if x is split into x_hi + x_lo where x_lo is the
// two bits corresponding to bit positions 2^(m-62) and 2^(m-63); then
-//
-// P_0 * x_hi
-//
+//
+// P_0 * x_hi
+//
// is also an integer of the form 4*K; and thus can also be ignored.
// Let M := P_0 * x_lo which is a small integer. The main part of the
// calculation is really the multiplication of x with the four pieces
// P_1, P_2, P_3, and P_4.
-//
+//
// Unless the reduced argument is extremely small in magnitude, it
// suffices to carry out the multiplication of x with P_1, P_2, and
-// P_3. x*P_4 will be carried out and added on as a correction only
+// P_3. x*P_4 will be carried out and added on as a correction only
// when it is found to be needed. Note also that x*P_4 need not be
// computed exactly. A straightforward multiplication suffices since
// the rounding error thus produced would be bounded by 2^(-3*64),
// that is 2^(-192) which is small enough as the reduced argument
// is bounded from below by 2^(-75).
-//
+//
// Now that we have four 64-bit data representing 2/pi and a
// 64-bit x. We first need to calculate a highly accurate product
// of x and P_1, P_2, P_3. This is best understood as integer
// multiplication.
-//
-//
+//
+//
// STEP 1. Multiplication
// ----------------------
-//
-//
+//
+//
// --------- --------- ---------
-// | P_1 | | P_2 | | P_3 |
-// --------- --------- ---------
-//
+// | P_1 | | P_2 | | P_3 |
+// --------- --------- ---------
+//
+// ---------
+// X | X |
// ---------
-// X | X |
-// ---------
// ----------------------------------------------------
//
// --------- ---------
-// | A_hi | | A_lo |
-// --------- ---------
+// | A_hi | | A_lo |
+// --------- ---------
//
//
// --------- ---------
-// | B_hi | | B_lo |
-// --------- ---------
+// | B_hi | | B_lo |
+// --------- ---------
//
//
-// --------- ---------
-// | C_hi | | C_lo |
-// --------- ---------
+// --------- ---------
+// | C_hi | | C_lo |
+// --------- ---------
//
// ====================================================
// --------- --------- --------- ---------
-// | S_0 | | S_1 | | S_2 | | S_3 |
-// --------- --------- --------- ---------
+// | S_0 | | S_1 | | S_2 | | S_3 |
+// --------- --------- --------- ---------
//
//
//
// STEP 2. Get N and f
// -------------------
-//
+//
// Conceptually, after the individual pieces S_0, S_1, ..., are obtained,
// we have to sum them and obtain an integer part, N, and a fraction, f.
// Here, |f| <= 1/2, and N is an integer. Note also that N need only to
// be known to module 2^k, k >= 2. In the case when |f| is small enough,
// we would need to add in the value x*P_4.
-//
-//
+//
+//
// STEP 3. Get reduced argument
// ----------------------------
-//
+//
// The value f is not yet the reduced argument that we seek. The
// equation
-//
-// x * 2/pi = 4K + N + f
-//
+//
+// x * 2/pi = 4K + N + f
+//
// says that
-//
+//
// x = 2*K*pi + N * pi/2 + f * (pi/2).
-//
+//
// Thus, the reduced argument is given by
-//
-// reduced argument = f * pi/2.
-//
+//
+// reduced argument = f * pi/2.
+//
// This multiplication must be performed to extra precision.
-//
+//
// IV. Implementation
// ==================
-//
+//
// Step 0. Initialization
// ----------------------
-//
+//
// Set sgn_x := sign(x); x := |x|; x_lo := 2 lsb of x.
-//
+//
// In memory, 2/pi is stored contigously as
-//
+//
// 0x00000000 0x00000000 0xA2F....
// ^
// |__ implied binary bit
-//
+//
// Given x = 2^m * 1.xxxx...xxx; we calculate L := 62 - m. Thus
// -1 <= L <= -16321. We fetch from memory 5 integer pieces of data.
-//
+//
// P_0 is the two bits corresponding to bit positions L+2 and L+1
// P_1 is the 64-bit starting at bit position L
// P_2 is the 64-bit starting at bit position L-64
// P_3 is the 64-bit starting at bit position L-128
// P_4 is the 64-bit starting at bit position L-192
-//
+//
// For example, if m = 63, P_0 would be 0 and P_1 would look like
// 0xA2F...
-//
+//
// If m = 65, P_0 would be the two msb of 0xA, thus, P_0 is 10 in binary.
-// P_1 in binary would be 1 0 0 0 1 0 1 1 1 1 ....
-//
+// P_1 in binary would be 1 0 0 0 1 0 1 1 1 1 ....
+//
// Step 1. Multiplication
// ----------------------
-//
+//
// At this point, P_1, P_2, P_3, P_4 are integers. They are
// supposed to be interpreted as
-//
+//
// 2^(L-63) * P_1;
// 2^(L-63-64) * P_2;
// 2^(L-63-128) * P_3;
// 2^(L-63-192) * P_4;
-//
+//
// Since each of them need to be multiplied to x, we would scale
// both x and the P_j's by some convenient factors: scale each
// of P_j's up by 2^(63-L), and scale x down by 2^(L-63).
-//
+//
// p_1 := fcvt.xf ( P_1 )
// p_2 := fcvt.xf ( P_2 ) * 2^(-64)
// p_3 := fcvt.xf ( P_3 ) * 2^(-128)
@@ -325,30 +331,30 @@
// x := replace exponent of x by -1
// because 2^m * 1.xxxx...xxx * 2^(L-63)
// is 2^(-1) * 1.xxxx...xxx
-//
+//
// We are now faced with the task of computing the following
-//
+//
// --------- --------- ---------
-// | P_1 | | P_2 | | P_3 |
-// --------- --------- ---------
-//
+// | P_1 | | P_2 | | P_3 |
+// --------- --------- ---------
+//
// ---------
-// X | X |
-// ---------
+// X | X |
+// ---------
// ----------------------------------------------------
-//
+//
// --------- ---------
-// | A_hi | | A_lo |
-// --------- ---------
-//
+// | A_hi | | A_lo |
+// --------- ---------
+//
// --------- ---------
-// | B_hi | | B_lo |
-// --------- ---------
-//
-// --------- ---------
-// | C_hi | | C_lo |
-// --------- ---------
-//
+// | B_hi | | B_lo |
+// --------- ---------
+//
+// --------- ---------
+// | C_hi | | C_lo |
+// --------- ---------
+//
// ====================================================
// ----------- --------- --------- ---------
// | S_0 | | S_1 | | S_2 | | S_3 |
@@ -357,108 +363,108 @@
// | |___ binary point
// |
// |___ possibly one more bit
-//
+//
// Let FPSR3 be set to round towards zero with widest precision
-// and exponent range. Unless an explicit FPSR is given,
+// and exponent range. Unless an explicit FPSR is given,
// round-to-nearest with widest precision and exponent range is
// used.
-//
+//
// Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_C := 2^(-65).
-//
+//
// Tmp_C := fmpy.fpsr3( x, p_1 );
// If Tmp_C >= sigma_C then
// C_hi := Tmp_C;
// C_lo := x*p_1 - C_hi ...fma, exact
// Else
// C_hi := fadd.fpsr3(sigma_C, Tmp_C) - sigma_C
-// ...subtraction is exact, regardless
-// ...of rounding direction
+// ...subtraction is exact, regardless
+// ...of rounding direction
// C_lo := x*p_1 - C_hi ...fma, exact
// End If
-//
+//
// Tmp_B := fmpy.fpsr3( x, p_2 );
// If Tmp_B >= sigma_B then
// B_hi := Tmp_B;
// B_lo := x*p_2 - B_hi ...fma, exact
// Else
// B_hi := fadd.fpsr3(sigma_B, Tmp_B) - sigma_B
-// ...subtraction is exact, regardless
-// ...of rounding direction
+// ...subtraction is exact, regardless
+// ...of rounding direction
// B_lo := x*p_2 - B_hi ...fma, exact
// End If
-//
+//
// Tmp_A := fmpy.fpsr3( x, p_3 );
// If Tmp_A >= sigma_A then
// A_hi := Tmp_A;
// A_lo := x*p_3 - A_hi ...fma, exact
// Else
// A_hi := fadd.fpsr3(sigma_A, Tmp_A) - sigma_A
-// ...subtraction is exact, regardless
-// ...of rounding direction
+// ...subtraction is exact, regardless
+// ...of rounding direction
// A_lo := x*p_3 - A_hi ...fma, exact
// End If
-//
+//
// ...Note that C_hi is of integer value. We need only the
-// ...last few bits. Thus we can ensure C_hi is never a big
+// ...last few bits. Thus we can ensure C_hi is never a big
// ...integer, freeing us from overflow worry.
-//
+//
// Tmp_C := fadd.fpsr3( C_hi, 2^(70) ) - 2^(70);
// ...Tmp_C is the upper portion of C_hi
// C_hi := C_hi - Tmp_C
// ...0 <= C_hi < 2^7
-//
+//
// Step 2. Get N and f
// -------------------
-//
-// At this point, we have all the components to obtain
+//
+// At this point, we have all the components to obtain
// S_0, S_1, S_2, S_3 and thus N and f. We start by adding
// C_lo and B_hi. This sum together with C_hi gives a good
-// estimation of N and f.
-//
+// estimation of N and f.
+//
// A := fadd.fpsr3( B_hi, C_lo )
// B := max( B_hi, C_lo )
// b := min( B_hi, C_lo )
-//
-// a := (B - A) + b ...exact. Note that a is either 0
-// ...or 2^(-64).
-//
+//
+// a := (B - A) + b ...exact. Note that a is either 0
+// ...or 2^(-64).
+//
// N := round_to_nearest_integer_value( A );
-// f := A - N; ...exact because lsb(A) >= 2^(-64)
-// ...and |f| <= 1/2.
-//
-// f := f + a ...exact because a is 0 or 2^(-64);
-// ...the msb of the sum is <= 1/2
-// ...lsb >= 2^(-64).
-//
+// f := A - N; ...exact because lsb(A) >= 2^(-64)
+// ...and |f| <= 1/2.
+//
+// f := f + a ...exact because a is 0 or 2^(-64);
+// ...the msb of the sum is <= 1/2
+// ...lsb >= 2^(-64).
+//
// N := convert to integer format( C_hi + N );
// M := P_0 * x_lo;
// N := N + M;
-//
+//
// If sgn_x == 1 (that is original x was negative)
// N := 2^10 - N
// ...this maintains N to be non-negative, but still
// ...equivalent to the (negated N) mod 4.
// End If
-//
+//
// If |f| >= 2^(-33)
-//
+//
// ...Case 1
// CASE := 1
// g := A_hi + B_lo;
// s_hi := f + g;
// s_lo := (f - s_hi) + g;
-//
+//
// Else
-//
+//
// ...Case 2
// CASE := 2
// A := fadd.fpsr3( A_hi, B_lo )
// B := max( A_hi, B_lo )
// b := min( A_hi, B_lo )
-//
-// a := (B - A) + b ...exact. Note that a is either 0
-// ...or 2^(-128).
-//
+//
+// a := (B - A) + b ...exact. Note that a is either 0
+// ...or 2^(-128).
+//
// f_hi := A + f;
// f_lo := (f - f_hi) + A;
// ...this is exact.
@@ -468,9 +474,9 @@
// ...If f = 2^(-64), f-f_hi involves cancellation and is
// ...exact. If f = -2^(-64), then A + f is exact. Hence
// ...f-f_hi is -A exactly, giving f_lo = 0.
-//
+//
// f_lo := f_lo + a;
-//
+//
// If |f| >= 2^(-50) then
// s_hi := f_hi;
// s_lo := f_lo;
@@ -479,117 +485,111 @@
// s_hi := f_hi + f_lo
// s_lo := (f_hi - s_hi) + f_lo
// End If
-//
+//
// End If
-//
+//
// Step 3. Get reduced argument
// ----------------------------
-//
+//
// If sgn_x == 0 (that is original x is positive)
-//
+//
// D_hi := Pi_by_2_hi
// D_lo := Pi_by_2_lo
// ...load from table
-//
+//
// Else
-//
+//
// D_hi := neg_Pi_by_2_hi
// D_lo := neg_Pi_by_2_lo
// ...load from table
// End If
-//
+//
// r_hi := s_hi*D_hi
-// r_lo := s_hi*D_hi - r_hi ...fma
+// r_lo := s_hi*D_hi - r_hi ...fma
// r_lo := (s_hi*D_lo + r_lo) + s_lo*D_hi
-//
-// Return CASE, N, r_hi, r_lo
-//
-
-#include "libm_support.h"
-
-FR_X = f32
-FR_N = f33
-FR_p_1 = f34
-FR_TWOM33 = f35
-FR_TWOM50 = f36
-FR_g = f37
-FR_p_2 = f38
-FR_f = f39
-FR_s_lo = f40
-FR_p_3 = f41
-FR_f_abs = f42
-FR_D_lo = f43
-FR_p_4 = f44
-FR_D_hi = f45
-FR_Tmp2_C = f46
-FR_s_hi = f47
-FR_sigma_A = f48
-FR_A = f49
-FR_sigma_B = f50
-FR_B = f51
-FR_sigma_C = f52
-FR_b = f53
-FR_ScaleP2 = f54
-FR_ScaleP3 = f55
-FR_ScaleP4 = f56
-FR_Tmp_A = f57
-FR_Tmp_B = f58
-FR_Tmp_C = f59
-FR_A_hi = f60
-FR_f_hi = f61
-FR_r_hi = f62
-FR_A_lo = f63
-FR_B_hi = f64
-FR_a = f65
-FR_B_lo = f66
+//
+// Return N, r_hi, r_lo
+//
+FR_input_X = f8
+FR_r_hi = f8
+FR_r_lo = f9
+
+FR_X = f32
+FR_N = f33
+FR_p_1 = f34
+FR_TWOM33 = f35
+FR_TWOM50 = f36
+FR_g = f37
+FR_p_2 = f38
+FR_f = f39
+FR_s_lo = f40
+FR_p_3 = f41
+FR_f_abs = f42
+FR_D_lo = f43
+FR_p_4 = f44
+FR_D_hi = f45
+FR_Tmp2_C = f46
+FR_s_hi = f47
+FR_sigma_A = f48
+FR_A = f49
+FR_sigma_B = f50
+FR_B = f51
+FR_sigma_C = f52
+FR_b = f53
+FR_ScaleP2 = f54
+FR_ScaleP3 = f55
+FR_ScaleP4 = f56
+FR_Tmp_A = f57
+FR_Tmp_B = f58
+FR_Tmp_C = f59
+FR_A_hi = f60
+FR_f_hi = f61
+FR_RSHF = f62
+FR_A_lo = f63
+FR_B_hi = f64
+FR_a = f65
+FR_B_lo = f66
FR_f_lo = f67
-FR_r_lo = f68
-FR_C_hi = f69
-FR_C_lo = f70
+FR_N_fix = f68
+FR_C_hi = f69
+FR_C_lo = f70
GR_N = r8
-GR_Address_of_Input = r32
-GR_Address_of_Outputs = r33
-GR_Exp_x = r36
-GR_Temp = r37
-GR_BIASL63 = r38
+GR_Exp_x = r36
+GR_Temp = r37
+GR_BIASL63 = r38
GR_CASE = r39
-GR_x_lo = r40
-GR_sgn_x = r41
+GR_x_lo = r40
+GR_sgn_x = r41
GR_M = r42
GR_BASE = r43
GR_LENGTH1 = r44
GR_LENGTH2 = r45
GR_ASUB = r46
GR_P_0 = r47
-GR_P_1 = r48
-GR_P_2 = r49
-GR_P_3 = r50
-GR_P_4 = r51
+GR_P_1 = r48
+GR_P_2 = r49
+GR_P_3 = r50
+GR_P_4 = r51
GR_START = r52
GR_SEGMENT = r53
GR_A = r54
-GR_B = r55
+GR_B = r55
GR_C = r56
GR_D = r57
GR_E = r58
-GR_TEMP1 = r59
-GR_TEMP2 = r60
-GR_TEMP3 = r61
-GR_TEMP4 = r62
+GR_TEMP1 = r59
+GR_TEMP2 = r60
+GR_TEMP3 = r61
+GR_TEMP4 = r62
GR_TEMP5 = r63
GR_TEMP6 = r64
+GR_rshf = r64
+RODATA
.align 64
-#ifdef _LIBC
-.rodata
-#else
-.data
-#endif
-
-Constants_Bits_of_2_by_pi:
-ASM_TYPE_DIRECTIVE(Constants_Bits_of_2_by_pi,@object)
+LOCAL_OBJECT_START(Constants_Bits_of_2_by_pi)
data8 0x0000000000000000,0xA2F9836E4E441529
data8 0xFC2757D1F534DDC0,0xDB6295993C439041
data8 0xFE5163ABDEBBC561,0xB7246E3A424DD2E0
@@ -721,34 +721,33 @@ data8 0xB5D6DF8261DD9602,0x36169F3AC4A1A283
data8 0x6DED727A8D39A9B8,0x825C326B5B2746ED
data8 0x34007700D255F4FC,0x4D59018071E0E13F
data8 0x89B295F364A8F1AE,0xA74B38FC4CEAB2BB
-ASM_SIZE_DIRECTIVE(Constants_Bits_of_2_by_pi)
+LOCAL_OBJECT_END(Constants_Bits_of_2_by_pi)
-Constants_Bits_of_pi_by_2:
-ASM_TYPE_DIRECTIVE(Constants_Bits_of_pi_by_2,@object)
-data4 0x2168C234,0xC90FDAA2,0x00003FFF,0x00000000
-data4 0x80DC1CD1,0xC4C6628B,0x00003FBF,0x00000000
-ASM_SIZE_DIRECTIVE(Constants_Bits_of_pi_by_2)
+LOCAL_OBJECT_START(Constants_Bits_of_pi_by_2)
+data8 0xC90FDAA22168C234,0x00003FFF
+data8 0xC4C6628B80DC1CD1,0x00003FBF
+LOCAL_OBJECT_END(Constants_Bits_of_pi_by_2)
.section .text
-.proc __libm_pi_by_2_reduce#
.global __libm_pi_by_2_reduce#
-.align 64
+.proc __libm_pi_by_2_reduce#
+.align 32
-__libm_pi_by_2_reduce:
+__libm_pi_by_2_reduce:
-// X is at the address in Address_of_Input
-// Place the two-piece result at the address in Address_of_Outputs
-// r followed by c
-// N is returned
+// X is in f8
+// Place the two-piece result r (r_hi) in f8 and c (r_lo) in f9
+// N is returned in r8
-{ .mmf
-alloc r34 = ar.pfs,2,34,0,0
-(p0) ldfe FR_X = [GR_Address_of_Input]
-(p0) fsetc.s3 0x00,0x7F ;;
+{ .mfi
+ alloc r34 = ar.pfs,2,34,0,0
+ fsetc.s3 0x00,0x7F // Set sf3 to round to zero, 82-bit prec, td, ftz
+ nop.i 999
}
-{ .mlx
- nop.m 999
-(p0) movl GR_BIASL63 = 0x1003E
+{ .mfi
+ addl GR_BASE = @ltoff(Constants_Bits_of_2_by_pi#), gp
+ nop.f 999
+ mov GR_BIASL63 = 0x1003E
}
;;
@@ -765,73 +764,61 @@ alloc r34 = ar.pfs,2,34,0,0
// Address_BASE = shladd(SEGMENT,3) + BASE
-
{ .mmi
- nop.m 999
-(p0) addl GR_BASE = @ltoff(Constants_Bits_of_2_by_pi#), gp
- nop.i 999
+ getf.exp GR_Exp_x = FR_input_X
+ ld8 GR_BASE = [GR_BASE]
+ mov GR_TEMP5 = 0x0FFFE
}
;;
+// Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_A := 2^(-65).
{ .mmi
- ld8 GR_BASE = [GR_BASE]
- nop.m 999
+ getf.sig GR_x_lo = FR_input_X
+ mov GR_TEMP6 = 0x0FFBE
nop.i 999
}
;;
-
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP5 = 0x000000000000FFFE
-}
-{ .mmi
- nop.m 999 ;;
-(p0) setf.exp FR_sigma_B = GR_TEMP5
- nop.i 999
-}
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP6 = 0x000000000000FFBE ;;
-}
-// Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_A := 2^(-65).
-{ .mfi
-(p0) setf.exp FR_sigma_A = GR_TEMP6
- nop.f 999
- nop.i 999 ;;
-}
-// Special Code for testing DE arguments
-// (p0) movl GR_BIASL63 = 0x0000000000013FFE
-// (p0) movl GR_x_lo = 0xFFFFFFFFFFFFFFFF
-// (p0) setf.exp FR_X = GR_BIASL63
-// (p0) setf.sig FR_ScaleP3 = GR_x_lo
-// (p0) fmerge.se FR_X = FR_X,FR_ScaleP3
+// Special Code for testing DE arguments
+// movl GR_BIASL63 = 0x0000000000013FFE
+// movl GR_x_lo = 0xFFFFFFFFFFFFFFFF
+// setf.exp FR_X = GR_BIASL63
+// setf.sig FR_ScaleP3 = GR_x_lo
+// fmerge.se FR_X = FR_X,FR_ScaleP3
// Set sgn_x := sign(x); x := |x|; x_lo := 2 lsb of x.
// 2/pi is stored contigously as
// 0x00000000 0x00000000.0xA2F....
// M = EXP - BIAS ( M >= 63)
// Given x = 2^m * 1.xxxx...xxx; we calculate L := 62 - m.
// Thus -1 <= L <= -16321.
-{ .mmf
-(p0) getf.exp GR_Exp_x = FR_X
-(p0) getf.sig GR_x_lo = FR_X
-(p0) fabs FR_X = FR_X ;;
+{ .mmi
+ setf.exp FR_sigma_B = GR_TEMP5
+ setf.exp FR_sigma_A = GR_TEMP6
+ extr.u GR_M = GR_Exp_x,0,17
}
+;;
+
{ .mii
-(p0) and GR_x_lo = 0x03,GR_x_lo
-(p0) extr.u GR_M = GR_Exp_x,0,17 ;;
-(p0) sub GR_START = GR_M,GR_BIASL63
+ and GR_x_lo = 0x03,GR_x_lo
+ sub GR_START = GR_M,GR_BIASL63
+ add GR_BASE = 8,GR_BASE // To effectively add 1 to SEGMENT
}
-{ .mmi
- nop.m 999 ;;
-(p0) and GR_LENGTH1 = 0x3F,GR_START
-(p0) shr.u GR_SEGMENT = GR_START,6
+;;
+
+{ .mii
+ and GR_LENGTH1 = 0x3F,GR_START
+ shr.u GR_SEGMENT = GR_START,6
+ nop.i 999
}
+;;
+
{ .mmi
- nop.m 999 ;;
-(p0) add GR_SEGMENT = 0x1,GR_SEGMENT
-(p0) sub GR_LENGTH2 = 0x40,GR_LENGTH1
+ shladd GR_BASE = GR_SEGMENT,3,GR_BASE
+ sub GR_LENGTH2 = 0x40,GR_LENGTH1
+ cmp.le p6,p7 = 0x2,GR_LENGTH1
}
+;;
+
// P_0 is the two bits corresponding to bit positions L+2 and L+1
// P_1 is the 64-bit starting at bit position L
// P_2 is the 64-bit starting at bit position L-64
@@ -849,13 +836,13 @@ alloc r34 = ar.pfs,2,34,0,0
// P_4 is made up of Clo and Dhi
// P_4 = deposit Dlo, position 0, length2 into P_4, position length1
// deposit Ehi, position length2, length1 into P_4, position 0
-{ .mmi
-(p0) cmp.le.unc p6,p7 = 0x2,GR_LENGTH1 ;;
-(p0) shladd GR_BASE = GR_SEGMENT,3,GR_BASE
-(p7) cmp.eq.unc p8,p9 = 0x1,GR_LENGTH1 ;;
+{ .mfi
+ ld8 GR_A = [GR_BASE],8
+ fabs FR_X = FR_input_X
+(p7) cmp.eq.unc p8,p9 = 0x1,GR_LENGTH1
}
-{ .mmi
- nop.m 999
+;;
+
// ld_64 A at Base and increment Base by 8
// ld_64 B at Base and increment Base by 8
// ld_64 C at Base and increment Base by 8
@@ -866,31 +853,35 @@ alloc r34 = ar.pfs,2,34,0,0
// A, B, C, D, and E look like | length1 | length2 |
// ---------------------
// hi lo
-(p0) ld8 GR_A = [GR_BASE],8
-(p0) extr.u GR_sgn_x = GR_Exp_x,17,1 ;;
-}
-{ .mmf
- nop.m 999
-(p0) ld8 GR_B = [GR_BASE],8
-(p0) fmerge.se FR_X = FR_sigma_B,FR_X ;;
+{ .mlx
+ ld8 GR_B = [GR_BASE],8
+ movl GR_rshf = 0x43e8000000000000 // 1.10000 2^63 for right shift N_fix
}
-{ .mii
-(p0) ld8 GR_C = [GR_BASE],8
-(p8) extr.u GR_Temp = GR_A,63,1 ;;
-(p0) shl GR_TEMP1 = GR_A,GR_LENGTH1
+;;
+
+{ .mmi
+ ld8 GR_C = [GR_BASE],8
+ nop.m 999
+(p8) extr.u GR_Temp = GR_A,63,1
}
-{ .mii
-(p0) ld8 GR_D = [GR_BASE],8
+;;
+
// If length1 >= 2,
// P_0 = deposit Ahi, position length2, 2 bit into P_0 at position 0.
-(p6) shr.u GR_P_0 = GR_A,GR_LENGTH2 ;;
-(p0) shl GR_TEMP2 = GR_B,GR_LENGTH1
+{ .mii
+ ld8 GR_D = [GR_BASE],8
+ shl GR_TEMP1 = GR_A,GR_LENGTH1 // MM instruction
+(p6) shr.u GR_P_0 = GR_A,GR_LENGTH2 // MM instruction
}
+;;
+
{ .mii
-(p0) ld8 GR_E = [GR_BASE],-40
-(p0) shr.u GR_P_1 = GR_B,GR_LENGTH2 ;;
-(p0) shr.u GR_P_2 = GR_C,GR_LENGTH2
+ ld8 GR_E = [GR_BASE],-40
+ shl GR_TEMP2 = GR_B,GR_LENGTH1 // MM instruction
+ shr.u GR_P_1 = GR_B,GR_LENGTH2 // MM instruction
}
+;;
+
// Else
// Load 16 bit of ASUB from (Base_Address_of_A - 2)
// P_0 = ASUB & 0x3
@@ -900,43 +891,56 @@ alloc r34 = ar.pfs,2,34,0,0
// Deposit element 63 from Ahi and place in element 0 of P_0.
// Endif
// Endif
+
{ .mii
(p7) ld2 GR_ASUB = [GR_BASE],8
-(p0) shl GR_TEMP3 = GR_C,GR_LENGTH1 ;;
-(p0) shl GR_TEMP4 = GR_D,GR_LENGTH1
+ shl GR_TEMP3 = GR_C,GR_LENGTH1 // MM instruction
+ shr.u GR_P_2 = GR_C,GR_LENGTH2 // MM instruction
}
+;;
+
{ .mii
- nop.m 999
-(p0) shr.u GR_P_3 = GR_D,GR_LENGTH2 ;;
-(p0) shr.u GR_P_4 = GR_E,GR_LENGTH2
+ setf.d FR_RSHF = GR_rshf // Form right shift const 1.100 * 2^63
+ shl GR_TEMP4 = GR_D,GR_LENGTH1 // MM instruction
+ shr.u GR_P_3 = GR_D,GR_LENGTH2 // MM instruction
}
-{ .mii
+;;
+
+{ .mmi
(p7) and GR_P_0 = 0x03,GR_ASUB
-(p6) and GR_P_0 = 0x03,GR_P_0 ;;
-(p0) or GR_P_1 = GR_P_1,GR_TEMP1
+(p6) and GR_P_0 = 0x03,GR_P_0
+ shr.u GR_P_4 = GR_E,GR_LENGTH2 // MM instruction
}
+;;
+
{ .mmi
-(p8) and GR_P_0 = 0x1,GR_P_0 ;;
-(p0) or GR_P_2 = GR_P_2,GR_TEMP2
-(p8) shl GR_P_0 = GR_P_0,0x1 ;;
-}
-{ .mii
- nop.m 999
-(p0) or GR_P_3 = GR_P_3,GR_TEMP3
-(p8) or GR_P_0 = GR_P_0,GR_Temp
+ nop.m 999
+ or GR_P_1 = GR_P_1,GR_TEMP1
+(p8) and GR_P_0 = 0x1,GR_P_0
}
+;;
+
{ .mmi
-(p0) setf.sig FR_p_1 = GR_P_1 ;;
-(p0) setf.sig FR_p_2 = GR_P_2
-(p0) or GR_P_4 = GR_P_4,GR_TEMP4 ;;
+ setf.sig FR_p_1 = GR_P_1
+ or GR_P_2 = GR_P_2,GR_TEMP2
+(p8) shladd GR_P_0 = GR_P_0,1,GR_Temp
}
+;;
+
+{ .mmf
+ setf.sig FR_p_2 = GR_P_2
+ or GR_P_3 = GR_P_3,GR_TEMP3
+ fmerge.se FR_X = FR_sigma_B,FR_X
+}
+;;
+
{ .mmi
- nop.m 999 ;;
-(p0) setf.sig FR_p_3 = GR_P_3
-(p0) pmpy2.r GR_M = GR_P_0,GR_x_lo
+ setf.sig FR_p_3 = GR_P_3
+ or GR_P_4 = GR_P_4,GR_TEMP4
+ pmpy2.r GR_M = GR_P_0,GR_x_lo
}
-{ .mlx
-(p0) setf.sig FR_p_4 = GR_P_4
+;;
+
// P_1, P_2, P_3, P_4 are integers. They should be
// 2^(L-63) * P_1;
// 2^(L-63-64) * P_2;
@@ -954,18 +958,18 @@ alloc r34 = ar.pfs,2,34,0,0
// | P_1 | | P_2 | | P_3 |
// --------- --------- ---------
// ---------
-// X | X |
-// ---------
+// X | X |
+// ---------
// ----------------------------------------------------
// --------- ---------
-// | A_hi | | A_lo |
-// --------- ---------
+// | A_hi | | A_lo |
+// --------- ---------
// --------- ---------
-// | B_hi | | B_lo |
-// --------- ---------
+// | B_hi | | B_lo |
+// --------- ---------
+// --------- ---------
+// | C_hi | | C_lo |
// --------- ---------
-// | C_hi | | C_lo |
-// --------- ---------
// ====================================================
// ----------- --------- --------- ---------
// | S_0 | | S_1 | | S_2 | | S_3 |
@@ -977,52 +981,55 @@ alloc r34 = ar.pfs,2,34,0,0
// and exponent range. Unless an explicit FPSR is given,
// round-to-nearest with widest precision and exponent range is
// used.
-(p0) movl GR_TEMP1 = 0x000000000000FFBF
-}
{ .mmi
- nop.m 999 ;;
-(p0) setf.exp FR_ScaleP2 = GR_TEMP1
- nop.i 999
-}
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP4 = 0x000000000001003E
+ setf.sig FR_p_4 = GR_P_4
+ mov GR_TEMP1 = 0x0FFBF
+ nop.i 999
}
+;;
+
{ .mmi
- nop.m 999 ;;
-(p0) setf.exp FR_sigma_C = GR_TEMP4
- nop.i 999
+ setf.exp FR_ScaleP2 = GR_TEMP1
+ mov GR_TEMP2 = 0x0FF7F
+ nop.i 999
}
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP2 = 0x000000000000FF7F ;;
+;;
+
+{ .mmi
+ setf.exp FR_ScaleP3 = GR_TEMP2
+ mov GR_TEMP4 = 0x1003E
+ nop.i 999
}
+;;
+
{ .mmf
- nop.m 999
-(p0) setf.exp FR_ScaleP3 = GR_TEMP2
-(p0) fcvt.xuf.s1 FR_p_1 = FR_p_1 ;;
+ setf.exp FR_sigma_C = GR_TEMP4
+ mov GR_Temp = 0x0FFDE
+ fcvt.xuf.s1 FR_p_1 = FR_p_1
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fcvt.xuf.s1 FR_p_2 = FR_p_2
- nop.i 999
-}
-{ .mlx
- nop.m 999
-(p0) movl GR_Temp = 0x000000000000FFDE ;;
-}
-{ .mmf
- nop.m 999
-(p0) setf.exp FR_TWOM33 = GR_Temp
-(p0) fcvt.xuf.s1 FR_p_3 = FR_p_3 ;;
+ setf.exp FR_TWOM33 = GR_Temp
+ fcvt.xuf.s1 FR_p_2 = FR_p_2
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fcvt.xuf.s1 FR_p_4 = FR_p_4
- nop.i 999 ;;
+ nop.m 999
+ fcvt.xuf.s1 FR_p_3 = FR_p_3
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
+ nop.m 999
+ fcvt.xuf.s1 FR_p_4 = FR_p_4
+ nop.i 999
+}
+;;
+
// Tmp_C := fmpy.fpsr3( x, p_1 );
// Tmp_B := fmpy.fpsr3( x, p_2 );
// Tmp_A := fmpy.fpsr3( x, p_3 );
@@ -1048,55 +1055,62 @@ alloc r34 = ar.pfs,2,34,0,0
// Exact, regardless ...of rounding direction
// A_lo := x*p_3 - A_hi ...fma, exact
// Endif
-(p0) fmpy.s3 FR_Tmp_C = FR_X,FR_p_1
- nop.i 999 ;;
-}
{ .mfi
- nop.m 999
-(p0) fmpy.s1 FR_p_2 = FR_p_2,FR_ScaleP2
- nop.i 999
-}
-{ .mlx
- nop.m 999
-(p0) movl GR_Temp = 0x0000000000000400
+ nop.m 999
+ fmpy.s3 FR_Tmp_C = FR_X,FR_p_1
+ nop.i 999
}
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP3 = 0x000000000000FF3F ;;
+;;
+
+{ .mfi
+ mov GR_TEMP3 = 0x0FF3F
+ fmpy.s1 FR_p_2 = FR_p_2,FR_ScaleP2
+ nop.i 999
}
+;;
+
{ .mmf
- nop.m 999
-(p0) setf.exp FR_ScaleP4 = GR_TEMP3
-(p0) fmpy.s1 FR_p_3 = FR_p_3,FR_ScaleP3 ;;
+ setf.exp FR_ScaleP4 = GR_TEMP3
+ mov GR_TEMP4 = 0x10045
+ fmpy.s1 FR_p_3 = FR_p_3,FR_ScaleP3
}
-{ .mlx
- nop.m 999
-(p0) movl GR_TEMP4 = 0x0000000000010045 ;;
+;;
+
+{ .mfi
+ nop.m 999
+ fadd.s3 FR_C_hi = FR_sigma_C,FR_Tmp_C // For Tmp_C < sigma_C case
+ nop.i 999
}
+;;
+
{ .mmf
- nop.m 999
-(p0) setf.exp FR_Tmp2_C = GR_TEMP4
-(p0) fmpy.s3 FR_Tmp_B = FR_X,FR_p_2 ;;
+ setf.exp FR_Tmp2_C = GR_TEMP4
+ nop.m 999
+ fmpy.s3 FR_Tmp_B = FR_X,FR_p_2
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fcmp.ge.unc.s1 p12, p9 = FR_Tmp_C,FR_sigma_C
- nop.i 999 ;;
+ addl GR_BASE = @ltoff(Constants_Bits_of_pi_by_2#), gp
+ fcmp.ge.s1 p12, p9 = FR_Tmp_C,FR_sigma_C
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fmpy.s3 FR_Tmp_A = FR_X,FR_p_3
- nop.i 999 ;;
+ nop.m 999
+ fmpy.s3 FR_Tmp_A = FR_X,FR_p_3
+ nop.i 99
}
+;;
+
{ .mfi
- nop.m 999
+ ld8 GR_BASE = [GR_BASE]
(p12) mov FR_C_hi = FR_Tmp_C
- nop.i 999 ;;
+ nop.i 999
}
{ .mfi
-(p0) addl GR_BASE = @ltoff(Constants_Bits_of_pi_by_2#), gp
-(p9) fadd.s3 FR_C_hi = FR_sigma_C,FR_Tmp_C
- nop.i 999
+ nop.m 999
+(p9) fsub.s1 FR_C_hi = FR_C_hi,FR_sigma_C
+ nop.i 999
}
;;
@@ -1114,97 +1128,106 @@ alloc r34 = ar.pfs,2,34,0,0
// Load from table
// End If
-
-{ .mmi
- ld8 GR_BASE = [GR_BASE]
+{ .mfi
nop.m 999
+ fmpy.s1 FR_p_4 = FR_p_4,FR_ScaleP4
nop.i 999
}
-;;
-
-
{ .mfi
-(p0) ldfe FR_D_hi = [GR_BASE],16
-(p0) fmpy.s1 FR_p_4 = FR_p_4,FR_ScaleP4
- nop.i 999 ;;
+ nop.m 999
+ fadd.s3 FR_B_hi = FR_sigma_B,FR_Tmp_B // For Tmp_B < sigma_B case
+ nop.i 999
}
+;;
+
{ .mfi
-(p0) ldfe FR_D_lo = [GR_BASE],0
-(p0) fcmp.ge.unc.s1 p13, p10 = FR_Tmp_B,FR_sigma_B
- nop.i 999 ;;
+ nop.m 999
+ fadd.s3 FR_A_hi = FR_sigma_A,FR_Tmp_A // For Tmp_A < sigma_A case
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p13) mov FR_B_hi = FR_Tmp_B
- nop.i 999
+ nop.m 999
+ fcmp.ge.s1 p13, p10 = FR_Tmp_B,FR_sigma_B
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p12) fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
- nop.i 999 ;;
+ nop.m 999
+ fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p10) fadd.s3 FR_B_hi = FR_sigma_B,FR_Tmp_B
- nop.i 999
+ ldfe FR_D_hi = [GR_BASE],16
+ fcmp.ge.s1 p14, p11 = FR_Tmp_A,FR_sigma_A
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p9) fsub.s1 FR_C_hi = FR_C_hi,FR_sigma_C
- nop.i 999 ;;
+ ldfe FR_D_lo = [GR_BASE]
+(p13) mov FR_B_hi = FR_Tmp_B
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fcmp.ge.unc.s1 p14, p11 = FR_Tmp_A,FR_sigma_A
- nop.i 999 ;;
+ nop.m 999
+(p10) fsub.s1 FR_B_hi = FR_B_hi,FR_sigma_B
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
+ nop.m 999
(p14) mov FR_A_hi = FR_Tmp_A
- nop.i 999 ;;
-}
-{ .mfi
- nop.m 999
-(p11) fadd.s3 FR_A_hi = FR_sigma_A,FR_Tmp_A
- nop.i 999 ;;
-}
-{ .mfi
- nop.m 999
-(p9) fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
-(p0) cmp.eq.unc p12,p9 = 0x1,GR_sgn_x
-}
-{ .mfi
- nop.m 999
-(p13) fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
- nop.i 999 ;;
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p10) fsub.s1 FR_B_hi = FR_B_hi,FR_sigma_B
- nop.i 999
+ nop.m 999
+(p11) fsub.s1 FR_A_hi = FR_A_hi,FR_sigma_A
+ nop.i 999
}
-{ .mfi
- nop.m 999
+;;
+
// Note that C_hi is of integer value. We need only the
// last few bits. Thus we can ensure C_hi is never a big
// integer, freeing us from overflow worry.
// Tmp_C := fadd.fpsr3( C_hi, 2^(70) ) - 2^(70);
// Tmp_C is the upper portion of C_hi
-(p0) fadd.s3 FR_Tmp_C = FR_C_hi,FR_Tmp2_C
- nop.i 999 ;;
+{ .mfi
+ nop.m 999
+ fadd.s3 FR_Tmp_C = FR_C_hi,FR_Tmp2_C
+ tbit.z p12,p9 = GR_Exp_x, 17
}
+;;
+
{ .mfi
- nop.m 999
-(p14) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
- nop.i 999
+ nop.m 999
+ fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p11) fsub.s1 FR_A_hi = FR_A_hi,FR_sigma_A
- nop.i 999 ;;
+ nop.m 999
+ fadd.s3 FR_A = FR_B_hi,FR_C_lo
+ nop.i 999
}
+;;
+
+{ .mfi
+ nop.m 999
+ fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
+ nop.i 999
+}
+;;
+
{ .mfi
- nop.m 999
+ nop.m 999
+ fsub.s1 FR_Tmp_C = FR_Tmp_C,FR_Tmp2_C
+ nop.i 999
+}
+;;
+
// *******************
// Step 2. Get N and f
// *******************
@@ -1215,168 +1238,213 @@ alloc r34 = ar.pfs,2,34,0,0
// A := fadd.fpsr3( B_hi, C_lo )
// B := max( B_hi, C_lo )
// b := min( B_hi, C_lo )
-(p0) fadd.s3 FR_A = FR_B_hi,FR_C_lo
- nop.i 999
-}
{ .mfi
- nop.m 999
-(p10) fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
- nop.i 999 ;;
+ nop.m 999
+ fmax.s1 FR_B = FR_B_hi,FR_C_lo
+ nop.i 999
}
+;;
+
+// We use a right-shift trick to get the integer part of A into the rightmost
+// bits of the significand by adding 1.1000..00 * 2^63. This operation is good
+// if |A| < 2^61, which it is in this case. We are doing this to save a few
+// cycles over using fcvt.fx followed by fnorm. The second step of the trick
+// is to subtract the same constant to float the rounded integer into a fp reg.
+
{ .mfi
- nop.m 999
-(p0) fsub.s1 FR_Tmp_C = FR_Tmp_C,FR_Tmp2_C
- nop.i 999 ;;
+ nop.m 999
+// N := round_to_nearest_integer_value( A );
+ fma.s1 FR_N_fix = FR_A, f1, FR_RSHF
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fmax.s1 FR_B = FR_B_hi,FR_C_lo
- nop.i 999 ;;
+ nop.m 999
+ fmin.s1 FR_b = FR_B_hi,FR_C_lo
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fmin.s1 FR_b = FR_B_hi,FR_C_lo
- nop.i 999
+ nop.m 999
+// C_hi := C_hi - Tmp_C ...0 <= C_hi < 2^7
+ fsub.s1 FR_C_hi = FR_C_hi,FR_Tmp_C
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p11) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
- nop.i 999 ;;
+ nop.m 999
+// a := (B - A) + b: Exact - note that a is either 0 or 2^(-64).
+ fsub.s1 FR_a = FR_B,FR_A
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-// N := round_to_nearest_integer_value( A );
-(p0) fcvt.fx.s1 FR_N = FR_A
- nop.i 999 ;;
+ nop.m 999
+ fms.s1 FR_N = FR_N_fix, f1, FR_RSHF
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-// C_hi := C_hi - Tmp_C ...0 <= C_hi < 2^7
-(p0) fsub.s1 FR_C_hi = FR_C_hi,FR_Tmp_C
- nop.i 999 ;;
+ nop.m 999
+ fadd.s1 FR_a = FR_a,FR_b
+ nop.i 999
}
+;;
+
+// f := A - N; Exact because lsb(A) >= 2^(-64) and |f| <= 1/2.
+// N := convert to integer format( C_hi + N );
+// M := P_0 * x_lo;
+// N := N + M;
{ .mfi
- nop.m 999
-// a := (B - A) + b: Exact - note that a is either 0 or 2^(-64).
-(p0) fsub.s1 FR_a = FR_B,FR_A
- nop.i 999 ;;
+ nop.m 999
+ fsub.s1 FR_f = FR_A,FR_N
+ nop.i 999
}
{ .mfi
- nop.m 999
-// f := A - N; Exact because lsb(A) >= 2^(-64) and |f| <= 1/2.
-(p0) fnorm.s1 FR_N = FR_N
- nop.i 999
+ nop.m 999
+ fadd.s1 FR_N = FR_N,FR_C_hi
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fadd.s1 FR_a = FR_a,FR_b
- nop.i 999 ;;
+ nop.m 999
+(p9) fsub.s1 FR_D_hi = f0, FR_D_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fsub.s1 FR_f = FR_A,FR_N
- nop.i 999
+ nop.m 999
+(p9) fsub.s1 FR_D_lo = f0, FR_D_lo
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-// N := convert to integer format( C_hi + N );
-// M := P_0 * x_lo;
-// N := N + M;
-(p0) fadd.s1 FR_N = FR_N,FR_C_hi
- nop.i 999 ;;
+ nop.m 999
+ fadd.s1 FR_g = FR_A_hi,FR_B_lo // For Case 1, g=A_hi+B_lo
+ nop.i 999
}
{ .mfi
- nop.m 999
-// f = f + a Exact because a is 0 or 2^(-64);
-// the msb of the sum is <= 1/2 and lsb >= 2^(-64).
-(p0) fadd.s1 FR_f = FR_f,FR_a
- nop.i 999
+ nop.m 999
+ fadd.s3 FR_A = FR_A_hi,FR_B_lo // For Case 2, A=A_hi+B_lo w/ sf3
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-//
-// Create 2**(-33)
-//
-(p0) fcvt.fx.s1 FR_N = FR_N
- nop.i 999 ;;
+ mov GR_Temp = 0x0FFCD // For Case 2, exponent of 2^-50
+ fmax.s1 FR_B = FR_A_hi,FR_B_lo // For Case 2, B=max(A_hi,B_lo)
+ nop.i 999
}
+;;
+
+// f = f + a Exact because a is 0 or 2^(-64);
+// the msb of the sum is <= 1/2 and lsb >= 2^(-64).
{ .mfi
- nop.m 999
-(p0) fabs FR_f_abs = FR_f
- nop.i 999 ;;
+ setf.exp FR_TWOM50 = GR_Temp // For Case 2, form 2^-50
+ fcvt.fx.s1 FR_N = FR_N
+ nop.i 999
}
{ .mfi
-(p0) getf.sig GR_N = FR_N
- nop.f 999
- nop.i 999 ;;
+ nop.m 999
+ fadd.s1 FR_f = FR_f,FR_a
+ nop.i 999
}
-{ .mii
- nop.m 999
- nop.i 999 ;;
-(p0) add GR_N = GR_N,GR_M ;;
+;;
+
+{ .mfi
+ nop.m 999
+ fmin.s1 FR_b = FR_A_hi,FR_B_lo // For Case 2, b=min(A_hi,B_lo)
+ nop.i 999
}
-// If sgn_x == 1 (that is original x was negative)
-// N := 2^10 - N
-// this maintains N to be non-negative, but still
-// equivalent to the (negated N) mod 4.
-// End If
-{ .mii
-(p12) sub GR_N = GR_Temp,GR_N
-(p0) cmp.eq.unc p12,p9 = 0x0,GR_sgn_x ;;
- nop.i 999
+;;
+
+{ .mfi
+ nop.m 999
+ fsub.s1 FR_a = FR_B,FR_A // For Case 2, a=B-A
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fcmp.ge.unc.s1 p13, p10 = FR_f_abs,FR_TWOM33
- nop.i 999 ;;
+ nop.m 999
+ fadd.s1 FR_s_hi = FR_f,FR_g // For Case 1, s_hi=f+g
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p9) fsub.s1 FR_D_hi = f0, FR_D_hi
- nop.i 999 ;;
+ nop.m 999
+ fadd.s1 FR_f_hi = FR_A,FR_f // For Case 2, f_hi=A+f
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p10) fadd.s3 FR_A = FR_A_hi,FR_B_lo
- nop.i 999
+ nop.m 999
+ fabs FR_f_abs = FR_f
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p13) fadd.s1 FR_g = FR_A_hi,FR_B_lo
- nop.i 999 ;;
+ getf.sig GR_N = FR_N
+ fsetc.s3 0x7F,0x40 // Reset sf3 to user settings + td
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p10) fmax.s1 FR_B = FR_A_hi,FR_B_lo
- nop.i 999
+ nop.m 999
+ fsub.s1 FR_s_lo = FR_f,FR_s_hi // For Case 1, s_lo=f-s_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p9) fsub.s1 FR_D_lo = f0, FR_D_lo
- nop.i 999 ;;
+ nop.m 999
+ fsub.s1 FR_f_lo = FR_f,FR_f_hi // For Case 2, f_lo=f-f_hi
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p10) fmin.s1 FR_b = FR_A_hi,FR_B_lo
- nop.i 999 ;;
+ nop.m 999
+ fmpy.s1 FR_r_hi = FR_s_hi,FR_D_hi // For Case 1, r_hi=s_hi*D_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fsetc.s3 0x7F,0x40
- nop.i 999
+ nop.m 999
+ fadd.s1 FR_a = FR_a,FR_b // For Case 2, a=a+b
+ nop.i 999
}
-{ .mlx
- nop.m 999
-(p10) movl GR_Temp = 0x000000000000FFCD ;;
+;;
+
+
+// If sgn_x == 1 (that is original x was negative)
+// N := 2^10 - N
+// this maintains N to be non-negative, but still
+// equivalent to the (negated N) mod 4.
+// End If
+{ .mfi
+ add GR_N = GR_N,GR_M
+ fcmp.ge.s1 p13, p10 = FR_f_abs,FR_TWOM33
+ mov GR_Temp = 0x00400
}
-{ .mmf
- nop.m 999
-(p10) setf.exp FR_TWOM50 = GR_Temp
-(p10) fadd.s1 FR_f_hi = FR_A,FR_f ;;
+;;
+
+{ .mfi
+(p9) sub GR_N = GR_Temp,GR_N
+ fadd.s1 FR_s_lo = FR_s_lo,FR_g // For Case 1, s_lo=s_lo+g
+ nop.i 999
}
{ .mfi
- nop.m 999
-// a := (B - A) + b Exact.
+ nop.m 999
+ fadd.s1 FR_f_lo = FR_f_lo,FR_A // For Case 2, f_lo=f_lo+A
+ nop.i 999
+}
+;;
+
+// a := (B - A) + b Exact.
// Note that a is either 0 or 2^(-128).
// f_hi := A + f;
// f_lo := (f - f_hi) + A
@@ -1387,68 +1455,32 @@ alloc r34 = ar.pfs,2,34,0,0
// exact. If f = -2^(-64), then A + f is exact. Hence
// f-f_hi is -A exactly, giving f_lo = 0.
// f_lo := f_lo + a;
-(p10) fsub.s1 FR_a = FR_B,FR_A
- nop.i 999
-}
-{ .mfi
- nop.m 999
-(p13) fadd.s1 FR_s_hi = FR_f,FR_g
- nop.i 999 ;;
-}
-{ .mlx
- nop.m 999
+
// If |f| >= 2^(-33)
// Case 1
// CASE := 1
// g := A_hi + B_lo;
// s_hi := f + g;
// s_lo := (f - s_hi) + g;
-(p13) movl GR_CASE = 0x1 ;;
-}
-{ .mlx
- nop.m 999
// Else
// Case 2
// CASE := 2
// A := fadd.fpsr3( A_hi, B_lo )
// B := max( A_hi, B_lo )
// b := min( A_hi, B_lo )
-(p10) movl GR_CASE = 0x2
-}
-{ .mfi
- nop.m 999
-(p10) fsub.s1 FR_f_lo = FR_f,FR_f_hi
- nop.i 999 ;;
-}
-{ .mfi
- nop.m 999
-(p10) fadd.s1 FR_a = FR_a,FR_b
- nop.i 999
-}
-{ .mfi
- nop.m 999
-(p13) fsub.s1 FR_s_lo = FR_f,FR_s_hi
- nop.i 999 ;;
-}
-{ .mfi
- nop.m 999
-(p13) fadd.s1 FR_s_lo = FR_s_lo,FR_g
- nop.i 999 ;;
-}
+
{ .mfi
- nop.m 999
-(p10) fcmp.ge.unc.s1 p14, p11 = FR_f_abs,FR_TWOM50
- nop.i 999 ;;
+ nop.m 999
+(p10) fcmp.ge.unc.s1 p14, p11 = FR_f_abs,FR_TWOM50
+ nop.i 999
}
{ .mfi
- nop.m 999
-//
-// Create 2**(-50)
-(p10) fadd.s1 FR_f_lo = FR_f_lo,FR_A
- nop.i 999 ;;
+ nop.m 999
+(p13) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi //For Case 1, r_lo=s_hi*D_hi+r_hi
+ nop.i 999
}
-{ .mfi
- nop.m 999
+;;
+
// If |f| >= 2^(-50) then
// s_hi := f_hi;
// s_lo := f_lo;
@@ -1457,84 +1489,90 @@ alloc r34 = ar.pfs,2,34,0,0
// s_hi := f_hi + f_lo
// s_lo := (f_hi - s_hi) + f_lo
// End If
-(p14) mov FR_s_hi = FR_f_hi
- nop.i 999 ;;
+{ .mfi
+ nop.m 999
+(p14) mov FR_s_hi = FR_f_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p10) fadd.s1 FR_f_lo = FR_f_lo,FR_a
- nop.i 999 ;;
+ nop.m 999
+(p10) fadd.s1 FR_f_lo = FR_f_lo,FR_a
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p14) mov FR_s_lo = FR_f_lo
- nop.i 999
+ nop.m 999
+(p14) mov FR_s_lo = FR_f_lo
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p11) fadd.s1 FR_f_lo = FR_f_lo,FR_A_lo
- nop.i 999 ;;
+ nop.m 999
+(p11) fadd.s1 FR_f_lo = FR_f_lo,FR_A_lo
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p11) fma.s1 FR_f_lo = FR_X,FR_p_4,FR_f_lo
- nop.i 999 ;;
+ nop.m 999
+(p11) fma.s1 FR_f_lo = FR_X,FR_p_4,FR_f_lo
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p11) fadd.s1 FR_s_hi = FR_f_hi,FR_f_lo
- nop.i 999 ;;
+ nop.m 999
+(p13) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo //For Case 1, r_lo=s_hi*D_lo+r_lo
+ nop.i 999
}
{ .mfi
- nop.m 999
+ nop.m 999
+(p11) fadd.s1 FR_s_hi = FR_f_hi,FR_f_lo
+ nop.i 999
+}
+;;
+
// r_hi := s_hi*D_hi
// r_lo := s_hi*D_hi - r_hi with fma
// r_lo := (s_hi*D_lo + r_lo) + s_lo*D_hi
-(p0) fmpy.s1 FR_r_hi = FR_s_hi,FR_D_hi
- nop.i 999
-}
{ .mfi
- nop.m 999
-(p11) fsub.s1 FR_s_lo = FR_f_hi,FR_s_hi
- nop.i 999 ;;
+ nop.m 999
+(p10) fmpy.s1 FR_r_hi = FR_s_hi,FR_D_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi
- nop.i 999
+ nop.m 999
+(p11) fsub.s1 FR_s_lo = FR_f_hi,FR_s_hi
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p11) fadd.s1 FR_s_lo = FR_s_lo,FR_f_lo
- nop.i 999 ;;
-}
-{ .mmi
- nop.m 999 ;;
-// Return N, r_hi, r_lo
-// We do not return CASE
-(p0) stfe [GR_Address_of_Outputs] = FR_r_hi,16
- nop.i 999 ;;
+ nop.m 999
+(p10) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi
+ nop.i 999
}
{ .mfi
- nop.m 999
-(p0) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo
- nop.i 999 ;;
+ nop.m 999
+(p11) fadd.s1 FR_s_lo = FR_s_lo,FR_f_lo
+ nop.i 999
}
+;;
+
{ .mfi
- nop.m 999
-(p0) fma.s1 FR_r_lo = FR_s_lo,FR_D_hi,FR_r_lo
- nop.i 999 ;;
-}
-{ .mmi
- nop.m 999 ;;
-(p0) stfe [GR_Address_of_Outputs] = FR_r_lo,-16
- nop.i 999
+ nop.m 999
+(p10) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo
+ nop.i 999
}
-{ .mib
- nop.m 999
- nop.i 999
-(p0) br.ret.sptk b0 ;;
+;;
+
+// Return N, r_hi, r_lo
+// We do not return CASE
+{ .mfb
+ nop.m 999
+ fma.s1 FR_r_lo = FR_s_lo,FR_D_hi,FR_r_lo
+ br.ret.sptk b0
}
+;;
-.endp __libm_pi_by_2_reduce
-ASM_SIZE_DIRECTIVE(__libm_pi_by_2_reduce)
+.endp __libm_pi_by_2_reduce#
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-28 14:07:41 +0000