📄 libm_lgamma.s
字号:
.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⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -