.file "libm_lgamma.s"// Copyright (c) 2002 - 2005, Intel Corporation// All rights reserved.//// Contributed 2002 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// 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://www.intel.com/software/products/opensource/libraries/num.htm.////*********************************************************************//// History:// 01/10/02  Initial version// 01/25/02  Corrected error tag numbers// 02/04/02  Added support of SIGN(GAMMA(x)) calculation// 05/20/02  Cleaned up namespace and sf0 syntax// 09/15/02  Fixed bug on the branch lgamma_negrecursion// 10/21/02  Now it returns SIGN(GAMMA(x))=-1 for negative zero// 02/10/03  Reordered header: .section, .global, .proc, .align// 07/22/03  Reformatted some data tables// 03/31/05  Reformatted delimiters between data tables////*********************************************************************////*********************************************************************//// Function: __libm_lgamma(double x, int* signgam, int szsigngam)// computes the principle value of the logarithm of the GAMMA function// of x. Signum of GAMMA(x) is stored to memory starting at the address// specified by the signgam.////*********************************************************************//// Resources Used:////    Floating-Point Registers: f6-f15//                              f32-f122////    General Purpose Registers://      r8-r11//      r14-r31//      r32-r36//      r37-r40 (Used to pass arguments to error handling routine)////    Predicate Registers:      p6-p15////*********************************************************************//// IEEE Special Conditions:////    __libm_lgamma(+inf) = +inf//    __libm_lgamma(-inf) = QNaN//    __libm_lgamma(+/-0) = +inf//    __libm_lgamma(x<0, x - integer) = +inf//    __libm_lgamma(SNaN) = QNaN//    __libm_lgamma(QNaN) = QNaN////*********************************************************************//// Overview//// The method consists of three cases.//// If      512 <= x < OVERFLOW_BOUNDARY   use case lgamma_pstirling;// else if 1 < x < 512                    use case lgamma_regular;// else if -17 < x < 1                    use case lgamma_negrecursion;// else if -512 <  x < -17                use case lgamma_negpoly;// else if x < -512                       use case lgamma_negstirling;// else if x is close to negative//         roots of ln(GAMMA(x))          use case lgamma_negroots;////// Case 512 <= x < OVERFLOW_BOUNDARY// ---------------------------------//   Here we use algorithm based on the Stirling formula://   ln(GAMMA(x)) = ln(sqrt(2*Pi)) + (x-0.5)ln(x) - x + (W2 + W4/x^2)/x//// Case 1 < x < 512// ----------------//   To calculate GAMMA(x) on this interval we use polynomial approximation//   on following intervals [0.875; 1.25), [1.25; 1.75), [1.75, 2.25),//   [2.25; 4), [2^i; 2^(i+1)), i=2..8////   Following variants of approximation and argument reduction are used://    1. [0.875; 1.25)//       ln(GAMMA(x)) ~ (x-1.0)*P17(x-1.0)////    2. [1.25; 1.75)//       ln(GAMMA(x)) ~ (x-LocalMinimun)*P17(x-LocalMinimun)////    3. [1.75, 2.25)//       ln(GAMMA(x)) ~ (x-2.0)*P17(x-2.0)////    4. [2.25; 4)//       ln(GAMMA(x)) ~ P22(x)////    5. [2^i; 2^(i+1)), i=2..8//       ln(GAMMA(x)) ~ P22((x-2^i)/2^i)//// Case -17 < x < 1// ----------------//   Here we use the recursive formula://   ln(GAMMA(x)) = ln(GAMMA(x+1)) - ln(x)////   Using this formula we reduce argument to base interval [1.0; 2.0]//// Case -512 <  x < -17// --------------------//   Here we use the formula://   ln(GAMMA(-x)) = ln(Pi/(x*GAMMA(x)*sin(Pi*x))) =//   = -ln(x) - ln((GAMMA(x)) - ln(sin(Pi*r)/(Pi*r)) - ln(|r|)//   where r = x - rounded_to_nearest(x), i.e |r| <= 0.5 and//   ln(sin(Pi*r)/(Pi*r)) is approximated by 14-degree polynomial of r^2////// Case x < -512// -------------//   Here we use algorithm based on the Stirling formula://   ln(GAMMA(-x)) = -ln(sqrt(2*Pi)) + (-x-0.5)ln(x) + x - (W2 + W4/x^2)/x -//   - ln(sin(Pi*r)/(Pi*r)) - ln(|r|)//   where r = x - rounded_to_nearest(x).//// Neighbourhoods of negative roots// --------------------------------//   Here we use polynomial approximation//   ln(GAMMA(x-x0)) = ln(GAMMA(x0)) + (x-x0)*P14(x-x0),//   where x0 is a root of ln(GAMMA(x)) rounded to nearest double//   precision number.////*********************************************************************FR_X                   = f10FR_Y                   = f1 // __libm_lgamma is single argument functionFR_RESULT              = f8FR_B11                 = f6FR_B10                 = f7FR_int_N               = f9FR_N                   = f10FR_P5                  = f11FR_P4                  = f12FR_P3                  = f13FR_P2                  = f14FR_NormX               = f15FR_Ln2                 = f32FR_C01                 = f33FR_A17                 = f33FR_C00                 = f34FR_Xp2                 = f34FR_A00                 = f34FR_A16                 = f34FR_C11                 = f35FR_A15                 = f35FR_C10                 = f36FR_Xp3                 = f36FR_A14                 = f36FR_B1                  = f36FR_C21                 = f37FR_A13                 = f37FR_PR01                = f37FR_C20                 = f38FR_Xp6                 = f38FR_A12                 = f38FR_C31                 = f39FR_Xp7                 = f39FR_B0                  = f39FR_A11                 = f39FR_C30                 = f40FR_Xp8                 = f40FR_A10                 = f40FR_PR00                = f40FR_C41                 = f41FR_Xp9                 = f41FR_A9                  = f41FR_PR11                = f41FR_C40                 = f42FR_A8                  = f42FR_C51                 = f43FR_Xp11                = f43FR_A7                  = f43FR_C50                 = f44FR_C                   = f44FR_Xp12                = f44FR_A6                  = f44FR_Xm2                 = f45FR_Xp13                = f45FR_A5                  = f45FR_PR10                = f45FR_C61                 = f46FR_Xp14                = f46FR_A4                  = f46FR_PR21                = f46FR_C60                 = f47FR_Xp15                = f47FR_A3                  = f47FR_PR20                = f47FR_C71                 = f48FR_Xp16                = f48FR_A2                  = f48FR_PR31                = f48FR_C70                 = f49FR_Xp17                = f49FR_A1                  = f49FR_PR30                = f49FR_C81                 = f50FR_B17                 = f50FR_A0                  = f50FR_C80                 = f51FR_B16                 = f51FR_C91                 = f52FR_B15                 = f52FR_C90                 = f53FR_B14                 = f53FR_CA1                 = f54FR_B13                 = f54FR_CA0                 = f55FR_B12                 = f55FR_CN                  = f56FR_Qlo                 = f56FR_PRN                 = f56FR_B7                  = f57FR_B6                  = f58FR_Qhi                 = f59FR_x                   = f60FR_x2                  = f61FR_TpNxLn2             = f62FR_W2                  = f63FR_x4                  = f64FR_r4                  = f64FR_x8                  = f65FR_r8                  = f65FR_r05                 = f66FR_Xm05                = f66FR_B5                  = f66FR_LnSqrt2Pi           = f67FR_B4                  = f67FR_InvX                = f68FR_B3                  = f68FR_InvX2               = f69FR_B2                  = f69FR_W4                  = f70FR_OvfBound            = f71FR_05                  = f72FR_LocalMin            = f73FR_tmp                 = f73FR_LnX                 = f74FR_Xf                  = f75FR_InvXf               = f76FR_rf                  = f77FR_rf2                 = f78FR_P54f                = f79FR_P32f                = f80FR_rf3                 = f81FR_P10f                = f82FR_TpNxLn2f            = f83FR_Nf                  = f84FR_LnXf                = f85FR_int_Nf              = f86FR_Tf                  = f87FR_Xf2                 = f88FR_Xp10                = f89FR_w3                  = f90FR_S28                 = f90FR_w2                  = f91FR_S26                 = f91FR_w6                  = f92FR_S24                 = f92FR_w4                  = f93FR_S22                 = f93FR_w                   = f94FR_S20                 = f94FR_Q8                  = f95FR_S18                 = f95FR_Q7                  = f96FR_S16                 = f96FR_Q4                  = f97FR_S14                 = f97FR_Q3                  = f98FR_S12                 = f98FR_Q6                  = f99FR_S10                 = f99FR_Q5                  = f100FR_S8                  = f100FR_Q2                  = f101FR_S6                  = f101FR_Root                = f101FR_S4                  = f102FR_Q1                  = f102FR_S2                  = f103FR_Xp1                 = f104FR_Xf4                 = f105FR_Xf8                 = f106FR_Xfr                 = f107FR_Xf6                 = f108FR_Ntrunc              = f109FR_B9                  = f110FR_2                   = f110FR_B8                  = f111FR_3                   = f111FR_5                   = f112FR_Xp4                 = f113FR_Xp5                 = f114FR_P54                 = f115FR_P32                 = f116FR_P10                 = f117FR_r                   = f118FR_r2                  = f119FR_r3                  = f120FR_T                   = f121FR_int_Ntrunc          = f122//===================================GR_TAG                 = r8GR_ExpMask             = r8GR_ExpBias             = r9GR_ad_Roots            = r9GR_Expf                = r10