about summary refs log tree commit diff
path: root/sysdeps/ia64/fpu/libm_reduce.S
diff options
context:
space:
mode:
authorUlrich Drepper <drepper@redhat.com>2004-12-22 20:10:10 +0000
committerUlrich Drepper <drepper@redhat.com>2004-12-22 20:10:10 +0000
commita334319f6530564d22e775935d9c91663623a1b4 (patch)
treeb5877475619e4c938e98757d518bb1e9cbead751 /sysdeps/ia64/fpu/libm_reduce.S
parent0ecb606cb6cf65de1d9fc8a919bceb4be476c602 (diff)
downloadglibc-a334319f6530564d22e775935d9c91663623a1b4.tar.gz
glibc-a334319f6530564d22e775935d9c91663623a1b4.tar.xz
glibc-a334319f6530564d22e775935d9c91663623a1b4.zip
(CFLAGS-tst-align.c): Add -mpreferred-stack-boundary=4.
Diffstat (limited to 'sysdeps/ia64/fpu/libm_reduce.S')
-rw-r--r--sysdeps/ia64/fpu/libm_reduce.S1492
1 files changed, 727 insertions, 765 deletions
diff --git a/sysdeps/ia64/fpu/libm_reduce.S b/sysdeps/ia64/fpu/libm_reduce.S
index 8bdf91d6de..1c7f4e1e88 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 - 2003, Intel Corporation
+// Copyright (C) 2000, 2001, Intel Corporation
 // All rights reserved.
-//
-// Contributed 2000 by the Intel Numerics Group, Intel Corporation
+// 
+// Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
+// and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
 //
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are
@@ -20,310 +20,304 @@
 // * 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://www.intel.com/software/products/opensource/libraries/num.htm.
+// problem reports or change requests be submitted to it directly at 
+// http://developer.intel.com/opensource.
 //
-// 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
+// History:  02/02/00 Initial Version
 //
-//*********************************************************************
-//*********************************************************************
+// *********************************************************************
+// *********************************************************************
 //
 // 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:
-//      f8  = Input x, return value r
-//      f9  = return value c
-//      f32-f70
+//    Floating-Point Registers: 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 N and r_hi, r_lo where
-//
+// |x| >= 2^63, this routine returns CASE, 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  |
-//                    ---------   ---------   ---------
-//
-//                                            ---------
-//             X                              |   X   |
+// 	             |  P_1  |   |  P_2  |   |  P_3  |
+// 	             ---------   ---------   ---------
+// 
 //                                            ---------
+// 	      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)
@@ -331,30 +325,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  |
@@ -363,108 +357,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.
@@ -474,9 +468,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;
@@ -485,111 +479,117 @@
 //    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  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
+// 
+// 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 
 FR_f_lo    = f67
-FR_N_fix   = f68
-FR_C_hi    = f69
-FR_C_lo    = f70
+FR_r_lo    = f68 
+FR_C_hi    = f69 
+FR_C_lo    = f70 
 
 GR_N       = r8
-GR_Exp_x   = r36
-GR_Temp    = r37
-GR_BIASL63 = r38
+GR_Address_of_Input  = r32 
+GR_Address_of_Outputs = r33 
+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
 
-LOCAL_OBJECT_START(Constants_Bits_of_2_by_pi)
+#ifdef _LIBC
+.rodata
+#else
+.data
+#endif
+
+Constants_Bits_of_2_by_pi:
+ASM_TYPE_DIRECTIVE(Constants_Bits_of_2_by_pi,@object)
 data8 0x0000000000000000,0xA2F9836E4E441529
 data8 0xFC2757D1F534DDC0,0xDB6295993C439041
 data8 0xFE5163ABDEBBC561,0xB7246E3A424DD2E0
@@ -721,33 +721,34 @@ data8 0xB5D6DF8261DD9602,0x36169F3AC4A1A283
 data8 0x6DED727A8D39A9B8,0x825C326B5B2746ED
 data8 0x34007700D255F4FC,0x4D59018071E0E13F
 data8 0x89B295F364A8F1AE,0xA74B38FC4CEAB2BB
-LOCAL_OBJECT_END(Constants_Bits_of_2_by_pi)
+ASM_SIZE_DIRECTIVE(Constants_Bits_of_2_by_pi)
 
-LOCAL_OBJECT_START(Constants_Bits_of_pi_by_2)
-data8 0xC90FDAA22168C234,0x00003FFF
-data8 0xC4C6628B80DC1CD1,0x00003FBF
-LOCAL_OBJECT_END(Constants_Bits_of_pi_by_2)
+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)
 
 .section .text
-.global __libm_pi_by_2_reduce#
 .proc __libm_pi_by_2_reduce#
-.align 32
+.global __libm_pi_by_2_reduce#
+.align 64 
 
-__libm_pi_by_2_reduce:
+__libm_pi_by_2_reduce: 
 
-//    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
+//    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
 
-{ .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
+{ .mmf
+alloc  r34 = ar.pfs,2,34,0,0
+(p0)  ldfe  FR_X = [GR_Address_of_Input]
+(p0)  fsetc.s3 0x00,0x7F ;;
 }
-{ .mfi
-      addl           GR_BASE   = @ltoff(Constants_Bits_of_2_by_pi#), gp
-      nop.f 999
-      mov GR_BIASL63 = 0x1003E
+{ .mlx
+	nop.m 999
+(p0)  movl GR_BIASL63 = 0x1003E
 }
 ;;
 
@@ -764,61 +765,73 @@ __libm_pi_by_2_reduce:
 //    Address_BASE = shladd(SEGMENT,3) + BASE
 
 
+
 { .mmi
-      getf.exp GR_Exp_x = FR_input_X
-      ld8 GR_BASE = [GR_BASE]
-      mov GR_TEMP5 = 0x0FFFE
+      nop.m 999
+(p0)  addl           GR_BASE   = @ltoff(Constants_Bits_of_2_by_pi#), gp
+      nop.i 999
 }
 ;;
 
-//    Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_A := 2^(-65).
 { .mmi
-      getf.sig GR_x_lo = FR_input_X
-      mov GR_TEMP6 = 0x0FFBE
+      ld8 GR_BASE = [GR_BASE]
+      nop.m 999
       nop.i 999
 }
 ;;
 
-//    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
+
+{ .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
 //    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.
-{ .mmi
-      setf.exp FR_sigma_B = GR_TEMP5
-      setf.exp FR_sigma_A = GR_TEMP6
-      extr.u GR_M = GR_Exp_x,0,17
+{ .mmf
+(p0)  getf.exp GR_Exp_x = FR_X
+(p0)  getf.sig GR_x_lo = FR_X
+(p0)  fabs FR_X = FR_X ;;
 }
-;;
-
 { .mii
-      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
+(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
 }
-;;
-
-{ .mii
-      and  GR_LENGTH1 = 0x3F,GR_START
-      shr.u  GR_SEGMENT = GR_START,6
-      nop.i 999
+{ .mmi
+	nop.m 999 ;;
+(p0)  and  GR_LENGTH1 = 0x3F,GR_START
+(p0)  shr.u  GR_SEGMENT = GR_START,6
 }
-;;
-
 { .mmi
-      shladd GR_BASE = GR_SEGMENT,3,GR_BASE
-      sub  GR_LENGTH2 = 0x40,GR_LENGTH1
-      cmp.le p6,p7 = 0x2,GR_LENGTH1
+	nop.m 999 ;;
+(p0)  add  GR_SEGMENT = 0x1,GR_SEGMENT
+(p0)  sub  GR_LENGTH2 = 0x40,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
@@ -836,13 +849,13 @@ __libm_pi_by_2_reduce:
 //    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
-{ .mfi
-      ld8 GR_A = [GR_BASE],8
-      fabs FR_X = FR_input_X
-(p7)  cmp.eq.unc p8,p9 = 0x1,GR_LENGTH1
+{ .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 ;;
 }
-;;
-
+{ .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
@@ -853,35 +866,31 @@ __libm_pi_by_2_reduce:
 //    A, B, C, D, and E look like    | length1 | length2   |
 //                                    ---------------------
 //                                       hi        lo
-{ .mlx
-      ld8 GR_B = [GR_BASE],8
-      movl GR_rshf = 0x43e8000000000000   // 1.10000 2^63 for right shift N_fix
+(p0)  ld8 GR_A = [GR_BASE],8
+(p0)  extr.u GR_sgn_x = GR_Exp_x,17,1 ;;
 }
-;;
-
-{ .mmi
-      ld8 GR_C = [GR_BASE],8
-      nop.m 999
-(p8)  extr.u GR_Temp = GR_A,63,1
+{ .mmf
+	nop.m 999
+(p0)  ld8 GR_B = [GR_BASE],8
+(p0)  fmerge.se FR_X = FR_sigma_B,FR_X ;;
 }
-;;
-
+{ .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
+}
+{ .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.
-{ .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
+(p6)     shr.u GR_P_0 = GR_A,GR_LENGTH2 ;;
+(p0)  shl GR_TEMP2 = GR_B,GR_LENGTH1
 }
-;;
-
 { .mii
-      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
+(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
 }
-;;
-
 //    Else
 //       Load 16 bit of ASUB from (Base_Address_of_A - 2)
 //       P_0 = ASUB & 0x3
@@ -891,56 +900,43 @@ __libm_pi_by_2_reduce:
 //          Deposit element 63 from Ahi and place in element 0 of P_0.
 //       Endif
 //    Endif
-
 { .mii
 (p7)  ld2 GR_ASUB = [GR_BASE],8
-      shl GR_TEMP3 = GR_C,GR_LENGTH1   // MM instruction
-      shr.u GR_P_2 = GR_C,GR_LENGTH2   // MM instruction
+(p0)  shl GR_TEMP3 = GR_C,GR_LENGTH1 ;;
+(p0)  shl GR_TEMP4 = GR_D,GR_LENGTH1
 }
-;;
-
 { .mii
-      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
+	nop.m 999
+(p0)  shr.u GR_P_3 = GR_D,GR_LENGTH2 ;;
+(p0)  shr.u GR_P_4 = GR_E,GR_LENGTH2
 }
-;;
-
-{ .mmi
+{ .mii
 (p7)  and GR_P_0 = 0x03,GR_ASUB
-(p6)  and GR_P_0 = 0x03,GR_P_0
-      shr.u GR_P_4 = GR_E,GR_LENGTH2   // MM instruction
+(p6)     and GR_P_0 = 0x03,GR_P_0 ;;
+(p0)  or GR_P_1 = GR_P_1,GR_TEMP1
 }
-;;
-
 { .mmi
-      nop.m 999
-      or GR_P_1 = GR_P_1,GR_TEMP1
-(p8)  and GR_P_0 = 0x1,GR_P_0
+(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 ;;
 }
-;;
-
-{ .mmi
-      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
+{ .mii
+	nop.m 999
+(p0)  or GR_P_3 = GR_P_3,GR_TEMP3
+(p8)  or GR_P_0 = GR_P_0,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
+(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 ;;
 }
-;;
-
 { .mmi
-      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
+	nop.m 999 ;;
+(p0)  setf.sig FR_p_3 = GR_P_3
+(p0)  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;
@@ -958,18 +954,18 @@ __libm_pi_by_2_reduce:
 //             |  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  |
@@ -981,55 +977,52 @@ __libm_pi_by_2_reduce:
 //    and exponent range. Unless an explicit FPSR is given,
 //    round-to-nearest with widest precision and exponent range is
 //    used.
-{ .mmi
-      setf.sig FR_p_4 = GR_P_4
-      mov GR_TEMP1 = 0x0FFBF
-      nop.i 999
+(p0)  movl GR_TEMP1 = 0x000000000000FFBF
 }
-;;
-
 { .mmi
-      setf.exp FR_ScaleP2 = GR_TEMP1
-      mov GR_TEMP2 = 0x0FF7F
-      nop.i 999
+	nop.m 999 ;;
+(p0)  setf.exp FR_ScaleP2 = GR_TEMP1
+	nop.i 999
+}
+{ .mlx
+	nop.m 999
+(p0)  movl GR_TEMP4 = 0x000000000001003E
 }
-;;
-
 { .mmi
-      setf.exp FR_ScaleP3 = GR_TEMP2
-      mov GR_TEMP4 = 0x1003E
-      nop.i 999
+	nop.m 999 ;;
+(p0)  setf.exp FR_sigma_C = GR_TEMP4
+	nop.i 999
+}
+{ .mlx
+	nop.m 999
+(p0)  movl GR_TEMP2 = 0x000000000000FF7F ;;
 }
-;;
-
 { .mmf
-      setf.exp FR_sigma_C = GR_TEMP4
-      mov GR_Temp = 0x0FFDE
-      fcvt.xuf.s1 FR_p_1 = FR_p_1
+	nop.m 999
+(p0)  setf.exp FR_ScaleP3 = GR_TEMP2
+(p0)  fcvt.xuf.s1 FR_p_1 = FR_p_1 ;;
 }
-;;
-
 { .mfi
-      setf.exp FR_TWOM33 = GR_Temp
-      fcvt.xuf.s1 FR_p_2 = FR_p_2
-      nop.i 999
+	nop.m 999
+(p0)  fcvt.xuf.s1 FR_p_2 = FR_p_2
+	nop.i 999
 }
-;;
-
-{ .mfi
-      nop.m 999
-      fcvt.xuf.s1 FR_p_3 = FR_p_3
-      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 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fcvt.xuf.s1 FR_p_4 = FR_p_4
-      nop.i 999
+	nop.m 999
+(p0)  fcvt.xuf.s1 FR_p_4 = FR_p_4
+	nop.i 999 ;;
 }
-;;
-
+{ .mfi
+	nop.m 999
 //    Tmp_C := fmpy.fpsr3( x, p_1 );
 //    Tmp_B := fmpy.fpsr3( x, p_2 );
 //    Tmp_A := fmpy.fpsr3( x, p_3 );
@@ -1055,62 +1048,55 @@ __libm_pi_by_2_reduce:
 //      Exact, regardless ...of rounding direction
 //      A_lo := x*p_3 - A_hi ...fma, exact
 //    Endif
-{ .mfi
-      nop.m 999
-      fmpy.s3 FR_Tmp_C = FR_X,FR_p_1
-      nop.i 999
+(p0)  fmpy.s3 FR_Tmp_C = FR_X,FR_p_1
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      mov GR_TEMP3 = 0x0FF3F
-      fmpy.s1 FR_p_2 = FR_p_2,FR_ScaleP2
-      nop.i 999
+	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
+}
+{ .mlx
+	nop.m 999
+(p0)  movl GR_TEMP3 = 0x000000000000FF3F ;;
 }
-;;
-
 { .mmf
-      setf.exp FR_ScaleP4 = GR_TEMP3
-      mov GR_TEMP4 = 0x10045
-      fmpy.s1 FR_p_3 = FR_p_3,FR_ScaleP3
+	nop.m 999
+(p0)  setf.exp FR_ScaleP4 = GR_TEMP3
+(p0)  fmpy.s1 FR_p_3 = FR_p_3,FR_ScaleP3 ;;
 }
-;;
-
-{ .mfi
-      nop.m 999
-      fadd.s3 FR_C_hi = FR_sigma_C,FR_Tmp_C   // For Tmp_C < sigma_C case
-      nop.i 999
+{ .mlx
+	nop.m 999
+(p0)  movl GR_TEMP4 = 0x0000000000010045 ;;
 }
-;;
-
 { .mmf
-      setf.exp FR_Tmp2_C = GR_TEMP4
-      nop.m 999
-      fmpy.s3 FR_Tmp_B = FR_X,FR_p_2
+	nop.m 999
+(p0)  setf.exp FR_Tmp2_C = GR_TEMP4
+(p0)  fmpy.s3 FR_Tmp_B = FR_X,FR_p_2 ;;
 }
-;;
-
 { .mfi
-      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
+	nop.m 999
+(p0)  fcmp.ge.unc.s1 p12,  p9 = FR_Tmp_C,FR_sigma_C
+	nop.i 999 ;;
 }
 { .mfi
-      nop.m 999
-      fmpy.s3 FR_Tmp_A = FR_X,FR_p_3
-      nop.i 99
+	nop.m 999
+(p0)  fmpy.s3 FR_Tmp_A = FR_X,FR_p_3
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      ld8 GR_BASE = [GR_BASE]
+	nop.m 999
 (p12) mov FR_C_hi = FR_Tmp_C
-      nop.i 999
+	nop.i 999 ;;
 }
 { .mfi
-      nop.m 999
-(p9)  fsub.s1 FR_C_hi = FR_C_hi,FR_sigma_C
-      nop.i 999
+(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
 }
 ;;
 
@@ -1128,106 +1114,97 @@ __libm_pi_by_2_reduce:
 //      Load from table
 //   End If
 
-{ .mfi
-      nop.m 999
-      fmpy.s1 FR_p_4 = FR_p_4,FR_ScaleP4
-      nop.i 999
-}
-{ .mfi
+
+{ .mmi
+      ld8 GR_BASE = [GR_BASE]
       nop.m 999
-      fadd.s3 FR_B_hi = FR_sigma_B,FR_Tmp_B     // For Tmp_B < sigma_B case
       nop.i 999
 }
 ;;
 
+
 { .mfi
-      nop.m 999
-      fadd.s3 FR_A_hi = FR_sigma_A,FR_Tmp_A     // For Tmp_A < sigma_A case
-      nop.i 999
+(p0) ldfe FR_D_hi = [GR_BASE],16
+(p0)  fmpy.s1 FR_p_4 = FR_p_4,FR_ScaleP4
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fcmp.ge.s1 p13, p10 = FR_Tmp_B,FR_sigma_B
-      nop.i 999
+(p0) ldfe FR_D_lo = [GR_BASE],0
+(p0)  fcmp.ge.unc.s1 p13, p10 = FR_Tmp_B,FR_sigma_B
+	nop.i 999 ;;
 }
 { .mfi
-      nop.m 999
-      fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
-      nop.i 999
+	nop.m 999
+(p13) mov FR_B_hi = FR_Tmp_B
+	nop.i 999
 }
-;;
-
 { .mfi
-      ldfe FR_D_hi = [GR_BASE],16
-      fcmp.ge.s1 p14, p11 = FR_Tmp_A,FR_sigma_A
-      nop.i 999
+	nop.m 999
+(p12) fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      ldfe FR_D_lo = [GR_BASE]
-(p13) mov FR_B_hi = FR_Tmp_B
-      nop.i 999
+	nop.m 999
+(p10) fadd.s3 FR_B_hi = FR_sigma_B,FR_Tmp_B
+	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
+(p9)  fsub.s1 FR_C_hi = FR_C_hi,FR_sigma_C
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
+	nop.m 999
+(p0)  fcmp.ge.unc.s1 p14, p11 = FR_Tmp_A,FR_sigma_A
+	nop.i 999 ;;
+}
+{ .mfi
+	nop.m 999
 (p14) mov FR_A_hi = FR_Tmp_A
-      nop.i 999
+	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
+(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 ;;
+}
+{ .mfi
+	nop.m 999
+(p10) fsub.s1 FR_B_hi = FR_B_hi,FR_sigma_B
+	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
-{ .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
-      fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
-      nop.i 999
+(p0)  fadd.s3 FR_Tmp_C = FR_C_hi,FR_Tmp2_C
+	nop.i 999 ;;
 }
 { .mfi
-      nop.m 999
-      fadd.s3 FR_A = FR_B_hi,FR_C_lo
-      nop.i 999
+	nop.m 999
+(p14) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
+	nop.i 999
 }
-;;
-
 { .mfi
-      nop.m 999
-      fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
-      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
-      fsub.s1 FR_Tmp_C = FR_Tmp_C,FR_Tmp2_C
-      nop.i 999
-}
-;;
-
+	nop.m 999
 //    *******************
 //    Step 2. Get N and f
 //    *******************
@@ -1238,213 +1215,168 @@ __libm_pi_by_2_reduce:
 //    A := fadd.fpsr3( B_hi, C_lo )
 //    B := max( B_hi, C_lo )
 //    b := min( B_hi, C_lo )
-{ .mfi
-      nop.m 999
-      fmax.s1 FR_B = FR_B_hi,FR_C_lo
-      nop.i 999
+(p0)  fadd.s3 FR_A = 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
-//    N := round_to_nearest_integer_value( A );
-      fma.s1 FR_N_fix = FR_A, f1, FR_RSHF
-      nop.i 999
+	nop.m 999
+(p10) fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fmin.s1 FR_b = FR_B_hi,FR_C_lo
-      nop.i 999
+	nop.m 999
+(p0)  fsub.s1 FR_Tmp_C = FR_Tmp_C,FR_Tmp2_C
+	nop.i 999 ;;
 }
 { .mfi
-      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
+	nop.m 999
+(p0)  fmax.s1 FR_B = FR_B_hi,FR_C_lo
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      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
+	nop.m 999
+(p0)  fmin.s1 FR_b = FR_B_hi,FR_C_lo
+	nop.i 999
 }
-;;
-
 { .mfi
-      nop.m 999
-      fms.s1 FR_N = FR_N_fix, f1, FR_RSHF
-      nop.i 999
+	nop.m 999
+(p11) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fadd.s1 FR_a = FR_a,FR_b
-      nop.i 999
+	nop.m 999
+//    N := round_to_nearest_integer_value( A );
+(p0)  fcvt.fx.s1 FR_N = FR_A
+	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
-      fsub.s1 FR_f = FR_A,FR_N
-      nop.i 999
+	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 ;;
 }
 { .mfi
-      nop.m 999
-      fadd.s1 FR_N = FR_N,FR_C_hi
-      nop.i 999
+	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 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-(p9)  fsub.s1 FR_D_hi = f0, FR_D_hi
-      nop.i 999
+	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
 }
 { .mfi
-      nop.m 999
-(p9)  fsub.s1 FR_D_lo = f0, FR_D_lo
-      nop.i 999
+	nop.m 999
+(p0)  fadd.s1 FR_a = FR_a,FR_b
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fadd.s1 FR_g = FR_A_hi,FR_B_lo          // For Case 1, g=A_hi+B_lo
-      nop.i 999
+	nop.m 999
+(p0)  fsub.s1 FR_f = FR_A,FR_N
+	nop.i 999
 }
 { .mfi
-      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
+	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 ;;
 }
-;;
-
 { .mfi
-      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);
+	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).
-{ .mfi
-      setf.exp FR_TWOM50 = GR_Temp            // For Case 2, form 2^-50
-      fcvt.fx.s1 FR_N = FR_N
-      nop.i 999
+(p0)  fadd.s1 FR_f = FR_f,FR_a
+	nop.i 999
 }
 { .mfi
-      nop.m 999
-      fadd.s1 FR_f = FR_f,FR_a
-      nop.i 999
+	nop.m 999
+//
+//    Create 2**(-33)
+//
+(p0)  fcvt.fx.s1 FR_N = FR_N
+	nop.i 999 ;;
 }
-;;
-
 { .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
+	nop.m 999
+(p0)  fabs FR_f_abs = FR_f
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fsub.s1 FR_a = FR_B,FR_A                // For Case 2, a=B-A
-      nop.i 999
+(p0)  getf.sig GR_N = FR_N
+	nop.f 999
+	nop.i 999 ;;
 }
-;;
-
-{ .mfi
-      nop.m 999
-      fadd.s1 FR_s_hi = FR_f,FR_g             // For Case 1, s_hi=f+g
-      nop.i 999
+{ .mii
+	nop.m 999
+	nop.i 999 ;;
+(p0)  add GR_N = GR_N,GR_M ;;
 }
-{ .mfi
-      nop.m 999
-      fadd.s1 FR_f_hi = FR_A,FR_f             // For Case 2, f_hi=A+f
-      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
-      fabs FR_f_abs = FR_f
-      nop.i 999
+	nop.m 999
+(p0)  fcmp.ge.unc.s1 p13, p10 = FR_f_abs,FR_TWOM33
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      getf.sig GR_N = FR_N
-      fsetc.s3 0x7F,0x40                 // Reset sf3 to user settings + td
-      nop.i 999
+	nop.m 999
+(p9) fsub.s1 FR_D_hi = f0, FR_D_hi
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-      fsub.s1 FR_s_lo = FR_f,FR_s_hi          // For Case 1, s_lo=f-s_hi
-      nop.i 999
+	nop.m 999
+(p10)    fadd.s3 FR_A = FR_A_hi,FR_B_lo
+	nop.i 999
 }
 { .mfi
-      nop.m 999
-      fsub.s1 FR_f_lo = FR_f,FR_f_hi          // For Case 2, f_lo=f-f_hi
-      nop.i 999
+	nop.m 999
+(p13)    fadd.s1 FR_g = FR_A_hi,FR_B_lo
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      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
+	nop.m 999
+(p10)    fmax.s1 FR_B = FR_A_hi,FR_B_lo
+	nop.i 999
 }
 { .mfi
-      nop.m 999
-      fadd.s1 FR_a = FR_a,FR_b                // For Case 2, a=a+b
-      nop.i 999
+	nop.m 999
+(p9) fsub.s1 FR_D_lo = f0, FR_D_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
 { .mfi
-      add GR_N = GR_N,GR_M
-      fcmp.ge.s1 p13, p10 = FR_f_abs,FR_TWOM33
-      mov GR_Temp = 0x00400
+	nop.m 999
+(p10)    fmin.s1 FR_b = FR_A_hi,FR_B_lo
+	nop.i 999 ;;
 }
-;;
-
 { .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
+	nop.m 999
+(p0) fsetc.s3 0x7F,0x40
+	nop.i 999
 }
-{ .mfi
-      nop.m 999
-      fadd.s1 FR_f_lo = FR_f_lo,FR_A           // For Case 2, f_lo=f_lo+A
-      nop.i 999
+{ .mlx
+	nop.m 999
+(p10)    movl GR_Temp = 0x000000000000FFCD ;;
 }
-;;
-
-//       a := (B - A) + b      Exact.
+{ .mmf
+	nop.m 999
+(p10)    setf.exp FR_TWOM50 = GR_Temp
+(p10)    fadd.s1 FR_f_hi = FR_A,FR_f ;;
+}
+{ .mfi
+	nop.m 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
@@ -1455,32 +1387,68 @@ __libm_pi_by_2_reduce:
 //       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) fcmp.ge.unc.s1 p14, p11 = FR_f_abs,FR_TWOM50
-      nop.i 999
+	nop.m 999
+(p10)    fsub.s1 FR_f_lo = FR_f,FR_f_hi
+	nop.i 999 ;;
 }
 { .mfi
-      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
+	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 ;;
+}
+{ .mfi
+	nop.m 999
+//
+//       Create 2**(-50)
+(p10)    fadd.s1 FR_f_lo = FR_f_lo,FR_A
+	nop.i 999 ;;
+}
+{ .mfi
+	nop.m 999
 //       If |f| >= 2^(-50) then
 //          s_hi := f_hi;
 //          s_lo := f_lo;
@@ -1489,90 +1457,84 @@ __libm_pi_by_2_reduce:
 //          s_hi := f_hi + f_lo
 //          s_lo := (f_hi - s_hi) + f_lo
 //       End If
-{ .mfi
-      nop.m 999
-(p14) mov FR_s_hi = FR_f_hi
-      nop.i 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
-(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
+	nop.m 999
+(p11)  fadd.s1 FR_s_hi = FR_f_hi,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
 //   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
-(p10) fmpy.s1 FR_r_hi = FR_s_hi,FR_D_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) fsub.s1 FR_s_lo = FR_f_hi,FR_s_hi
-      nop.i 999
+	nop.m 999
+(p0) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi
+	nop.i 999
 }
-;;
-
 { .mfi
-      nop.m 999
-(p10) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi
-      nop.i 999
+	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 ;;
 }
 { .mfi
-      nop.m 999
-(p11) fadd.s1 FR_s_lo = FR_s_lo,FR_f_lo
-      nop.i 999
+	nop.m 999
+(p0) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo
+	nop.i 999 ;;
 }
-;;
-
 { .mfi
-      nop.m 999
-(p10) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo
-      nop.i 999
+	nop.m 999
+(p0) fma.s1 FR_r_lo = FR_s_lo,FR_D_hi,FR_r_lo
+	nop.i 999 ;;
 }
-;;
-
-//   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
+{ .mmi
+	nop.m 999 ;;
+(p0) stfe [GR_Address_of_Outputs] = FR_r_lo,-16
+	nop.i 999
+}
+{ .mib
+	nop.m 999
+	nop.i 999
+(p0) br.ret.sptk   b0 ;;
 }
-;;
 
-.endp __libm_pi_by_2_reduce#
+.endp __libm_pi_by_2_reduce
+ASM_SIZE_DIRECTIVE(__libm_pi_by_2_reduce)