2.5-18.1

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#