.file "libm_tan.s"// Copyright (C) 2000, 2001, 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.//// Redistribution and use in source and binary forms, with or without// modification, are permitted provided that the following conditions are// met://// * Redistributions of source code must retain the above copyright// notice, this list of conditions and the following disclaimer.//// * Redistributions in binary form must reproduce the above copyright// notice, this list of conditions and the following disclaimer in the// documentation and/or other materials provided with the distribution.//// * The name of Intel Corporation may not be used to endorse or promote// products derived from this software without specific prior written// permission.//// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR// 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 // 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. // // 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.//// *********************************************************************//// History:  // 02/02/00 Initial Version // 4/04/00  Unwind support added// 12/28/00 Fixed false invalid flags//// *********************************************************************//// Function:   tan(x) = tangent(x), for double precision x values//// *********************************************************************//// Accuracy:       Very accurate for double-precision values  //// *********************************************************************//// Resources Used:////    Floating-Point Registers: f8 (Input and Return Value)//                              f9-f15//                              f32-f112////    General Purpose Registers://      r32-r48//      r49-r50 (Used to pass arguments to pi_by_2 reduce routine)////    Predicate Registers:      p6-p15//// *********************************************************************//// IEEE Special Conditions:////    Denormal  fault raised on denormal inputs//    Overflow exceptions do not occur//    Underflow exceptions raised when appropriate for tan //    (No specialized error handling for this routine)//    Inexact raised when appropriate by algorithm////    tan(SNaN) = QNaN//    tan(QNaN) = QNaN//    tan(inf) = QNaN//    tan(+/-0) = +/-0//// *********************************************************************//// Mathematical Description//// We consider the computation of FPTAN of Arg. Now, given////      Arg = N pi/2  + alpha,          |alpha| <= pi/4,//// basic mathematical relationship shows that////      tan( Arg ) =  tan( alpha )     if N is even;//                 = -cot( alpha )      otherwise.//// The value of alpha is obtained by argument reduction and// represented by two working precision numbers r and c where////      alpha =  r  +  c     accurately.//// The reduction method is described in a previous write up.// The argument reduction scheme identifies 4 cases. For Cases 2// and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be// computed very easily by 2 or 3 terms of the Taylor series// expansion as follows://// Case 2:// -------////      tan(r + c) = r + c + r^3/3          ...accurately//        -cot(r + c) = -1/(r+c) + r/3          ...accurately//// Case 4:// -------////      tan(r + c) = r + c + r^3/3 + 2r^5/15     ...accurately//        -cot(r + c) = -1/(r+c) + r/3 + r^3/45     ...accurately////// The only cases left are Cases 1 and 3 of the argument reduction// procedure. These two cases will be merged since after the// argument is reduced in either cases, we have the reduced argument// represented as r + c and that the magnitude |r + c| is not small// enough to allow the usage of a very short approximation.//// The greatest challenge of this task is that the second terms of// the Taylor series for tan(r) and -cot(r)////      r + r^3/3 + 2 r^5/15 + ...//// and////      -1/r + r/3 + r^3/45 + ...//// are not very small when |r| is close to pi/4 and the rounding// errors will be a concern if simple polynomial accumulation is// used. When |r| < 2^(-2), however, the second terms will be small// enough (5 bits or so of right shift) that a normal Horner// recurrence suffices. Hence there are two cases that we consider// in the accurate computation of tan(r) and cot(r), |r| <= pi/4.//// Case small_r: |r| < 2^(-2)// --------------------------//// Since Arg = N pi/4 + r + c accurately, we have////      tan(Arg) =  tan(r+c)            for N even,//            = -cot(r+c)          otherwise.//// Here for this case, both tan(r) and -cot(r) can be approximated// by simple polynomials:////      tan(r) =    r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19//        -cot(r) = -1/r + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13//// accurately. Since |r| is relatively small, tan(r+c) and// -cot(r+c) can be accurately approximated by replacing r with// r+c only in the first two terms of the corresponding polynomials.//// Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to// almost 64 sig. bits, thus////      P1_1 (r+c)^3 =  P1_1 r^3 + c * r^2     accurately.//// Hence,////      tan(r+c) =    r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19//                     + c*(1 + r^2)////        -cot(r+c) = -1/(r+c) + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13//               + Q1_1*c////// Case normal_r: 2^(-2) <= |r| <= pi/4// ------------------------------------//// This case is more likely than the previous one if one considers// r to be uniformly distributed in [-pi/4 pi/4].//// The required calculation is either////      tan(r + c)  =  tan(r)  +  correction,  or//        -cot(r + c)  = -cot(r)  +  correction.//// Specifically,////      tan(r + c) =  tan(r) + c tan'(r)  + O(c^2)//              =  tan(r) + c sec^2(r) + O(c^2)//              =  tan(r) + c SEC_sq     ...accurately//                as long as SEC_sq approximates sec^2(r)//                to, say, 5 bits or so.//// Similarly,////        -cot(r + c) = -cot(r) - c cot'(r)  + O(c^2)//              = -cot(r) + c csc^2(r) + O(c^2)//              = -cot(r) + c CSC_sq     ...accurately//                as long as CSC_sq approximates csc^2(r)//                to, say, 5 bits or so.//// We therefore concentrate on accurately calculating tan(r) and// cot(r) for a working-precision number r, |r| <= pi/4 to within// 0.1% or so.//// We will employ a table-driven approach. Let////      r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63//        = sgn_r * ( B + x )//// where////      B = 2^k * 1.b_1 b_2 ... b_5 1//         x = |r| - B//// Now,//                   tan(B)  +   tan(x)//      tan( B + x ) =  ------------------------//                   1 -  tan(B)*tan(x)////               /                         \ //               |   tan(B)  +   tan(x)          |//      = tan(B) +  | ------------------------ - tan(B) |//               |     1 -  tan(B)*tan(x)          |//               \                         /////                 sec^2(B) * tan(x)//      = tan(B) + ------------------------//                 1 -  tan(B)*tan(x)////                (1/[sin(B)*cos(B)]) * tan(x)//      = tan(B) + --------------------------------//                      cot(B)  -  tan(x)////// Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are// calculated beforehand and stored in a table. Since////      |x| <= 2^k * 2^(-6)  <= 2^(-7)  (because k = -1, -2)//// a very short polynomial will be sufficient to approximate tan(x)// accurately. The details involved in computing the last expression// will be given in the next section on algorithm description.////// Now, we turn to the case where cot( B + x ) is needed.//////                   1 - tan(B)*tan(x)//      cot( B + x ) =  ------------------------//                   tan(B)  +  tan(x)////               /                           \ //               |   1 - tan(B)*tan(x)              |//      = cot(B) +  | ----------------------- - cot(B) |//               |     tan(B)  +  tan(x)            |//               \                           /////               [tan(B) + cot(B)] * tan(x)//      = cot(B) - ----------------------------//                   tan(B)  +  tan(x)////                (1/[sin(B)*cos(B)]) * tan(x)//      = cot(B) - --------------------------------//                      tan(B)  +  tan(x)////// Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that// are needed are the same set of values needed in the previous// case.//// Finally, we can put all the ingredients together as follows:////      Arg = N * pi/2 +  r + c          ...accurately////      tan(Arg) =  tan(r) + correction    if N is even;//            = -cot(r) + correction    otherwise.//// For Cases 2 and 4,////     Case 2://     tan(Arg) =  tan(r + c) = r + c + r^3/3           N even//              = -cot(r + c) = -1/(r+c) + r/3           N odd//     Case 4://     tan(Arg) =  tan(r + c) = r + c + r^3/3 + 2r^5/15  N even//              = -cot(r + c) = -1/(r+c) + r/3 + r^3/45  N odd////// For Cases 1 and 3,////     Case small_r: |r| < 2^(-2)////      tan(Arg) =  r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19//                     + c*(1 + r^2)               N even////                  = -1/(r+c) + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13//               + Q1_1*c                    N odd////     Case normal_r: 2^(-2) <= |r| <= pi/4////      tan(Arg) =  tan(r) + c * sec^2(r)     N even//               = -cot(r) + c * csc^2(r)     otherwise////     For N even,////      tan(Arg) = tan(r) + c*sec^2(r)//               = tan( sgn_r * (B+x) ) + c * sec^2(|r|)//                  = sgn_r * ( tan(B+x)  + sgn_r*c*sec^2(|r|) )//                  = sgn_r * ( tan(B+x)  + sgn_r*c*sec^2(B) )//// since B approximates |r| to 2^(-6) in relative accuracy.////                 /            (1/[sin(B)*cos(B)]) * tan(x)//    tan(Arg) = sgn_r * | tan(B) + --------------------------------//                 \                     cot(B)  -  tan(x)//                                        \ //                       + CORR  |//                                     /// where////    CORR = sgn_r*c*tan(B)*SC_inv(B);  SC_inv(B) = 1/(sin(B)*cos(B)).//// For N odd,////      tan(Arg) = -cot(r) + c*csc^2(r)//               = -cot( sgn_r * (B+x) ) + c * csc^2(|r|)//                  = sgn_r * ( -cot(B+x)  + sgn_r*c*csc^2(|r|) )//                  = sgn_r * ( -cot(B+x)  + sgn_r*c*csc^2(B) )//// since B approximates |r| to 2^(-6) in relative accuracy.////                 /            (1/[sin(B)*cos(B)]) * tan(x)//    tan(Arg) = sgn_r * | -cot(B) + --------------------------------//                 \                     tan(B)  +  tan(x)//                                        \ //                       + CORR  |//                                     /// where////    CORR = sgn_r*c*cot(B)*SC_inv(B);  SC_inv(B) = 1/(sin(B)*cos(B)).////// The actual algorithm prescribes how all the mathematical formulas// are calculated.////// 2. Algorithmic Description// ==========================//// 2.1 Computation for Cases 2 and 4.// ----------------------------------//// For Case 2, we use two-term polynomials.////    For N even,////    rsq := r * r//    Result := c + r * rsq * P1_1//    Result := r + Result          ...in user-defined rounding////    For N odd,//    S_hi  := -frcpa(r)               ...8 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits//    S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )//    ...S_hi + S_lo is -1/(r+c) to extra precision//    S_lo  := S_lo + Q1_1*r////    Result := S_hi + S_lo     ...in user-defined rounding//// For Case 4, we use three-term polynomials////    For N even,////    rsq := r * r//    Result := c + r * rsq * (P1_1 + rsq * P1_2)//    Result := r + Result          ...in user-defined rounding////    For N odd,//    S_hi  := -frcpa(r)               ...8 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits//    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits//    S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )//    ...S_hi + S_lo is -1/(r+c) to extra precision//    rsq   := r * r//    P      := Q1_1 + rsq*Q1_2//    S_lo  := S_lo + r*P////    Result := S_hi + S_lo     ...in user-defined rounding////// Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are// the same as those used in the small_r case of Cases 1 and 3// below.////// 2.2 Computation for Cases 1 and 3.// ----------------------------------// This is further divided into the case of small_r,// where |r| < 2^(-2), and the case of normal_r, where |r| lies between// 2^(-2) and pi/4.//// Algorithm for the case of small_r// ---------------------------------//// For N even,//      rsq   := r * r//      Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3))//      r_to_the_8    := rsq * rsq//      r_to_the_8    := r_to_the_8 * r_to_the_8//      Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9))//      CORR  := c * ( 1 + rsq )//      Poly  := Poly1 + r_to_the_8*Poly2//      Result := r*Poly + CORR//      Result := r + Result     ...in user-defined rounding//      ...note that Poly1 and r_to_the_8 can be computed in parallel//      ...with Poly2 (Poly1 is intentionally set to be much//      ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden)//// For N odd,//      S_hi  := -frcpa(r)               ...8 bits//      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits//      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits//      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits//      S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )//      ...S_hi + S_lo is -1/(r+c) to extra precision//      S_lo  := S_lo + Q1_1*c////      ...S_hi and S_lo are computed in parallel with//      ...the following//      rsq := r*r//      P   := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7))////      Result :=  r*P + S_lo//      Result :=  S_hi  +  Result      ...in user-defined rounding////// Algorithm for the case of normal_r// ----------------------------------//// Here, we first consider the computation of tan( r + c ). As// presented in the previous section,////      tan( r + c )  =  tan(r) + c * sec^2(r)//                 =  sgn_r * [ tan(B+x) + CORR ]//      CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)]//// because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits.////      tan( r + c ) =//           /           (1/[sin(B)*cos(B)]) * tan(x)//      sgn_r * | tan(B) + --------------------------------  +//           \                     cot(B)  -  tan(x)//                                \ //                          CORR  |//                                ///// The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are// calculated beforehand and stored in a table. Specifically,// the table values are////      tan(B)                as  T_hi  +  T_lo;//      cot(B)             as  C_hi  +  C_lo;//      1/[sin(B)*cos(B)]  as  SC_inv//// T_hi, C_hi are in  double-precision  memory format;// T_lo, C_lo are in  single-precision  memory format;// SC_inv     is  in extended-precision memory format.//// The value of tan(x) will be approximated by a short polynomial of// the form////      tan(x)  as  x  +  x * P, where//           P  =   x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))//// Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x)// to a relative accuracy better than 2^(-20). Thus, a good// initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative// division is:////      1/(cot(B) - tan(x))      is approximately//      1/(cot(B) -   x)         is//      tan(B)/(1 - x*tan(B))    is approximately//      T_hi / ( 1 - T_hi * x )  is approximately////      T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ]//// The calculation of tan(r+c) therefore proceed as follows:////      Tx     := T_hi * x//      xsq     := x * x////      V_hi     := T_hi*(1 + Tx*(1 + Tx))//      P     := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))//      ...V_hi serves as an initial guess of 1/(cot(B) - tan(x))//         ...good to about 20 bits of accuracy////      tanx     := x + x*P//      D     := C_hi - tanx//      ...D is a double precision denominator: cot(B) - tan(x)////      V_hi     := V_hi + V_hi*(1 - V_hi*D)//      ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits////      V_lo     := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ]//                           - V_hi*C_lo )   ...observe all order//         ...V_hi + V_lo approximates 1/(cot(B) - tan(x))//      ...to extra accuracy////      ...               SC_inv(B) * (x + x*P)//      ...   tan(B) +      ------------------------- + CORR//         ...                cot(B) - (x + x*P)//      ...//      ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR//      ...////      Sx     := SC_inv * x//      CORR     := sgn_r * c * SC_inv * T_hi////      ...put the ingredients together to compute//      ...               SC_inv(B) * (x + x*P)//      ...   tan(B) +      ------------------------- + CORR//         ...                cot(B) - (x + x*P)//      ...//      ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR//      ...//      ... = T_hi + T_lo + CORR +//      ...    Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)////      CORR := CORR + T_lo//      tail := V_lo + P*(V_hi + V_lo)//         tail := Sx * tail  +  CORR//      tail := Sx * V_hi  +  tail//         T_hi := sgn_r * T_hi////         ...T_hi + sgn_r*tail  now approximate//      ...sgn_r*(tan(B+x) + CORR) accurately////      Result :=  T_hi + sgn_r*tail  ...in user-defined//                           ...rounding control//      ...It is crucial that independent paths be fully//      ...exploited for performance's sake.////// Next, we consider the computation of -cot( r + c ). As// presented in the previous section,////        -cot( r + c )  =  -cot(r) + c * csc^2(r)//                 =  sgn_r * [ -cot(B+x) + CORR ]//      CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)]//// because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits.