1.file "libm_lgammaf.s" 2 3 4// Copyright (c) 2002 - 2005, Intel Corporation 5// All rights reserved. 6// 7// 8// Redistribution and use in source and binary forms, with or without 9// modification, are permitted provided that the following conditions are 10// met: 11// 12// * Redistributions of source code must retain the above copyright 13// notice, this list of conditions and the following disclaimer. 14// 15// * Redistributions in binary form must reproduce the above copyright 16// notice, this list of conditions and the following disclaimer in the 17// documentation and/or other materials provided with the distribution. 18// 19// * The name of Intel Corporation may not be used to endorse or promote 20// products derived from this software without specific prior written 21// permission. 22 23// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 24// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,INCLUDING,BUT NOT 25// LIMITED TO,THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 26// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS 27// CONTRIBUTORS BE LIABLE FOR ANY DIRECT,INDIRECT,INCIDENTAL,SPECIAL, 28// EXEMPLARY,OR CONSEQUENTIAL DAMAGES (INCLUDING,BUT NOT LIMITED TO, 29// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,DATA,OR 30// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 31// OF LIABILITY,WHETHER IN CONTRACT,STRICT LIABILITY OR TORT (INCLUDING 32// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 33// SOFTWARE,EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 34// 35// Intel Corporation is the author of this code,and requests that all 36// problem reports or change requests be submitted to it directly at 37// http://www.intel.com/software/products/opensource/libraries/num.htm. 38// 39//********************************************************************* 40// 41// History: 42// 01/10/02 Initial version 43// 01/25/02 Corrected parameter store, load, and tag for __libm_error_support 44// 02/01/02 Added support of SIGN(GAMMA(x)) calculation 45// 05/20/02 Cleaned up namespace and sf0 syntax 46// 09/16/02 Improved accuracy on intervals reduced to [1;1.25] 47// 10/21/02 Now it returns SIGN(GAMMA(x))=-1 for negative zero 48// 02/10/03 Reordered header: .section, .global, .proc, .align 49// 07/22/03 Reformatted some data tables 50// 03/31/05 Reformatted delimiters between data tables 51// 52//********************************************************************* 53// 54//********************************************************************* 55// 56// Function: __libm_lgammaf(float x, int* signgam, int szsigngam) 57// computes the principle value of the logarithm of the GAMMA function 58// of x. Signum of GAMMA(x) is stored to memory starting at the address 59// specified by the signgam. 60// 61//********************************************************************* 62// 63// Resources Used: 64// 65// Floating-Point Registers: f6-f15 66// f32-f97 67// 68// General Purpose Registers: 69// r8-r11 70// r14-r30 71// r32-r36 72// r37-r40 (Used to pass arguments to error handling routine) 73// 74// Predicate Registers: p6-p15 75// 76//********************************************************************* 77// 78// IEEE Special Conditions: 79// 80// lgamma(+inf) = +inf 81// lgamma(-inf) = +inf 82// lgamma(+/-0) = +inf 83// lgamma(x<0, x - integer) = +inf 84// lgamma(SNaN) = QNaN 85// lgamma(QNaN) = QNaN 86// 87//********************************************************************* 88// 89// Overview 90// 91// The method consists of three cases. 92// 93// If 2^13 <= x < OVERFLOW_BOUNDARY use case lgammaf_pstirling; 94// else if 1 < x < 2^13 use case lgammaf_regular; 95// else if -9 < x < 1 use case lgammaf_negrecursion; 96// else if -2^13 < x < -9 use case lgammaf_negpoly; 97// else if x < -2^13 use case lgammaf_negstirling; 98// else if x is close to negative 99// roots of ln(GAMMA(x)) use case lgammaf_negroots; 100// 101// 102// Case 2^13 <= x < OVERFLOW_BOUNDARY 103// ---------------------------------- 104// Here we use algorithm based on the Stirling formula: 105// ln(GAMMA(x)) = ln(sqrt(2*Pi)) + (x-0.5)*ln(x) - x 106// 107// Case 1 < x < 2^13 108// ----------------- 109// To calculate ln(GAMMA(x)) for such arguments we use polynomial 110// approximation on following intervals: [1.0; 1.25), [1.25; 1.5), 111// [1.5, 1.75), [1.75; 2), [2; 4), [2^i; 2^(i+1)), i=1..8 112// 113// Following variants of approximation and argument reduction are used: 114// 1. [1.0; 1.25) 115// ln(GAMMA(x)) ~ (x-1.0)*P7(x) 116// 117// 2. [1.25; 1.5) 118// ln(GAMMA(x)) ~ ln(GAMMA(x0))+(x-x0)*P8(x-x0), 119// where x0 - point of local minimum on [1;2] rounded to nearest double 120// precision number. 121// 122// 3. [1.5; 1.75) 123// ln(GAMMA(x)) ~ P8(x) 124// 125// 4. [1.75; 2.0) 126// ln(GAMMA(x)) ~ (x-2)*P7(x) 127// 128// 5. [2; 4) 129// ln(GAMMA(x)) ~ (x-2)*P10(x) 130// 131// 6. [2^i; 2^(i+1)), i=2..8 132// ln(GAMMA(x)) ~ P10((x-2^i)/2^i) 133// 134// Case -9 < x < 1 135// --------------- 136// Here we use the recursive formula: 137// ln(GAMMA(x)) = ln(GAMMA(x+1)) - ln(x) 138// 139// Using this formula we reduce argument to base interval [1.0; 2.0] 140// 141// Case -2^13 < x < -9 142// -------------------- 143// Here we use the formula: 144// ln(GAMMA(x)) = ln(Pi/(|x|*GAMMA(|x|)*sin(Pi*|x|))) = 145// = -ln(|x|) - ln((GAMMA(|x|)) - ln(sin(Pi*r)/(Pi*r)) - ln(|r|) 146// where r = x - rounded_to_nearest(x), i.e |r| <= 0.5 and 147// ln(sin(Pi*r)/(Pi*r)) is approximated by 8-degree polynomial of r^2 148// 149// Case x < -2^13 150// -------------- 151// Here we use algorithm based on the Stirling formula: 152// ln(GAMMA(x)) = -ln(sqrt(2*Pi)) + (|x|-0.5)ln(x) - |x| - 153// - ln(sin(Pi*r)/(Pi*r)) - ln(|r|) 154// where r = x - rounded_to_nearest(x). 155// 156// Neighbourhoods of negative roots 157// -------------------------------- 158// Here we use polynomial approximation 159// ln(GAMMA(x-x0)) = ln(GAMMA(x0)) + (x-x0)*P14(x-x0), 160// where x0 is a root of ln(GAMMA(x)) rounded to nearest double 161// precision number. 162// 163// 164// Claculation of logarithm 165// ------------------------ 166// Consider x = 2^N * xf so 167// ln(x) = ln(frcpa(x)*x/frcpa(x)) 168// = ln(1/frcpa(x)) + ln(frcpa(x)*x) 169// 170// frcpa(x) = 2^(-N) * frcpa(xf) 171// 172// ln(1/frcpa(x)) = -ln(2^(-N)) - ln(frcpa(xf)) 173// = N*ln(2) - ln(frcpa(xf)) 174// = N*ln(2) + ln(1/frcpa(xf)) 175// 176// ln(x) = ln(1/frcpa(x)) + ln(frcpa(x)*x) = 177// = N*ln(2) + ln(1/frcpa(xf)) + ln(frcpa(x)*x) 178// = N*ln(2) + T + ln(frcpa(x)*x) 179// 180// Let r = 1 - frcpa(x)*x, note that r is quite small by 181// absolute value so 182// 183// ln(x) = N*ln(2) + T + ln(1+r) ~ N*ln(2) + T + Series(r), 184// where T - is precomputed tabular value, 185// Series(r) = (P3*r + P2)*r^2 + (P1*r + 1) 186// 187//********************************************************************* 188 189GR_TAG = r8 190GR_ad_Data = r8 191GR_ad_Co = r9 192GR_ad_SignGam = r10 193GR_ad_Ce = r10 194GR_SignExp = r11 195 196GR_ad_C650 = r14 197GR_ad_RootCo = r14 198GR_ad_C0 = r15 199GR_Dx = r15 200GR_Ind = r16 201GR_Offs = r17 202GR_IntNum = r17 203GR_ExpBias = r18 204GR_ExpMask = r19 205GR_Ind4T = r20 206GR_RootInd = r20 207GR_Sig = r21 208GR_Exp = r22 209GR_PureExp = r23 210GR_ad_C43 = r24 211GR_StirlBound = r25 212GR_ad_T = r25 213GR_IndX8 = r25 214GR_Neg2 = r25 215GR_2xDx = r25 216GR_SingBound = r26 217GR_IndX2 = r26 218GR_Neg4 = r26 219GR_ad_RootCe = r26 220GR_Arg = r27 221GR_ExpOf2 = r28 222GR_fff7 = r28 223GR_Root = r28 224GR_ReqBound = r28 225GR_N = r29 226GR_ad_Root = r30 227GR_ad_OvfBound = r30 228GR_SignOfGamma = r31 229 230GR_SAVE_B0 = r33 231GR_SAVE_PFS = r34 232GR_SAVE_GP = r35 233GR_SAVE_SP = r36 234 235GR_Parameter_X = r37 236GR_Parameter_Y = r38 237GR_Parameter_RESULT = r39 238GR_Parameter_TAG = r40 239 240//********************************************************************* 241 242FR_X = f10 243FR_Y = f1 // lgammaf is single argument function 244FR_RESULT = f8 245 246FR_x = f6 247FR_x2 = f7 248 249FR_x3 = f9 250FR_x4 = f10 251FR_xm2 = f11 252FR_w = f11 253FR_w2 = f12 254FR_Q32 = f13 255FR_Q10 = f14 256FR_InvX = f15 257 258FR_NormX = f32 259 260FR_A0 = f33 261FR_A1 = f34 262FR_A2 = f35 263FR_A3 = f36 264FR_A4 = f37 265FR_A5 = f38 266FR_A6 = f39 267FR_A7 = f40 268FR_A8 = f41 269FR_A9 = f42 270FR_A10 = f43 271 272FR_int_N = f44 273FR_P3 = f45 274FR_P2 = f46 275FR_P1 = f47 276FR_LocalMin = f48 277FR_Ln2 = f49 278FR_05 = f50 279FR_LnSqrt2Pi = f51 280FR_3 = f52 281FR_r = f53 282FR_r2 = f54 283FR_T = f55 284FR_N = f56 285FR_xm05 = f57 286FR_int_Ln = f58 287FR_P32 = f59 288FR_P10 = f60 289 290FR_Xf = f61 291FR_InvXf = f62 292FR_rf = f63 293FR_rf2 = f64 294FR_Tf = f65 295FR_Nf = f66 296FR_xm05f = f67 297FR_P32f = f68 298FR_P10f = f69 299FR_Lnf = f70 300FR_Xf2 = f71 301FR_Xf4 = f72 302FR_Xf8 = f73 303FR_Ln = f74 304FR_xx = f75 305FR_Root = f75 306FR_Req = f76 307FR_1pXf = f77 308 309FR_S16 = f78 310FR_R3 = f78 311FR_S14 = f79 312FR_R2 = f79 313FR_S12 = f80 314FR_R1 = f80 315FR_S10 = f81 316FR_R0 = f81 317FR_S8 = f82 318FR_rx = f82 319FR_S6 = f83 320FR_rx2 = f84 321FR_S4 = f84 322FR_S2 = f85 323 324FR_Xp1 = f86 325FR_Xp2 = f87 326FR_Xp3 = f88 327FR_Xp4 = f89 328FR_Xp5 = f90 329FR_Xp6 = f91 330FR_Xp7 = f92 331FR_Xp8 = f93 332FR_OverflowBound = f93 333 334FR_2 = f94 335FR_tmp = f95 336FR_int_Ntrunc = f96 337FR_Ntrunc = f97 338 339//********************************************************************* 340 341RODATA 342.align 32 343LOCAL_OBJECT_START(lgammaf_data) 344log_table_1: 345data8 0xbfd0001008f39d59 // P3 346data8 0x3fd5556073e0c45a // P2 347data8 0x3fe62e42fefa39ef // ln(2) 348data8 0x3fe0000000000000 // 0.5 349// 350data8 0x3F60040155D5889E //ln(1/frcpa(1+ 0/256) 351data8 0x3F78121214586B54 //ln(1/frcpa(1+ 1/256) 352data8 0x3F841929F96832F0 //ln(1/frcpa(1+ 2/256) 353data8 0x3F8C317384C75F06 //ln(1/frcpa(1+ 3/256) 354data8 0x3F91A6B91AC73386 //ln(1/frcpa(1+ 4/256) 355data8 0x3F95BA9A5D9AC039 //ln(1/frcpa(1+ 5/256) 356data8 0x3F99D2A8074325F4 //ln(1/frcpa(1+ 6/256) 357data8 0x3F9D6B2725979802 //ln(1/frcpa(1+ 7/256) 358data8 0x3FA0C58FA19DFAAA //ln(1/frcpa(1+ 8/256) 359data8 0x3FA2954C78CBCE1B //ln(1/frcpa(1+ 9/256) 360data8 0x3FA4A94D2DA96C56 //ln(1/frcpa(1+ 10/256) 361data8 0x3FA67C94F2D4BB58 //ln(1/frcpa(1+ 11/256) 362data8 0x3FA85188B630F068 //ln(1/frcpa(1+ 12/256) 363data8 0x3FAA6B8ABE73AF4C //ln(1/frcpa(1+ 13/256) 364data8 0x3FAC441E06F72A9E //ln(1/frcpa(1+ 14/256) 365data8 0x3FAE1E6713606D07 //ln(1/frcpa(1+ 15/256) 366data8 0x3FAFFA6911AB9301 //ln(1/frcpa(1+ 16/256) 367data8 0x3FB0EC139C5DA601 //ln(1/frcpa(1+ 17/256) 368data8 0x3FB1DBD2643D190B //ln(1/frcpa(1+ 18/256) 369data8 0x3FB2CC7284FE5F1C //ln(1/frcpa(1+ 19/256) 370data8 0x3FB3BDF5A7D1EE64 //ln(1/frcpa(1+ 20/256) 371data8 0x3FB4B05D7AA012E0 //ln(1/frcpa(1+ 21/256) 372data8 0x3FB580DB7CEB5702 //ln(1/frcpa(1+ 22/256) 373data8 0x3FB674F089365A7A //ln(1/frcpa(1+ 23/256) 374data8 0x3FB769EF2C6B568D //ln(1/frcpa(1+ 24/256) 375data8 0x3FB85FD927506A48 //ln(1/frcpa(1+ 25/256) 376data8 0x3FB9335E5D594989 //ln(1/frcpa(1+ 26/256) 377data8 0x3FBA2B0220C8E5F5 //ln(1/frcpa(1+ 27/256) 378data8 0x3FBB0004AC1A86AC //ln(1/frcpa(1+ 28/256) 379data8 0x3FBBF968769FCA11 //ln(1/frcpa(1+ 29/256) 380data8 0x3FBCCFEDBFEE13A8 //ln(1/frcpa(1+ 30/256) 381data8 0x3FBDA727638446A2 //ln(1/frcpa(1+ 31/256) 382data8 0x3FBEA3257FE10F7A //ln(1/frcpa(1+ 32/256) 383data8 0x3FBF7BE9FEDBFDE6 //ln(1/frcpa(1+ 33/256) 384data8 0x3FC02AB352FF25F4 //ln(1/frcpa(1+ 34/256) 385data8 0x3FC097CE579D204D //ln(1/frcpa(1+ 35/256) 386data8 0x3FC1178E8227E47C //ln(1/frcpa(1+ 36/256) 387data8 0x3FC185747DBECF34 //ln(1/frcpa(1+ 37/256) 388data8 0x3FC1F3B925F25D41 //ln(1/frcpa(1+ 38/256) 389data8 0x3FC2625D1E6DDF57 //ln(1/frcpa(1+ 39/256) 390data8 0x3FC2D1610C86813A //ln(1/frcpa(1+ 40/256) 391data8 0x3FC340C59741142E //ln(1/frcpa(1+ 41/256) 392data8 0x3FC3B08B6757F2A9 //ln(1/frcpa(1+ 42/256) 393data8 0x3FC40DFB08378003 //ln(1/frcpa(1+ 43/256) 394data8 0x3FC47E74E8CA5F7C //ln(1/frcpa(1+ 44/256) 395data8 0x3FC4EF51F6466DE4 //ln(1/frcpa(1+ 45/256) 396data8 0x3FC56092E02BA516 //ln(1/frcpa(1+ 46/256) 397data8 0x3FC5D23857CD74D5 //ln(1/frcpa(1+ 47/256) 398data8 0x3FC6313A37335D76 //ln(1/frcpa(1+ 48/256) 399data8 0x3FC6A399DABBD383 //ln(1/frcpa(1+ 49/256) 400data8 0x3FC70337DD3CE41B //ln(1/frcpa(1+ 50/256) 401data8 0x3FC77654128F6127 //ln(1/frcpa(1+ 51/256) 402data8 0x3FC7E9D82A0B022D //ln(1/frcpa(1+ 52/256) 403data8 0x3FC84A6B759F512F //ln(1/frcpa(1+ 53/256) 404data8 0x3FC8AB47D5F5A310 //ln(1/frcpa(1+ 54/256) 405data8 0x3FC91FE49096581B //ln(1/frcpa(1+ 55/256) 406data8 0x3FC981634011AA75 //ln(1/frcpa(1+ 56/256) 407data8 0x3FC9F6C407089664 //ln(1/frcpa(1+ 57/256) 408data8 0x3FCA58E729348F43 //ln(1/frcpa(1+ 58/256) 409data8 0x3FCABB55C31693AD //ln(1/frcpa(1+ 59/256) 410data8 0x3FCB1E104919EFD0 //ln(1/frcpa(1+ 60/256) 411data8 0x3FCB94EE93E367CB //ln(1/frcpa(1+ 61/256) 412data8 0x3FCBF851C067555F //ln(1/frcpa(1+ 62/256) 413data8 0x3FCC5C0254BF23A6 //ln(1/frcpa(1+ 63/256) 414data8 0x3FCCC000C9DB3C52 //ln(1/frcpa(1+ 64/256) 415data8 0x3FCD244D99C85674 //ln(1/frcpa(1+ 65/256) 416data8 0x3FCD88E93FB2F450 //ln(1/frcpa(1+ 66/256) 417data8 0x3FCDEDD437EAEF01 //ln(1/frcpa(1+ 67/256) 418data8 0x3FCE530EFFE71012 //ln(1/frcpa(1+ 68/256) 419data8 0x3FCEB89A1648B971 //ln(1/frcpa(1+ 69/256) 420data8 0x3FCF1E75FADF9BDE //ln(1/frcpa(1+ 70/256) 421data8 0x3FCF84A32EAD7C35 //ln(1/frcpa(1+ 71/256) 422data8 0x3FCFEB2233EA07CD //ln(1/frcpa(1+ 72/256) 423data8 0x3FD028F9C7035C1C //ln(1/frcpa(1+ 73/256) 424data8 0x3FD05C8BE0D9635A //ln(1/frcpa(1+ 74/256) 425data8 0x3FD085EB8F8AE797 //ln(1/frcpa(1+ 75/256) 426data8 0x3FD0B9C8E32D1911 //ln(1/frcpa(1+ 76/256) 427data8 0x3FD0EDD060B78081 //ln(1/frcpa(1+ 77/256) 428data8 0x3FD122024CF0063F //ln(1/frcpa(1+ 78/256) 429data8 0x3FD14BE2927AECD4 //ln(1/frcpa(1+ 79/256) 430data8 0x3FD180618EF18ADF //ln(1/frcpa(1+ 80/256) 431data8 0x3FD1B50BBE2FC63B //ln(1/frcpa(1+ 81/256) 432data8 0x3FD1DF4CC7CF242D //ln(1/frcpa(1+ 82/256) 433data8 0x3FD214456D0EB8D4 //ln(1/frcpa(1+ 83/256) 434data8 0x3FD23EC5991EBA49 //ln(1/frcpa(1+ 84/256) 435data8 0x3FD2740D9F870AFB //ln(1/frcpa(1+ 85/256) 436data8 0x3FD29ECDABCDFA04 //ln(1/frcpa(1+ 86/256) 437data8 0x3FD2D46602ADCCEE //ln(1/frcpa(1+ 87/256) 438data8 0x3FD2FF66B04EA9D4 //ln(1/frcpa(1+ 88/256) 439data8 0x3FD335504B355A37 //ln(1/frcpa(1+ 89/256) 440data8 0x3FD360925EC44F5D //ln(1/frcpa(1+ 90/256) 441data8 0x3FD38BF1C3337E75 //ln(1/frcpa(1+ 91/256) 442data8 0x3FD3C25277333184 //ln(1/frcpa(1+ 92/256) 443data8 0x3FD3EDF463C1683E //ln(1/frcpa(1+ 93/256) 444data8 0x3FD419B423D5E8C7 //ln(1/frcpa(1+ 94/256) 445data8 0x3FD44591E0539F49 //ln(1/frcpa(1+ 95/256) 446data8 0x3FD47C9175B6F0AD //ln(1/frcpa(1+ 96/256) 447data8 0x3FD4A8B341552B09 //ln(1/frcpa(1+ 97/256) 448data8 0x3FD4D4F3908901A0 //ln(1/frcpa(1+ 98/256) 449data8 0x3FD501528DA1F968 //ln(1/frcpa(1+ 99/256) 450data8 0x3FD52DD06347D4F6 //ln(1/frcpa(1+ 100/256) 451data8 0x3FD55A6D3C7B8A8A //ln(1/frcpa(1+ 101/256) 452data8 0x3FD5925D2B112A59 //ln(1/frcpa(1+ 102/256) 453data8 0x3FD5BF406B543DB2 //ln(1/frcpa(1+ 103/256) 454data8 0x3FD5EC433D5C35AE //ln(1/frcpa(1+ 104/256) 455data8 0x3FD61965CDB02C1F //ln(1/frcpa(1+ 105/256) 456data8 0x3FD646A84935B2A2 //ln(1/frcpa(1+ 106/256) 457data8 0x3FD6740ADD31DE94 //ln(1/frcpa(1+ 107/256) 458data8 0x3FD6A18DB74A58C5 //ln(1/frcpa(1+ 108/256) 459data8 0x3FD6CF31058670EC //ln(1/frcpa(1+ 109/256) 460data8 0x3FD6F180E852F0BA //ln(1/frcpa(1+ 110/256) 461data8 0x3FD71F5D71B894F0 //ln(1/frcpa(1+ 111/256) 462data8 0x3FD74D5AEFD66D5C //ln(1/frcpa(1+ 112/256) 463data8 0x3FD77B79922BD37E //ln(1/frcpa(1+ 113/256) 464data8 0x3FD7A9B9889F19E2 //ln(1/frcpa(1+ 114/256) 465data8 0x3FD7D81B037EB6A6 //ln(1/frcpa(1+ 115/256) 466data8 0x3FD8069E33827231 //ln(1/frcpa(1+ 116/256) 467data8 0x3FD82996D3EF8BCB //ln(1/frcpa(1+ 117/256) 468data8 0x3FD85855776DCBFB //ln(1/frcpa(1+ 118/256) 469data8 0x3FD8873658327CCF //ln(1/frcpa(1+ 119/256) 470data8 0x3FD8AA75973AB8CF //ln(1/frcpa(1+ 120/256) 471data8 0x3FD8D992DC8824E5 //ln(1/frcpa(1+ 121/256) 472data8 0x3FD908D2EA7D9512 //ln(1/frcpa(1+ 122/256) 473data8 0x3FD92C59E79C0E56 //ln(1/frcpa(1+ 123/256) 474data8 0x3FD95BD750EE3ED3 //ln(1/frcpa(1+ 124/256) 475data8 0x3FD98B7811A3EE5B //ln(1/frcpa(1+ 125/256) 476data8 0x3FD9AF47F33D406C //ln(1/frcpa(1+ 126/256) 477data8 0x3FD9DF270C1914A8 //ln(1/frcpa(1+ 127/256) 478data8 0x3FDA0325ED14FDA4 //ln(1/frcpa(1+ 128/256) 479data8 0x3FDA33440224FA79 //ln(1/frcpa(1+ 129/256) 480data8 0x3FDA57725E80C383 //ln(1/frcpa(1+ 130/256) 481data8 0x3FDA87D0165DD199 //ln(1/frcpa(1+ 131/256) 482data8 0x3FDAAC2E6C03F896 //ln(1/frcpa(1+ 132/256) 483data8 0x3FDADCCC6FDF6A81 //ln(1/frcpa(1+ 133/256) 484data8 0x3FDB015B3EB1E790 //ln(1/frcpa(1+ 134/256) 485data8 0x3FDB323A3A635948 //ln(1/frcpa(1+ 135/256) 486data8 0x3FDB56FA04462909 //ln(1/frcpa(1+ 136/256) 487data8 0x3FDB881AA659BC93 //ln(1/frcpa(1+ 137/256) 488data8 0x3FDBAD0BEF3DB165 //ln(1/frcpa(1+ 138/256) 489data8 0x3FDBD21297781C2F //ln(1/frcpa(1+ 139/256) 490data8 0x3FDC039236F08819 //ln(1/frcpa(1+ 140/256) 491data8 0x3FDC28CB1E4D32FD //ln(1/frcpa(1+ 141/256) 492data8 0x3FDC4E19B84723C2 //ln(1/frcpa(1+ 142/256) 493data8 0x3FDC7FF9C74554C9 //ln(1/frcpa(1+ 143/256) 494data8 0x3FDCA57B64E9DB05 //ln(1/frcpa(1+ 144/256) 495data8 0x3FDCCB130A5CEBB0 //ln(1/frcpa(1+ 145/256) 496data8 0x3FDCF0C0D18F326F //ln(1/frcpa(1+ 146/256) 497data8 0x3FDD232075B5A201 //ln(1/frcpa(1+ 147/256) 498data8 0x3FDD490246DEFA6B //ln(1/frcpa(1+ 148/256) 499data8 0x3FDD6EFA918D25CD //ln(1/frcpa(1+ 149/256) 500data8 0x3FDD9509707AE52F //ln(1/frcpa(1+ 150/256) 501data8 0x3FDDBB2EFE92C554 //ln(1/frcpa(1+ 151/256) 502data8 0x3FDDEE2F3445E4AF //ln(1/frcpa(1+ 152/256) 503data8 0x3FDE148A1A2726CE //ln(1/frcpa(1+ 153/256) 504data8 0x3FDE3AFC0A49FF40 //ln(1/frcpa(1+ 154/256) 505data8 0x3FDE6185206D516E //ln(1/frcpa(1+ 155/256) 506data8 0x3FDE882578823D52 //ln(1/frcpa(1+ 156/256) 507data8 0x3FDEAEDD2EAC990C //ln(1/frcpa(1+ 157/256) 508data8 0x3FDED5AC5F436BE3 //ln(1/frcpa(1+ 158/256) 509data8 0x3FDEFC9326D16AB9 //ln(1/frcpa(1+ 159/256) 510data8 0x3FDF2391A2157600 //ln(1/frcpa(1+ 160/256) 511data8 0x3FDF4AA7EE03192D //ln(1/frcpa(1+ 161/256) 512data8 0x3FDF71D627C30BB0 //ln(1/frcpa(1+ 162/256) 513data8 0x3FDF991C6CB3B379 //ln(1/frcpa(1+ 163/256) 514data8 0x3FDFC07ADA69A910 //ln(1/frcpa(1+ 164/256) 515data8 0x3FDFE7F18EB03D3E //ln(1/frcpa(1+ 165/256) 516data8 0x3FE007C053C5002E //ln(1/frcpa(1+ 166/256) 517data8 0x3FE01B942198A5A1 //ln(1/frcpa(1+ 167/256) 518data8 0x3FE02F74400C64EB //ln(1/frcpa(1+ 168/256) 519data8 0x3FE04360BE7603AD //ln(1/frcpa(1+ 169/256) 520data8 0x3FE05759AC47FE34 //ln(1/frcpa(1+ 170/256) 521data8 0x3FE06B5F1911CF52 //ln(1/frcpa(1+ 171/256) 522data8 0x3FE078BF0533C568 //ln(1/frcpa(1+ 172/256) 523data8 0x3FE08CD9687E7B0E //ln(1/frcpa(1+ 173/256) 524data8 0x3FE0A10074CF9019 //ln(1/frcpa(1+ 174/256) 525data8 0x3FE0B5343A234477 //ln(1/frcpa(1+ 175/256) 526data8 0x3FE0C974C89431CE //ln(1/frcpa(1+ 176/256) 527data8 0x3FE0DDC2305B9886 //ln(1/frcpa(1+ 177/256) 528data8 0x3FE0EB524BAFC918 //ln(1/frcpa(1+ 178/256) 529data8 0x3FE0FFB54213A476 //ln(1/frcpa(1+ 179/256) 530data8 0x3FE114253DA97D9F //ln(1/frcpa(1+ 180/256) 531data8 0x3FE128A24F1D9AFF //ln(1/frcpa(1+ 181/256) 532data8 0x3FE1365252BF0865 //ln(1/frcpa(1+ 182/256) 533data8 0x3FE14AE558B4A92D //ln(1/frcpa(1+ 183/256) 534data8 0x3FE15F85A19C765B //ln(1/frcpa(1+ 184/256) 535data8 0x3FE16D4D38C119FA //ln(1/frcpa(1+ 185/256) 536data8 0x3FE18203C20DD133 //ln(1/frcpa(1+ 186/256) 537data8 0x3FE196C7BC4B1F3B //ln(1/frcpa(1+ 187/256) 538data8 0x3FE1A4A738B7A33C //ln(1/frcpa(1+ 188/256) 539data8 0x3FE1B981C0C9653D //ln(1/frcpa(1+ 189/256) 540data8 0x3FE1CE69E8BB106B //ln(1/frcpa(1+ 190/256) 541data8 0x3FE1DC619DE06944 //ln(1/frcpa(1+ 191/256) 542data8 0x3FE1F160A2AD0DA4 //ln(1/frcpa(1+ 192/256) 543data8 0x3FE2066D7740737E //ln(1/frcpa(1+ 193/256) 544data8 0x3FE2147DBA47A394 //ln(1/frcpa(1+ 194/256) 545data8 0x3FE229A1BC5EBAC3 //ln(1/frcpa(1+ 195/256) 546data8 0x3FE237C1841A502E //ln(1/frcpa(1+ 196/256) 547data8 0x3FE24CFCE6F80D9A //ln(1/frcpa(1+ 197/256) 548data8 0x3FE25B2C55CD5762 //ln(1/frcpa(1+ 198/256) 549data8 0x3FE2707F4D5F7C41 //ln(1/frcpa(1+ 199/256) 550data8 0x3FE285E0842CA384 //ln(1/frcpa(1+ 200/256) 551data8 0x3FE294294708B773 //ln(1/frcpa(1+ 201/256) 552data8 0x3FE2A9A2670AFF0C //ln(1/frcpa(1+ 202/256) 553data8 0x3FE2B7FB2C8D1CC1 //ln(1/frcpa(1+ 203/256) 554data8 0x3FE2C65A6395F5F5 //ln(1/frcpa(1+ 204/256) 555data8 0x3FE2DBF557B0DF43 //ln(1/frcpa(1+ 205/256) 556data8 0x3FE2EA64C3F97655 //ln(1/frcpa(1+ 206/256) 557data8 0x3FE3001823684D73 //ln(1/frcpa(1+ 207/256) 558data8 0x3FE30E97E9A8B5CD //ln(1/frcpa(1+ 208/256) 559data8 0x3FE32463EBDD34EA //ln(1/frcpa(1+ 209/256) 560data8 0x3FE332F4314AD796 //ln(1/frcpa(1+ 210/256) 561data8 0x3FE348D90E7464D0 //ln(1/frcpa(1+ 211/256) 562data8 0x3FE35779F8C43D6E //ln(1/frcpa(1+ 212/256) 563data8 0x3FE36621961A6A99 //ln(1/frcpa(1+ 213/256) 564data8 0x3FE37C299F3C366A //ln(1/frcpa(1+ 214/256) 565data8 0x3FE38AE2171976E7 //ln(1/frcpa(1+ 215/256) 566data8 0x3FE399A157A603E7 //ln(1/frcpa(1+ 216/256) 567data8 0x3FE3AFCCFE77B9D1 //ln(1/frcpa(1+ 217/256) 568data8 0x3FE3BE9D503533B5 //ln(1/frcpa(1+ 218/256) 569data8 0x3FE3CD7480B4A8A3 //ln(1/frcpa(1+ 219/256) 570data8 0x3FE3E3C43918F76C //ln(1/frcpa(1+ 220/256) 571data8 0x3FE3F2ACB27ED6C7 //ln(1/frcpa(1+ 221/256) 572data8 0x3FE4019C2125CA93 //ln(1/frcpa(1+ 222/256) 573data8 0x3FE4181061389722 //ln(1/frcpa(1+ 223/256) 574data8 0x3FE42711518DF545 //ln(1/frcpa(1+ 224/256) 575data8 0x3FE436194E12B6BF //ln(1/frcpa(1+ 225/256) 576data8 0x3FE445285D68EA69 //ln(1/frcpa(1+ 226/256) 577data8 0x3FE45BCC464C893A //ln(1/frcpa(1+ 227/256) 578data8 0x3FE46AED21F117FC //ln(1/frcpa(1+ 228/256) 579data8 0x3FE47A1527E8A2D3 //ln(1/frcpa(1+ 229/256) 580data8 0x3FE489445EFFFCCC //ln(1/frcpa(1+ 230/256) 581data8 0x3FE4A018BCB69835 //ln(1/frcpa(1+ 231/256) 582data8 0x3FE4AF5A0C9D65D7 //ln(1/frcpa(1+ 232/256) 583data8 0x3FE4BEA2A5BDBE87 //ln(1/frcpa(1+ 233/256) 584data8 0x3FE4CDF28F10AC46 //ln(1/frcpa(1+ 234/256) 585data8 0x3FE4DD49CF994058 //ln(1/frcpa(1+ 235/256) 586data8 0x3FE4ECA86E64A684 //ln(1/frcpa(1+ 236/256) 587data8 0x3FE503C43CD8EB68 //ln(1/frcpa(1+ 237/256) 588data8 0x3FE513356667FC57 //ln(1/frcpa(1+ 238/256) 589data8 0x3FE522AE0738A3D8 //ln(1/frcpa(1+ 239/256) 590data8 0x3FE5322E26867857 //ln(1/frcpa(1+ 240/256) 591data8 0x3FE541B5CB979809 //ln(1/frcpa(1+ 241/256) 592data8 0x3FE55144FDBCBD62 //ln(1/frcpa(1+ 242/256) 593data8 0x3FE560DBC45153C7 //ln(1/frcpa(1+ 243/256) 594data8 0x3FE5707A26BB8C66 //ln(1/frcpa(1+ 244/256) 595data8 0x3FE587F60ED5B900 //ln(1/frcpa(1+ 245/256) 596data8 0x3FE597A7977C8F31 //ln(1/frcpa(1+ 246/256) 597data8 0x3FE5A760D634BB8B //ln(1/frcpa(1+ 247/256) 598data8 0x3FE5B721D295F10F //ln(1/frcpa(1+ 248/256) 599data8 0x3FE5C6EA94431EF9 //ln(1/frcpa(1+ 249/256) 600data8 0x3FE5D6BB22EA86F6 //ln(1/frcpa(1+ 250/256) 601data8 0x3FE5E6938645D390 //ln(1/frcpa(1+ 251/256) 602data8 0x3FE5F673C61A2ED2 //ln(1/frcpa(1+ 252/256) 603data8 0x3FE6065BEA385926 //ln(1/frcpa(1+ 253/256) 604data8 0x3FE6164BFA7CC06B //ln(1/frcpa(1+ 254/256) 605data8 0x3FE62643FECF9743 //ln(1/frcpa(1+ 255/256) 606// 607// [2;4) 608data8 0xBEB2CC7A38B9355F,0x3F035F2D1833BF4C // A10,A9 609data8 0xBFF51BAA7FD27785,0x3FFC9D5D5B6CDEFF // A2,A1 610data8 0xBF421676F9CB46C7,0x3F7437F2FA1436C6 // A8,A7 611data8 0xBFD7A7041DE592FE,0x3FE9F107FEE8BD29 // A4,A3 612// [4;8) 613data8 0x3F6BBBD68451C0CD,0xBF966EC3272A16F7 // A10,A9 614data8 0x40022A24A39AD769,0x4014190EDF49C8C5 // A2,A1 615data8 0x3FB130FD016EE241,0xBFC151B46E635248 // A8,A7 616data8 0x3FDE8F611965B5FE,0xBFEB5110EB265E3D // A4,A3 617// [8;16) 618data8 0x3F736EF93508626A,0xBF9FE5DBADF58AF1 // A10,A9 619data8 0x40110A9FC5192058,0x40302008A6F96B29 // A2,A1 620data8 0x3FB8E74E0CE1E4B5,0xBFC9B5DA78873656 // A8,A7 621data8 0x3FE99D0DF10022DC,0xBFF829C0388F9484 // A4,A3 622// [16;32) 623data8 0x3F7FFF9D6D7E9269,0xBFAA780A249AEDB1 // A10,A9 624data8 0x402082A807AEA080,0x4045ED9868408013 // A2,A1 625data8 0x3FC4E1E54C2F99B7,0xBFD5DE2D6FFF1490 // A8,A7 626data8 0x3FF75FC89584AE87,0xC006B4BADD886CAE // A4,A3 627// [32;64) 628data8 0x3F8CE54375841A5F,0xBFB801ABCFFA1BE2 // A10,A9 629data8 0x403040A8B1815BDA,0x405B99A917D24B7A // A2,A1 630data8 0x3FD30CAB81BFFA03,0xBFE41AEF61ECF48B // A8,A7 631data8 0x400650CC136BEC43,0xC016022046E8292B // A4,A3 632// [64;128) 633data8 0x3F9B69BD22CAA8B8,0xBFC6D48875B7A213 // A10,A9 634data8 0x40402028CCAA2F6D,0x40709AACEB3CBE0F // A2,A1 635data8 0x3FE22C6A5924761E,0xBFF342F5F224523D // A8,A7 636data8 0x4015CD405CCA331F,0xC025AAD10482C769 // A4,A3 637// [128;256) 638data8 0x3FAAAD9CD0E40D06,0xBFD63FC8505D80CB // A10,A9 639data8 0x40501008D56C2648,0x408364794B0F4376 // A2,A1 640data8 0x3FF1BE0126E00284,0xC002D8E3F6F7F7CA // A8,A7 641data8 0x40258C757E95D860,0xC0357FA8FD398011 // A4,A3 642// [256;512) 643data8 0x3FBA4DAC59D49FEB,0xBFE5F476D1C43A77 // A10,A9 644data8 0x40600800D890C7C6,0x40962C42AAEC8EF0 // A2,A1 645data8 0x40018680ECF19B89,0xC012A3EB96FB7BA4 // A8,A7 646data8 0x40356C4CDD3B60F9,0xC0456A34BF18F440 // A4,A3 647// [512;1024) 648data8 0x3FCA1B54F6225A5A,0xBFF5CD67BA10E048 // A10,A9 649data8 0x407003FED94C58C2,0x40A8F30B4ACBCD22 // A2,A1 650data8 0x40116A135EB66D8C,0xC022891B1CED527E // A8,A7 651data8 0x40455C4617FDD8BC,0xC0555F82729E59C4 // A4,A3 652// [1024;2048) 653data8 0x3FD9FFF9095C6EC9,0xC005B88CB25D76C9 // A10,A9 654data8 0x408001FE58FA734D,0x40BBB953BAABB0F3 // A2,A1 655data8 0x40215B2F9FEB5D87,0xC0327B539DEA5058 // A8,A7 656data8 0x40555444B3E8D64D,0xC0655A2B26F9FC8A // A4,A3 657// [2048;4096) 658data8 0x3FE9F065A1C3D6B1,0xC015ACF6FAE8D78D // A10,A9 659data8 0x409000FE383DD2B7,0x40CE7F5C1E8BCB8B // A2,A1 660data8 0x40315324E5DB2EBE,0xC04274194EF70D18 // A8,A7 661data8 0x4065504353FF2207,0xC075577FE1BFE7B6 // A4,A3 662// [4096;8192) 663data8 0x3FF9E6FBC6B1C70D,0xC025A62DAF76F85D // A10,A9 664data8 0x40A0007E2F61EBE8,0x40E0A2A23FB5F6C3 // A2,A1 665data8 0x40414E9BC0A0141A,0xC0527030F2B69D43 // A8,A7 666data8 0x40754E417717B45B,0xC085562A447258E5 // A4,A3 667// 668data8 0xbfdffffffffaea15 // P1 669data8 0x3FDD8B618D5AF8FE // point of local minimum on [1;2] 670data8 0x3FED67F1C864BEB5 // ln(sqrt(2*Pi)) 671data8 0x4008000000000000 // 3.0 672// 673data8 0xBF9E1C289FB224AB,0x3FBF7422445C9460 // A6,A5 674data8 0xBFF01E76D66F8D8A // A0 675data8 0xBFE2788CFC6F91DA // A1 [1.0;1.25) 676data8 0x3FCB8CC69000EB5C,0xBFD41997A0C2C641 // A6,A5 677data8 0x3FFCAB0BFA0EA462 // A0 678data8 0xBFBF19B9BCC38A42 // A0 [1.25;1.5) 679data8 0x3FD51EE4DE0A364C,0xBFE00D7F98A16E4B // A6,A5 680data8 0x40210CE1F327E9E4 // A0 681data8 0x4001DB08F9DFA0CC // A0 [1.5;1.75) 682data8 0x3FE24F606742D252,0xBFEC81D7D12574EC // A6,A5 683data8 0x403BE636A63A9C27 // A0 684data8 0x4000A0CB38D6CF0A // A0 [1.75;2.0) 685data8 0x3FF1029A9DD542B4,0xBFFAD37C209D3B25 // A6,A5 686data8 0x405385E6FD9BE7EA // A0 687data8 0x478895F1C0000000 // Overflow boundary 688data8 0x400062D97D26B523,0xC00A03E1529FF023 // A6,A5 689data8 0x4069204C51E566CE // A0 690data8 0x0000000000000000 // pad 691data8 0x40101476B38FD501,0xC0199DE7B387C0FC // A6,A5 692data8 0x407EB8DAEC83D759 // A0 693data8 0x0000000000000000 // pad 694data8 0x401FDB008D65125A,0xC0296B506E665581 // A6,A5 695data8 0x409226D93107EF66 // A0 696data8 0x0000000000000000 // pad 697data8 0x402FB3EAAF3E7B2D,0xC039521142AD8E0D // A6,A5 698data8 0x40A4EFA4F072792E // A0 699data8 0x0000000000000000 // pad 700data8 0x403FA024C66B2563,0xC0494569F250E691 // A6,A5 701data8 0x40B7B747C9235BB8 // A0 702data8 0x0000000000000000 // pad 703data8 0x404F9607D6DA512C,0xC0593F0B2EDDB4BC // A6,A5 704data8 0x40CA7E29C5F16DE2 // A0 705data8 0x0000000000000000 // pad 706data8 0x405F90C5F613D98D,0xC0693BD130E50AAF // A6,A5 707data8 0x40DD4495238B190C // A0 708data8 0x0000000000000000 // pad 709// 710// polynomial approximation of ln(sin(Pi*x)/(Pi*x)), |x| <= 0.5 711data8 0xBFD58731A486E820,0xBFA4452CC28E15A9 // S16,S14 712data8 0xBFD013F6E1B86C4F,0xBFD5B3F19F7A341F // S8,S6 713data8 0xBFC86A0D5252E778,0xBFC93E08C9EE284B // S12,S10 714data8 0xBFE15132555C9EDD,0xBFFA51A662480E35 // S4,S2 715// 716// [1.0;1.25) 717data8 0xBFA697D6775F48EA,0x3FB9894B682A98E7 // A9,A8 718data8 0xBFCA8969253CFF55,0x3FD15124EFB35D9D // A5,A4 719data8 0xBFC1B00158AB719D,0x3FC5997D04E7F1C1 // A7,A6 720data8 0xBFD9A4D50BAFF989,0x3FEA51A661F5176A // A3,A2 721// [1.25;1.5) 722data8 0x3F838E0D35A6171A,0xBF831BBBD61313B7 // A8,A7 723data8 0x3FB08B40196425D0,0xBFC2E427A53EB830 // A4,A3 724data8 0x3F9285DDDC20D6C3,0xBFA0C90C9C223044 // A6,A5 725data8 0x3FDEF72BC8F5287C,0x3D890B3DAEBC1DFC // A2,A1 726// [1.5;1.75) 727data8 0x3F65D5A7EB31047F,0xBFA44EAC9BFA7FDE // A8,A7 728data8 0x40051FEFE7A663D8,0xC012A5CFE00A2522 // A4,A3 729data8 0x3FD0E1583AB00E08,0xBFF084AF95883BA5 // A6,A5 730data8 0x40185982877AE0A2,0xC015F83DB73B57B7 // A2,A1 731// [1.75;2.0) 732data8 0x3F4A9222032EB39A,0xBF8CBC9587EEA5A3 // A8,A7 733data8 0x3FF795400783BE49,0xC00851BC418B8A25 // A4,A3 734data8 0x3FBBC992783E8C5B,0xBFDFA67E65E89B29 // A6,A5 735data8 0x4012B408F02FAF88,0xC013284CE7CB0C39 // A2,A1 736// 737// roots 738data8 0xC003A7FC9600F86C // -2.4570247382208005860 739data8 0xC009260DBC9E59AF // -3.1435808883499798405 740data8 0xC005FB410A1BD901 // -2.7476826467274126919 741data8 0xC00FA471547C2FE5 // -3.9552942848585979085 742// 743// polynomial approximation of ln(GAMMA(x)) near roots 744// near -2.4570247382208005860 745data8 0x3FF694A6058D9592,0x40136EEBB003A92B // R3,R2 746data8 0x3FF83FE966AF5360,0x3C90323B6D1FE86D // R1,R0 747// near -3.1435808883499798405 748data8 0x405C11371268DA38,0x4039D4D2977D2C23 // R3,R2 749data8 0x401F20A65F2FAC62,0x3CDE9605E3AE7A62 // R1,R0 750// near -2.7476826467274126919 751data8 0xC034185AC31314FF,0x4023267F3C28DFE3 // R3,R2 752data8 0xBFFEA12DA904B194,0x3CA8FB8530BA7689 // R1,R0 753// near -2.7476826467274126919 754data8 0xC0AD25359E70C888,0x406F76DEAEA1B8C6 // R3,R2 755data8 0xC034B99D966C5644,0xBCBDDC0336980B58 // R1,R0 756LOCAL_OBJECT_END(lgammaf_data) 757 758//********************************************************************* 759 760.section .text 761GLOBAL_LIBM_ENTRY(__libm_lgammaf) 762{ .mfi 763 getf.exp GR_SignExp = f8 764 frcpa.s1 FR_InvX,p0 = f1,f8 765 mov GR_ExpOf2 = 0x10000 766} 767{ .mfi 768 addl GR_ad_Data = @ltoff(lgammaf_data),gp 769 fcvt.fx.s1 FR_int_N = f8 770 mov GR_ExpMask = 0x1ffff 771};; 772{ .mfi 773 getf.sig GR_Sig = f8 774 fclass.m p13,p0 = f8,0x1EF // is x NaTVal, NaN, 775 // +/-0, +/-INF or +/-deno? 776 mov GR_ExpBias = 0xffff 777} 778{ .mfi 779 ld8 GR_ad_Data = [GR_ad_Data] 780 fma.s1 FR_Xp1 = f8,f1,f1 781 mov GR_StirlBound = 0x1000C 782};; 783{ .mfi 784 setf.exp FR_2 = GR_ExpOf2 785 fmerge.se FR_x = f1,f8 786 dep.z GR_Ind = GR_SignExp,3,4 787} 788{ .mfi 789 cmp.eq p8,p0 = GR_SignExp,GR_ExpBias 790 fcvt.fx.trunc.s1 FR_int_Ntrunc = f8 791 and GR_Exp = GR_ExpMask,GR_SignExp 792};; 793{ .mfi 794 add GR_ad_C650 = 0xB20,GR_ad_Data 795 fcmp.lt.s1 p14,p15 = f8,f0 796 extr.u GR_Ind4T = GR_Sig,55,8 797} 798{ .mfb 799 sub GR_PureExp = GR_Exp,GR_ExpBias 800 fnorm.s1 FR_NormX = f8 801 // jump if x is NaTVal, NaN, +/-0, +/-INF or +/-deno 802(p13) br.cond.spnt lgammaf_spec 803};; 804lgammaf_core: 805{ .mfi 806 ldfpd FR_P1,FR_LocalMin = [GR_ad_C650],16 807 fms.s1 FR_xm2 = f8,f1,f1 808 add GR_ad_Co = 0x820,GR_ad_Data 809} 810{ .mib 811 ldfpd FR_P3,FR_P2 = [GR_ad_Data],16 812 cmp.ltu p9,p0 = GR_SignExp,GR_ExpBias 813 // jump if x is from the interval [1; 2) 814(p8) br.cond.spnt lgammaf_1_2 815};; 816{ .mfi 817 setf.sig FR_int_Ln = GR_PureExp 818 fms.s1 FR_r = FR_InvX,f8,f1 819 shladd GR_ad_Co = GR_Ind,3,GR_ad_Co 820} 821{ .mib 822 ldfpd FR_LnSqrt2Pi,FR_3 = [GR_ad_C650],16 823 cmp.lt p13,p12 = GR_Exp,GR_StirlBound 824 // jump if x is from the interval (0; 1) 825(p9) br.cond.spnt lgammaf_0_1 826};; 827{ .mfi 828 ldfpd FR_Ln2,FR_05 = [GR_ad_Data],16 829 fma.s1 FR_Xp2 = f1,f1,FR_Xp1 // (x+2) 830 shladd GR_ad_C650 = GR_Ind,2,GR_ad_C650 831} 832{ .mfi 833 add GR_ad_Ce = 0x20,GR_ad_Co 834 nop.f 0 835 add GR_ad_C43 = 0x30,GR_ad_Co 836};; 837{ .mfi 838 // load coefficients of polynomial approximation 839 // of ln(GAMMA(x)), 2 <= x < 2^13 840(p13) ldfpd FR_A10,FR_A9 = [GR_ad_Co],16 841 fcvt.xf FR_N = FR_int_N 842 cmp.eq.unc p6,p7 = GR_ExpOf2,GR_SignExp 843} 844{ .mib 845(p13) ldfpd FR_A8,FR_A7 = [GR_ad_Ce] 846(p14) cmp.le.unc p9,p0 = GR_StirlBound,GR_Exp 847 // jump if x is less or equal to -2^13 848(p9) br.cond.spnt lgammaf_negstirling 849};; 850.pred.rel "mutex",p6,p7 851{ .mfi 852(p13) ldfpd FR_A6,FR_A5 = [GR_ad_C650],16 853(p6) fma.s1 FR_x = f0,f0,FR_NormX 854 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 855} 856{ .mfi 857(p13) ldfpd FR_A4,FR_A3 = [GR_ad_C43] 858(p7) fms.s1 FR_x = FR_x,f1,f1 859(p14) mov GR_ReqBound = 0x20005 860};; 861{ .mfi 862(p13) ldfpd FR_A2,FR_A1 = [GR_ad_Co],16 863 fms.s1 FR_xm2 = FR_xm2,f1,f1 864(p14) extr.u GR_Arg = GR_Sig,60,4 865} 866{ .mfi 867 mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1 868 fcvt.xf FR_Ntrunc = FR_int_Ntrunc 869 nop.i 0 870};; 871{ .mfi 872 ldfd FR_T = [GR_ad_T] 873 fma.s1 FR_r2 = FR_r,FR_r,f0 874 shl GR_ReqBound = GR_ReqBound,3 875} 876{ .mfi 877 add GR_ad_Co = 0xCA0,GR_ad_Data 878 fnma.s1 FR_Req = FR_Xp1,FR_NormX,f0 // -x*(x+1) 879(p14) shladd GR_Arg = GR_Exp,4,GR_Arg 880};; 881{ .mfi 882(p13) ldfd FR_A0 = [GR_ad_C650] 883 fma.s1 FR_Xp3 = FR_2,f1,FR_Xp1 // (x+3) 884(p14) cmp.le.unc p9,p0 = GR_Arg,GR_ReqBound 885} 886{ .mfi 887(p14) add GR_ad_Ce = 0x20,GR_ad_Co 888 fma.s1 FR_Xp4 = FR_2,FR_2,FR_NormX // (x+4) 889(p15) add GR_ad_OvfBound = 0xBB8,GR_ad_Data 890};; 891{ .mfi 892 // load coefficients of polynomial approximation 893 // of ln(sin(Pi*xf)/(Pi*xf)), |xf| <= 0.5 894(p14) ldfpd FR_S16,FR_S14 = [GR_ad_Co],16 895(p14) fms.s1 FR_Xf = FR_NormX,f1,FR_N // xf = x - [x] 896(p14) sub GR_SignOfGamma = r0,GR_SignOfGamma // set sign of 897 // gamma(x) to -1 898} 899{ .mfb 900(p14) ldfpd FR_S12,FR_S10 = [GR_ad_Ce],16 901 fma.s1 FR_Xp5 = FR_2,FR_2,FR_Xp1 // (x+5) 902 // jump if x is from the interval (-9; 0) 903(p9) br.cond.spnt lgammaf_negrecursion 904};; 905{ .mfi 906(p14) ldfpd FR_S8,FR_S6 = [GR_ad_Co],16 907 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2 908 nop.i 0 909} 910{ .mfb 911(p14) ldfpd FR_S4,FR_S2 = [GR_ad_Ce],16 912 fma.s1 FR_x2 = FR_x,FR_x,f0 913 // jump if x is from the interval (-2^13; -9) 914(p14) br.cond.spnt lgammaf_negpoly 915};; 916{ .mfi 917 ldfd FR_OverflowBound = [GR_ad_OvfBound] 918(p12) fcvt.xf FR_N = FR_int_Ln 919 // set p9 if signgum is 32-bit int 920 // set p10 if signgum is 64-bit int 921 cmp.eq p10,p9 = 8,r34 922} 923{ .mfi 924 nop.m 0 925(p12) fma.s1 FR_P10 = FR_P1,FR_r,f1 926 nop.i 0 927};; 928.pred.rel "mutex",p6,p7 929.pred.rel "mutex",p9,p10 930{ .mfi 931 // store sign of gamma(x) as 32-bit int 932(p9) st4 [r33] = GR_SignOfGamma 933(p6) fma.s1 FR_xx = FR_x,FR_xm2,f0 934 nop.i 0 935} 936{ .mfi 937 // store sign of gamma(x) as 64-bit int 938(p10) st8 [r33] = GR_SignOfGamma 939(p7) fma.s1 FR_xx = f0,f0,FR_x 940 nop.i 0 941};; 942{ .mfi 943 nop.m 0 944(p13) fma.s1 FR_A9 = FR_A10,FR_x,FR_A9 945 nop.i 0 946} 947{ .mfi 948 nop.m 0 949(p13) fma.s1 FR_A7 = FR_A8,FR_x,FR_A7 950 nop.i 0 951};; 952{ .mfi 953 nop.m 0 954(p13) fma.s1 FR_A5 = FR_A6,FR_x,FR_A5 955 nop.i 0 956} 957{ .mfi 958 nop.m 0 959(p13) fma.s1 FR_A3 = FR_A4,FR_x,FR_A3 960 nop.i 0 961};; 962{ .mfi 963 nop.m 0 964(p15) fcmp.eq.unc.s1 p8,p0 = FR_NormX,FR_2 // is input argument 2.0? 965 nop.i 0 966} 967{ .mfi 968 nop.m 0 969(p13) fma.s1 FR_A1 = FR_A2,FR_x,FR_A1 970 nop.i 0 971};; 972{ .mfi 973 nop.m 0 974(p12) fma.s1 FR_T = FR_N,FR_Ln2,FR_T 975 nop.i 0 976} 977{ .mfi 978 nop.m 0 979(p12) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10 980 nop.i 0 981};; 982{ .mfi 983 nop.m 0 984(p13) fma.s1 FR_x4 = FR_x2,FR_x2,f0 985 nop.i 0 986} 987{ .mfi 988 nop.m 0 989(p13) fma.s1 FR_x3 = FR_x2,FR_xx,f0 990 nop.i 0 991};; 992{ .mfi 993 nop.m 0 994(p13) fma.s1 FR_A7 = FR_A9,FR_x2,FR_A7 995 nop.i 0 996} 997{ .mfb 998 nop.m 0 999(p8) fma.s.s0 f8 = f0,f0,f0 1000(p8) br.ret.spnt b0 // fast exit for 2.0 1001};; 1002{ .mfi 1003 nop.m 0 1004(p6) fma.s1 FR_A0 = FR_A0,FR_xm2,f0 1005 nop.i 0 1006} 1007{ .mfi 1008 nop.m 0 1009(p13) fma.s1 FR_A3 = FR_A5,FR_x2,FR_A3 1010 nop.i 0 1011};; 1012{ .mfi 1013 nop.m 0 1014(p15) fcmp.le.unc.s1 p8,p0 = FR_OverflowBound,FR_NormX // overflow test 1015 nop.i 0 1016} 1017{ .mfi 1018 nop.m 0 1019(p12) fms.s1 FR_xm05 = FR_NormX,f1,FR_05 1020 nop.i 0 1021};; 1022{ .mfi 1023 nop.m 0 1024(p12) fma.s1 FR_Ln = FR_P32,FR_r,FR_T 1025 nop.i 0 1026} 1027{ .mfi 1028 nop.m 0 1029(p12) fms.s1 FR_LnSqrt2Pi = FR_LnSqrt2Pi,f1,FR_NormX 1030 nop.i 0 1031};; 1032{ .mfi 1033 nop.m 0 1034(p13) fma.s1 FR_A0 = FR_A1,FR_xx,FR_A0 1035 nop.i 0 1036} 1037{ .mfb 1038 nop.m 0 1039(p13) fma.s1 FR_A3 = FR_A7,FR_x4,FR_A3 1040 // jump if result overflows 1041(p8) br.cond.spnt lgammaf_overflow 1042};; 1043.pred.rel "mutex",p12,p13 1044{ .mfi 1045 nop.m 0 1046(p12) fma.s.s0 f8 = FR_Ln,FR_xm05,FR_LnSqrt2Pi 1047 nop.i 0 1048} 1049{ .mfb 1050 nop.m 0 1051(p13) fma.s.s0 f8 = FR_A3,FR_x3,FR_A0 1052 br.ret.sptk b0 1053};; 1054// branch for calculating of ln(GAMMA(x)) for 0 < x < 1 1055//--------------------------------------------------------------------- 1056.align 32 1057lgammaf_0_1: 1058{ .mfi 1059 getf.sig GR_Ind = FR_Xp1 1060 fma.s1 FR_r2 = FR_r,FR_r,f0 1061 mov GR_fff7 = 0xFFF7 1062} 1063{ .mfi 1064 ldfpd FR_Ln2,FR_05 = [GR_ad_Data],16 1065 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2 1066 // input argument can't be equal to 1.0 1067 cmp.eq p0,p14 = r0,r0 1068};; 1069{ .mfi 1070 getf.exp GR_Exp = FR_w 1071 fcvt.xf FR_N = FR_int_Ln 1072 add GR_ad_Co = 0xCE0,GR_ad_Data 1073} 1074{ .mfi 1075 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 1076 fma.s1 FR_P10 = FR_P1,FR_r,f1 1077 add GR_ad_Ce = 0xD00,GR_ad_Data 1078};; 1079{ .mfi 1080 ldfd FR_T = [GR_ad_T] 1081 fma.s1 FR_w2 = FR_w,FR_w,f0 1082 extr.u GR_Ind = GR_Ind,61,2 1083} 1084{ .mfi 1085 nop.m 0 1086 fma.s1 FR_Q32 = FR_P3,FR_w,FR_P2 1087//// add GR_ad_C0 = 0xB30,GR_ad_Data 1088 add GR_ad_C0 = 0xB38,GR_ad_Data 1089};; 1090{ .mfi 1091 and GR_Exp = GR_Exp,GR_ExpMask 1092 nop.f 0 1093 shladd GR_IndX8 = GR_Ind,3,r0 1094} 1095{ .mfi 1096 shladd GR_IndX2 = GR_Ind,1,r0 1097 fma.s1 FR_Q10 = FR_P1,FR_w,f1 1098 cmp.eq p6,p15 = 0,GR_Ind 1099};; 1100{ .mfi 1101 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co 1102(p6) fma.s1 FR_x = f0,f0,FR_NormX 1103 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0 1104} 1105{ .mfi 1106 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce 1107 nop.f 0 1108(p15) cmp.eq.unc p7,p8 = 1,GR_Ind 1109};; 1110.pred.rel "mutex",p7,p8 1111{ .mfi 1112 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16 1113(p7) fms.s1 FR_x = FR_NormX,f1,FR_LocalMin 1114 cmp.ge p10,p11 = GR_Exp,GR_fff7 1115} 1116{ .mfb 1117 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16 1118(p8) fma.s1 FR_x = f1,f1,FR_NormX 1119 br.cond.sptk lgamma_0_2_core 1120};; 1121// branch for calculating of ln(GAMMA(x)) for 1 <= x < 2 1122//--------------------------------------------------------------------- 1123.align 32 1124lgammaf_1_2: 1125{ .mfi 1126 add GR_ad_Co = 0xCF0,GR_ad_Data 1127 fcmp.eq.s1 p14,p0 = f1,FR_NormX // is input argument 1.0? 1128 extr.u GR_Ind = GR_Sig,61,2 1129} 1130{ .mfi 1131 add GR_ad_Ce = 0xD10,GR_ad_Data 1132 nop.f 0 1133//// add GR_ad_C0 = 0xB40,GR_ad_Data 1134 add GR_ad_C0 = 0xB48,GR_ad_Data 1135};; 1136{ .mfi 1137 shladd GR_IndX8 = GR_Ind,3,r0 1138 nop.f 0 1139 shladd GR_IndX2 = GR_Ind,1,r0 1140} 1141{ .mfi 1142 cmp.eq p6,p15 = 0,GR_Ind // p6 <- x from [1;1.25) 1143 nop.f 0 1144 cmp.ne p9,p0 = r0,r0 1145};; 1146{ .mfi 1147 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co 1148(p6) fms.s1 FR_x = FR_NormX,f1,f1 // reduced x for [1;1.25) 1149 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0 1150} 1151{ .mfi 1152 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce 1153(p14) fma.s.s0 f8 = f0,f0,f0 1154(p15) cmp.eq.unc p7,p8 = 1,GR_Ind // p7 <- x from [1.25;1.5) 1155};; 1156.pred.rel "mutex",p7,p8 1157{ .mfi 1158 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16 1159(p7) fms.s1 FR_x = FR_xm2,f1,FR_LocalMin 1160 nop.i 0 1161} 1162{ .mfi 1163 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16 1164(p8) fma.s1 FR_x = f0,f0,FR_NormX 1165(p9) cmp.eq.unc p10,p11 = r0,r0 1166};; 1167lgamma_0_2_core: 1168{ .mmi 1169 ldfpd FR_A4,FR_A3 = [GR_ad_Co],16 1170 ldfpd FR_A2,FR_A1 = [GR_ad_Ce],16 1171 mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1 1172};; 1173{ .mfi 1174// add GR_ad_C0 = 8,GR_ad_C0 1175 ldfd FR_A0 = [GR_ad_C0] 1176 nop.f 0 1177 // set p13 if signgum is 32-bit int 1178 // set p15 if signgum is 64-bit int 1179 cmp.eq p15,p13 = 8,r34 1180};; 1181.pred.rel "mutex",p13,p15 1182{ .mmf 1183 // store sign of gamma(x) 1184(p13) st4 [r33] = GR_SignOfGamma // as 32-bit int 1185(p15) st8 [r33] = GR_SignOfGamma // as 64-bit int 1186(p11) fma.s1 FR_Q32 = FR_Q32,FR_w2,FR_Q10 1187};; 1188{ .mfb 1189 nop.m 0 1190(p10) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10 1191(p14) br.ret.spnt b0 // fast exit for 1.0 1192};; 1193{ .mfi 1194 nop.m 0 1195(p10) fma.s1 FR_T = FR_N,FR_Ln2,FR_T 1196 cmp.eq p6,p7 = 0,GR_Ind // p6 <- x from [1;1.25) 1197} 1198{ .mfi 1199 nop.m 0 1200 fma.s1 FR_x2 = FR_x,FR_x,f0 1201 cmp.eq p8,p0 = r0,r0 // set p8 to 1 that means we on [1;2] 1202};; 1203{ .mfi 1204 nop.m 0 1205(p11) fma.s1 FR_Ln = FR_Q32,FR_w,f0 1206 nop.i 0 1207} 1208{ .mfi 1209 nop.m 0 1210 nop.f 0 1211 nop.i 0 1212};; 1213.pred.rel "mutex",p6,p7 1214{ .mfi 1215 nop.m 0 1216(p6) fma.s1 FR_xx = f0,f0,FR_x 1217 nop.i 0 1218} 1219{ .mfi 1220 nop.m 0 1221(p7) fma.s1 FR_xx = f0,f0,f1 1222 nop.i 0 1223};; 1224{ .mfi 1225 nop.m 0 1226 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7 1227 nop.i 0 1228} 1229{ .mfi 1230 nop.m 0 1231 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5 1232(p9) cmp.ne p8,p0 = r0,r0 // set p8 to 0 that means we on [0;1] 1233};; 1234{ .mfi 1235 nop.m 0 1236 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3 1237 nop.i 0 1238} 1239{ .mfi 1240 nop.m 0 1241 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1 1242 nop.i 0 1243};; 1244{ .mfi 1245 nop.m 0 1246 fma.s1 FR_x4 = FR_x2,FR_x2,f0 1247 nop.i 0 1248} 1249{ .mfi 1250 nop.m 0 1251(p10) fma.s1 FR_Ln = FR_P32,FR_r,FR_T 1252 nop.i 0 1253};; 1254{ .mfi 1255 nop.m 0 1256 fma.s1 FR_A5 = FR_A7,FR_x2,FR_A5 1257 nop.i 0 1258} 1259{ .mfi 1260 nop.m 0 1261 fma.s1 FR_A1 = FR_A3,FR_x2,FR_A1 1262 nop.i 0 1263};; 1264.pred.rel "mutex",p9,p8 1265{ .mfi 1266 nop.m 0 1267(p9) fms.d.s1 FR_A0 = FR_A0,FR_xx,FR_Ln 1268 nop.i 0 1269} 1270{ .mfi 1271 nop.m 0 1272(p8) fms.s1 FR_A0 = FR_A0,FR_xx,f0 1273 nop.i 0 1274};; 1275{ .mfi 1276 nop.m 0 1277 fma.d.s1 FR_A1 = FR_A5,FR_x4,FR_A1 1278 nop.i 0 1279} 1280{ .mfi 1281 nop.m 0 1282 nop.f 0 1283 nop.i 0 1284};; 1285.pred.rel "mutex",p6,p7 1286{ .mfi 1287 nop.m 0 1288(p6) fma.s.s0 f8 = FR_A1,FR_x2,FR_A0 1289 nop.i 0 1290} 1291{ .mfb 1292 nop.m 0 1293(p7) fma.s.s0 f8 = FR_A1,FR_x,FR_A0 1294 br.ret.sptk b0 1295};; 1296// branch for calculating of ln(GAMMA(x)) for -9 < x < 1 1297//--------------------------------------------------------------------- 1298.align 32 1299lgammaf_negrecursion: 1300{ .mfi 1301 getf.sig GR_N = FR_int_Ntrunc 1302 fms.s1 FR_1pXf = FR_Xp2,f1,FR_Ntrunc // 1 + (x+1) - [x] 1303 mov GR_Neg2 = 2 1304} 1305{ .mfi 1306 add GR_ad_Co = 0xCE0,GR_ad_Data 1307 fms.s1 FR_Xf = FR_Xp1,f1,FR_Ntrunc // (x+1) - [x] 1308 mov GR_Neg4 = 4 1309};; 1310{ .mfi 1311 add GR_ad_Ce = 0xD00,GR_ad_Data 1312 fma.s1 FR_Xp6 = FR_2,FR_2,FR_Xp2 // (x+6) 1313 add GR_ad_C0 = 0xB30,GR_ad_Data 1314} 1315{ .mfi 1316 sub GR_Neg2 = r0,GR_Neg2 1317 fma.s1 FR_Xp7 = FR_2,FR_3,FR_Xp1 // (x+7) 1318 sub GR_Neg4 = r0,GR_Neg4 1319};; 1320{ .mfi 1321 cmp.ne p8,p0 = r0,GR_N 1322 fcmp.eq.s1 p13,p0 = FR_NormX,FR_Ntrunc 1323 and GR_IntNum = 0xF,GR_N 1324} 1325{ .mfi 1326 cmp.lt p6,p0 = GR_N,GR_Neg2 1327 fma.s1 FR_Xp8 = FR_2,FR_3,FR_Xp2 // (x+8) 1328 cmp.lt p7,p0 = GR_N,GR_Neg4 1329};; 1330{ .mfi 1331 getf.d GR_Arg = FR_NormX 1332(p6) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp3,f0 1333(p8) tbit.z.unc p14,p15 = GR_IntNum,0 1334} 1335{ .mfi 1336 sub GR_RootInd = 0xE,GR_IntNum 1337(p7) fma.s1 FR_Xp4 = FR_Xp4,FR_Xp5,f0 1338 add GR_ad_Root = 0xDE0,GR_ad_Data 1339};; 1340{ .mfi 1341 shladd GR_ad_Root = GR_RootInd,3,GR_ad_Root 1342 fms.s1 FR_x = FR_Xp1,f1,FR_Ntrunc // (x+1) - [x] 1343 nop.i 0 1344} 1345{ .mfb 1346 nop.m 0 1347 nop.f 0 1348(p13) br.cond.spnt lgammaf_singularity 1349};; 1350.pred.rel "mutex",p14,p15 1351{ .mfi 1352 cmp.gt p6,p0 = 0xA,GR_IntNum 1353(p14) fma.s1 FR_Req = FR_Req,FR_Xf,f0 1354 cmp.gt p7,p0 = 0xD,GR_IntNum 1355} 1356{ .mfi 1357(p15) mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1 1358(p15) fnma.s1 FR_Req = FR_Req,FR_Xf,f0 1359 cmp.leu p0,p13 = 2,GR_RootInd 1360};; 1361{ .mfi 1362 nop.m 0 1363(p6) fma.s1 FR_Xp6 = FR_Xp6,FR_Xp7,f0 1364(p13) add GR_ad_RootCo = 0xE00,GR_ad_Data 1365};; 1366{ .mfi 1367 nop.m 0 1368 fcmp.eq.s1 p12,p11 = FR_1pXf,FR_2 1369 nop.i 0 1370};; 1371{ .mfi 1372 getf.sig GR_Sig = FR_1pXf 1373 fcmp.le.s1 p9,p0 = FR_05,FR_Xf 1374 nop.i 0 1375} 1376{ .mfi 1377(p13) shladd GR_RootInd = GR_RootInd,4,r0 1378(p7) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp4,f0 1379(p8) cmp.gt.unc p10,p0 = 0x9,GR_IntNum 1380};; 1381.pred.rel "mutex",p11,p12 1382{ .mfi 1383 nop.m 0 1384(p10) fma.s1 FR_Req = FR_Req,FR_Xp8,f0 1385(p11) extr.u GR_Ind = GR_Sig,61,2 1386} 1387{ .mfi 1388(p13) add GR_RootInd = GR_RootInd,GR_RootInd 1389 nop.f 0 1390(p12) mov GR_Ind = 3 1391};; 1392{ .mfi 1393 shladd GR_IndX2 = GR_Ind,1,r0 1394 nop.f 0 1395 cmp.gt p14,p0 = 2,GR_Ind 1396} 1397{ .mfi 1398 shladd GR_IndX8 = GR_Ind,3,r0 1399 nop.f 0 1400 cmp.eq p6,p0 = 1,GR_Ind 1401};; 1402.pred.rel "mutex",p6,p9 1403{ .mfi 1404 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co 1405(p6) fms.s1 FR_x = FR_Xf,f1,FR_LocalMin 1406 cmp.gt p10,p0 = 0xB,GR_IntNum 1407} 1408{ .mfi 1409 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce 1410(p9) fma.s1 FR_x = f0,f0,FR_1pXf 1411 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0 1412};; 1413{ .mfi 1414 // load coefficients of polynomial approximation 1415 // of ln(GAMMA(x)), 1 <= x < 2 1416 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16 1417(p10) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp6,f0 1418 add GR_ad_C0 = 8,GR_ad_C0 1419} 1420{ .mfi 1421 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16 1422 nop.f 0 1423(p14) add GR_ad_Root = 0x10,GR_ad_Root 1424};; 1425{ .mfi 1426 ldfpd FR_A4,FR_A3 = [GR_ad_Co],16 1427 nop.f 0 1428 add GR_ad_RootCe = 0xE10,GR_ad_Data 1429} 1430{ .mfi 1431 ldfpd FR_A2,FR_A1 = [GR_ad_Ce],16 1432 nop.f 0 1433(p14) add GR_RootInd = 0x40,GR_RootInd 1434};; 1435{ .mmi 1436 ldfd FR_A0 = [GR_ad_C0] 1437(p13) add GR_ad_RootCo = GR_ad_RootCo,GR_RootInd 1438(p13) add GR_ad_RootCe = GR_ad_RootCe,GR_RootInd 1439};; 1440{ .mmi 1441(p13) ld8 GR_Root = [GR_ad_Root] 1442(p13) ldfd FR_Root = [GR_ad_Root] 1443 mov GR_ExpBias = 0xffff 1444};; 1445{ .mfi 1446 nop.m 0 1447 fma.s1 FR_x2 = FR_x,FR_x,f0 1448 nop.i 0 1449} 1450{ .mlx 1451(p8) cmp.gt.unc p10,p0 = 0xF,GR_IntNum 1452 movl GR_Dx = 0x000000014F8B588E 1453};; 1454{ .mfi 1455 // load coefficients of polynomial approximation 1456 // of ln(GAMMA(x)), x is close to one of negative roots 1457(p13) ldfpd FR_R3,FR_R2 = [GR_ad_RootCo] 1458 // arguments for logarithm 1459(p10) fma.s1 FR_Req = FR_Req,FR_Xp2,f0 1460 mov GR_ExpMask = 0x1ffff 1461} 1462{ .mfi 1463(p13) ldfpd FR_R1,FR_R0 = [GR_ad_RootCe] 1464 nop.f 0 1465 // set p9 if signgum is 32-bit int 1466 // set p8 if signgum is 64-bit int 1467 cmp.eq p8,p9 = 8,r34 1468};; 1469.pred.rel "mutex",p9,p8 1470{ .mfi 1471(p9) st4 [r33] = GR_SignOfGamma // as 32-bit int 1472 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7 1473(p13) sub GR_Root = GR_Arg,GR_Root 1474} 1475{ .mfi 1476(p8) st8 [r33] = GR_SignOfGamma // as 64-bit int 1477 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5 1478 nop.i 0 1479};; 1480{ .mfi 1481 nop.m 0 1482 fms.s1 FR_w = FR_Req,f1,f1 1483(p13) add GR_Root = GR_Root,GR_Dx 1484} 1485{ .mfi 1486 nop.m 0 1487 nop.f 0 1488(p13) add GR_2xDx = GR_Dx,GR_Dx 1489};; 1490{ .mfi 1491 nop.m 0 1492 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3 1493 nop.i 0 1494} 1495{ .mfi 1496 nop.m 0 1497 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1 1498(p13) cmp.leu.unc p10,p0 = GR_Root,GR_2xDx 1499};; 1500{ .mfi 1501 nop.m 0 1502 frcpa.s1 FR_InvX,p0 = f1,FR_Req 1503 nop.i 0 1504} 1505{ .mfi 1506 nop.m 0 1507(p10) fms.s1 FR_rx = FR_NormX,f1,FR_Root 1508 nop.i 0 1509};; 1510{ .mfi 1511 getf.exp GR_SignExp = FR_Req 1512 fma.s1 FR_x4 = FR_x2,FR_x2,f0 1513 nop.i 0 1514};; 1515{ .mfi 1516 getf.sig GR_Sig = FR_Req 1517 fma.s1 FR_A5 = FR_A7,FR_x2,FR_A5 1518 nop.i 0 1519};; 1520{ .mfi 1521 sub GR_PureExp = GR_SignExp,GR_ExpBias 1522 fma.s1 FR_w2 = FR_w,FR_w,f0 1523 nop.i 0 1524} 1525{ .mfi 1526 nop.m 0 1527 fma.s1 FR_Q32 = FR_P3,FR_w,FR_P2 1528 nop.i 0 1529};; 1530{ .mfi 1531 setf.sig FR_int_Ln = GR_PureExp 1532 fma.s1 FR_A1 = FR_A3,FR_x2,FR_A1 1533 extr.u GR_Ind4T = GR_Sig,55,8 1534} 1535{ .mfi 1536 nop.m 0 1537 fma.s1 FR_Q10 = FR_P1,FR_w,f1 1538 nop.i 0 1539};; 1540{ .mfi 1541 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 1542 fms.s1 FR_r = FR_InvX,FR_Req,f1 1543 nop.i 0 1544} 1545{ .mfi 1546 nop.m 0 1547(p10) fms.s1 FR_rx2 = FR_rx,FR_rx,f0 1548 nop.i 0 1549};; 1550{ .mfi 1551 ldfd FR_T = [GR_ad_T] 1552(p10) fma.s1 FR_R2 = FR_R3,FR_rx,FR_R2 1553 nop.i 0 1554} 1555{ .mfi 1556 nop.m 0 1557(p10) fma.s1 FR_R0 = FR_R1,FR_rx,FR_R0 1558 nop.i 0 1559};; 1560{ .mfi 1561 getf.exp GR_Exp = FR_w 1562 fma.s1 FR_A1 = FR_A5,FR_x4,FR_A1 1563 mov GR_ExpMask = 0x1ffff 1564} 1565{ .mfi 1566 nop.m 0 1567 fma.s1 FR_Q32 = FR_Q32, FR_w2,FR_Q10 1568 nop.i 0 1569};; 1570{ .mfi 1571 nop.m 0 1572 fma.s1 FR_r2 = FR_r,FR_r,f0 1573 mov GR_fff7 = 0xFFF7 1574} 1575{ .mfi 1576 nop.m 0 1577 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2 1578 nop.i 0 1579};; 1580{ .mfi 1581 nop.m 0 1582 fma.s1 FR_P10 = FR_P1,FR_r,f1 1583 and GR_Exp = GR_ExpMask,GR_Exp 1584} 1585{ .mfb 1586 nop.m 0 1587(p10) fma.s.s0 f8 = FR_R2,FR_rx2,FR_R0 1588(p10) br.ret.spnt b0 // exit for arguments close to negative roots 1589};; 1590{ .mfi 1591 nop.m 0 1592 fcvt.xf FR_N = FR_int_Ln 1593 nop.i 0 1594} 1595{ .mfi 1596 cmp.ge p14,p15 = GR_Exp,GR_fff7 1597 nop.f 0 1598 nop.i 0 1599};; 1600{ .mfi 1601 nop.m 0 1602 fma.s1 FR_A0 = FR_A1,FR_x,FR_A0 1603 nop.i 0 1604} 1605{ .mfi 1606 nop.m 0 1607(p15) fma.s1 FR_Ln = FR_Q32,FR_w,f0 1608 nop.i 0 1609};; 1610{ .mfi 1611 nop.m 0 1612(p14) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10 1613 cmp.eq p6,p7 = 0,GR_Ind 1614};; 1615{ .mfi 1616 nop.m 0 1617(p14) fma.s1 FR_T = FR_N,FR_Ln2,FR_T 1618 nop.i 0 1619};; 1620{ .mfi 1621 nop.m 0 1622(p14) fma.s1 FR_Ln = FR_P32,FR_r,FR_T 1623 nop.i 0 1624};; 1625.pred.rel "mutex",p6,p7 1626{ .mfi 1627 nop.m 0 1628(p6) fms.s.s0 f8 = FR_A0,FR_x,FR_Ln 1629 nop.i 0 1630} 1631{ .mfb 1632 nop.m 0 1633(p7) fms.s.s0 f8 = FR_A0,f1,FR_Ln 1634 br.ret.sptk b0 1635};; 1636 1637// branch for calculating of ln(GAMMA(x)) for x < -2^13 1638//--------------------------------------------------------------------- 1639.align 32 1640lgammaf_negstirling: 1641{ .mfi 1642 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 1643 fms.s1 FR_Xf = FR_NormX,f1,FR_N // xf = x - [x] 1644 mov GR_SingBound = 0x10016 1645} 1646{ .mfi 1647 add GR_ad_Co = 0xCA0,GR_ad_Data 1648 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2 1649 nop.i 0 1650};; 1651{ .mfi 1652 ldfd FR_T = [GR_ad_T] 1653 fcvt.xf FR_int_Ln = FR_int_Ln 1654 cmp.le p6,p0 = GR_SingBound,GR_Exp 1655} 1656{ .mfb 1657 add GR_ad_Ce = 0x20,GR_ad_Co 1658 fma.s1 FR_r2 = FR_r,FR_r,f0 1659(p6) br.cond.spnt lgammaf_singularity 1660};; 1661{ .mfi 1662 // load coefficients of polynomial approximation 1663 // of ln(sin(Pi*xf)/(Pi*xf)), |xf| <= 0.5 1664 ldfpd FR_S16,FR_S14 = [GR_ad_Co],16 1665 fma.s1 FR_P10 = FR_P1,FR_r,f1 1666 nop.i 0 1667} 1668{ .mfi 1669 ldfpd FR_S12,FR_S10 = [GR_ad_Ce],16 1670 fms.s1 FR_xm05 = FR_NormX,f1,FR_05 1671 nop.i 0 1672};; 1673{ .mmi 1674 ldfpd FR_S8,FR_S6 = [GR_ad_Co],16 1675 ldfpd FR_S4,FR_S2 = [GR_ad_Ce],16 1676 nop.i 0 1677};; 1678{ .mfi 1679 getf.sig GR_N = FR_int_Ntrunc // signgam calculation 1680 fma.s1 FR_Xf2 = FR_Xf,FR_Xf,f0 1681 nop.i 0 1682};; 1683{ .mfi 1684 nop.m 0 1685 frcpa.s1 FR_InvXf,p0 = f1,FR_Xf 1686 nop.i 0 1687};; 1688{ .mfi 1689 getf.d GR_Arg = FR_Xf 1690 fcmp.eq.s1 p6,p0 = FR_NormX,FR_N 1691 mov GR_ExpBias = 0x3FF 1692};; 1693{ .mfi 1694 nop.m 0 1695 fma.s1 FR_T = FR_int_Ln,FR_Ln2,FR_T 1696 extr.u GR_Exp = GR_Arg,52,11 1697} 1698{ .mfi 1699 nop.m 0 1700 fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10 1701 nop.i 0 1702};; 1703{ .mfi 1704 sub GR_PureExp = GR_Exp,GR_ExpBias 1705 fma.s1 FR_S14 = FR_S16,FR_Xf2,FR_S14 1706 extr.u GR_Ind4T = GR_Arg,44,8 1707} 1708{ .mfb 1709 mov GR_SignOfGamma = 1 // set signgam to -1 1710 fma.s1 FR_S10 = FR_S12,FR_Xf2,FR_S10 1711(p6) br.cond.spnt lgammaf_singularity 1712};; 1713{ .mfi 1714 setf.sig FR_int_Ln = GR_PureExp 1715 fms.s1 FR_rf = FR_InvXf,FR_Xf,f1 1716 // set p14 if GR_N is even 1717 tbit.z p14,p0 = GR_N,0 1718} 1719{ .mfi 1720 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 1721 fma.s1 FR_Xf4 = FR_Xf2,FR_Xf2,f0 1722 nop.i 0 1723};; 1724{ .mfi 1725(p14) sub GR_SignOfGamma = r0,GR_SignOfGamma // set signgam to -1 1726 fma.s1 FR_S6 = FR_S8,FR_Xf2,FR_S6 1727 nop.i 0 1728} 1729{ .mfi 1730 // set p9 if signgum is 32-bit int 1731 // set p10 if signgum is 64-bit int 1732 cmp.eq p10,p9 = 8,r34 1733 fma.s1 FR_S2 = FR_S4,FR_Xf2,FR_S2 1734 nop.i 0 1735};; 1736{ .mfi 1737 ldfd FR_Tf = [GR_ad_T] 1738 fma.s1 FR_Ln = FR_P32,FR_r,FR_T 1739 nop.i 0 1740} 1741{ .mfi 1742 nop.m 0 1743 fma.s1 FR_LnSqrt2Pi = FR_LnSqrt2Pi,f1,FR_NormX 1744 nop.i 0 1745};; 1746.pred.rel "mutex",p9,p10 1747{ .mfi 1748(p9) st4 [r33] = GR_SignOfGamma // as 32-bit int 1749 fma.s1 FR_rf2 = FR_rf,FR_rf,f0 1750 nop.i 0 1751} 1752{ .mfi 1753(p10) st8 [r33] = GR_SignOfGamma // as 64-bit int 1754 fma.s1 FR_S10 = FR_S14,FR_Xf4,FR_S10 1755 nop.i 0 1756};; 1757{ .mfi 1758 nop.m 0 1759 fma.s1 FR_P32f = FR_P3,FR_rf,FR_P2 1760 nop.i 0 1761} 1762{ .mfi 1763 nop.m 0 1764 fma.s1 FR_Xf8 = FR_Xf4,FR_Xf4,f0 1765 nop.i 0 1766};; 1767{ .mfi 1768 nop.m 0 1769 fma.s1 FR_P10f = FR_P1,FR_rf,f1 1770 nop.i 0 1771} 1772{ .mfi 1773 nop.m 0 1774 fma.s1 FR_S2 = FR_S6,FR_Xf4,FR_S2 1775 nop.i 0 1776};; 1777{ .mfi 1778 nop.m 0 1779 fms.s1 FR_Ln = FR_Ln,FR_xm05,FR_LnSqrt2Pi 1780 nop.i 0 1781};; 1782{ .mfi 1783 nop.m 0 1784 fcvt.xf FR_Nf = FR_int_Ln 1785 nop.i 0 1786};; 1787{ .mfi 1788 nop.m 0 1789 fma.s1 FR_S2 = FR_S10,FR_Xf8,FR_S2 1790 nop.i 0 1791};; 1792{ .mfi 1793 nop.m 0 1794 fma.s1 FR_Tf = FR_Nf,FR_Ln2,FR_Tf 1795 nop.i 0 1796} 1797{ .mfi 1798 nop.m 0 1799 fma.s1 FR_P32f = FR_P32f,FR_rf2,FR_P10f // ?????? 1800 nop.i 0 1801};; 1802{ .mfi 1803 nop.m 0 1804 fnma.s1 FR_Ln = FR_S2,FR_Xf2,FR_Ln 1805 nop.i 0 1806};; 1807{ .mfi 1808 nop.m 0 1809 fma.s1 FR_Lnf = FR_P32f,FR_rf,FR_Tf 1810 nop.i 0 1811};; 1812{ .mfb 1813 nop.m 0 1814 fms.s.s0 f8 = FR_Ln,f1,FR_Lnf 1815 br.ret.sptk b0 1816};; 1817// branch for calculating of ln(GAMMA(x)) for -2^13 < x < -9 1818//--------------------------------------------------------------------- 1819.align 32 1820lgammaf_negpoly: 1821{ .mfi 1822 getf.d GR_Arg = FR_Xf 1823 frcpa.s1 FR_InvXf,p0 = f1,FR_Xf 1824 mov GR_ExpBias = 0x3FF 1825} 1826{ .mfi 1827 nop.m 0 1828 fma.s1 FR_Xf2 = FR_Xf,FR_Xf,f0 1829 nop.i 0 1830};; 1831{ .mfi 1832 getf.sig GR_N = FR_int_Ntrunc 1833 fcvt.xf FR_N = FR_int_Ln 1834 mov GR_SignOfGamma = 1 1835} 1836{ .mfi 1837 nop.m 0 1838 fma.s1 FR_A9 = FR_A10,FR_x,FR_A9 1839 nop.i 0 1840};; 1841{ .mfi 1842 nop.m 0 1843 fma.s1 FR_P10 = FR_P1,FR_r,f1 1844 extr.u GR_Exp = GR_Arg,52,11 1845} 1846{ .mfi 1847 nop.m 0 1848 fma.s1 FR_x4 = FR_x2,FR_x2,f0 1849 nop.i 0 1850};; 1851{ .mfi 1852 sub GR_PureExp = GR_Exp,GR_ExpBias 1853 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7 1854 tbit.z p14,p0 = GR_N,0 1855} 1856{ .mfi 1857 nop.m 0 1858 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5 1859 nop.i 0 1860};; 1861{ .mfi 1862 setf.sig FR_int_Ln = GR_PureExp 1863 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3 1864 nop.i 0 1865} 1866{ .mfi 1867 nop.m 0 1868 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1 1869(p14) sub GR_SignOfGamma = r0,GR_SignOfGamma 1870};; 1871{ .mfi 1872 nop.m 0 1873 fms.s1 FR_rf = FR_InvXf,FR_Xf,f1 1874 nop.i 0 1875} 1876{ .mfi 1877 nop.m 0 1878 fma.s1 FR_Xf4 = FR_Xf2,FR_Xf2,f0 1879 nop.i 0 1880};; 1881{ .mfi 1882 nop.m 0 1883 fma.s1 FR_S14 = FR_S16,FR_Xf2,FR_S14 1884 nop.i 0 1885} 1886{ .mfi 1887 nop.m 0 1888 fma.s1 FR_S10 = FR_S12,FR_Xf2,FR_S10 1889 nop.i 0 1890};; 1891{ .mfi 1892 nop.m 0 1893 fma.s1 FR_T = FR_N,FR_Ln2,FR_T 1894 nop.i 0 1895} 1896{ .mfi 1897 nop.m 0 1898 fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10 1899 nop.i 0 1900};; 1901{ .mfi 1902 nop.m 0 1903 fma.s1 FR_S6 = FR_S8,FR_Xf2,FR_S6 1904 extr.u GR_Ind4T = GR_Arg,44,8 1905} 1906{ .mfi 1907 nop.m 0 1908 fma.s1 FR_S2 = FR_S4,FR_Xf2,FR_S2 1909 nop.i 0 1910};; 1911{ .mfi 1912 nop.m 0 1913 fma.s1 FR_A7 = FR_A9,FR_x2,FR_A7 1914 nop.i 0 1915} 1916{ .mfi 1917 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data 1918 fma.s1 FR_A3 = FR_A5,FR_x2,FR_A3 1919 nop.i 0 1920};; 1921{ .mfi 1922 nop.m 0 1923 fma.s1 FR_Xf8 = FR_Xf4,FR_Xf4,f0 1924 nop.i 0 1925} 1926{ .mfi 1927 nop.m 0 1928 fma.s1 FR_rf2 = FR_rf,FR_rf,f0 1929 nop.i 0 1930};; 1931{ .mfi 1932 nop.m 0 1933 fma.s1 FR_P32f = FR_P3,FR_rf,FR_P2 1934 nop.i 0 1935} 1936{ .mfi 1937 nop.m 0 1938 fma.s1 FR_P10f = FR_P1,FR_rf,f1 1939 nop.i 0 1940};; 1941{ .mfi 1942 ldfd FR_Tf = [GR_ad_T] 1943 fma.s1 FR_Ln = FR_P32,FR_r,FR_T 1944 nop.i 0 1945} 1946{ .mfi 1947 nop.m 0 1948 fma.s1 FR_A0 = FR_A1,FR_x,FR_A0 1949 nop.i 0 1950};; 1951{ .mfi 1952 nop.m 0 1953 fma.s1 FR_S10 = FR_S14,FR_Xf4,FR_S10 1954 nop.i 0 1955} 1956{ .mfi 1957 nop.m 0 1958 fma.s1 FR_S2 = FR_S6,FR_Xf4,FR_S2 1959 nop.i 0 1960};; 1961{ .mfi 1962 nop.m 0 1963 fcvt.xf FR_Nf = FR_int_Ln 1964 nop.i 0 1965} 1966{ .mfi 1967 nop.m 0 1968 fma.s1 FR_A3 = FR_A7,FR_x4,FR_A3 1969 nop.i 0 1970};; 1971{ .mfi 1972 nop.m 0 1973 fcmp.eq.s1 p13,p0 = FR_NormX,FR_Ntrunc 1974 nop.i 0 1975} 1976{ .mfi 1977 nop.m 0 1978 fnma.s1 FR_x3 = FR_x2,FR_x,f0 // -x^3 1979 nop.i 0 1980};; 1981{ .mfi 1982 nop.m 0 1983 fma.s1 FR_P32f = FR_P32f,FR_rf2,FR_P10f 1984 nop.i 0 1985};; 1986{ .mfb 1987 // set p9 if signgum is 32-bit int 1988 // set p10 if signgum is 64-bit int 1989 cmp.eq p10,p9 = 8,r34 1990 fma.s1 FR_S2 = FR_S10,FR_Xf8,FR_S2 1991(p13) br.cond.spnt lgammaf_singularity 1992};; 1993.pred.rel "mutex",p9,p10 1994{ .mmf 1995(p9) st4 [r33] = GR_SignOfGamma // as 32-bit int 1996(p10) st8 [r33] = GR_SignOfGamma // as 64-bit int 1997 fms.s1 FR_A0 = FR_A3,FR_x3,FR_A0 // -A3*x^3-A0 1998};; 1999{ .mfi 2000 nop.m 0 2001 fma.s1 FR_Tf = FR_Nf,FR_Ln2,FR_Tf 2002 nop.i 0 2003};; 2004{ .mfi 2005 nop.m 0 2006 fma.s1 FR_Ln = FR_S2,FR_Xf2,FR_Ln // S2*Xf^2+Ln 2007 nop.i 0 2008};; 2009{ .mfi 2010 nop.m 0 2011 fma.s1 FR_Lnf = FR_P32f,FR_rf,FR_Tf 2012 nop.i 0 2013};; 2014{ .mfi 2015 nop.m 0 2016 fms.s1 FR_Ln = FR_A0,f1,FR_Ln 2017 nop.i 0 2018};; 2019{ .mfb 2020 nop.m 0 2021 fms.s.s0 f8 = FR_Ln,f1,FR_Lnf 2022 br.ret.sptk b0 2023};; 2024// branch for handling +/-0, NaT, QNaN, +/-INF and denormalised numbers 2025//--------------------------------------------------------------------- 2026.align 32 2027lgammaf_spec: 2028{ .mfi 2029 getf.exp GR_SignExp = FR_NormX 2030 fclass.m p6,p0 = f8,0x21 // is arg +INF? 2031 mov GR_SignOfGamma = 1 // set signgam to 1 2032};; 2033{ .mfi 2034 getf.sig GR_Sig = FR_NormX 2035 fclass.m p7,p0 = f8,0xB // is x deno? 2036 // set p11 if signgum is 32-bit int 2037 // set p12 if signgum is 64-bit int 2038 cmp.eq p12,p11 = 8,r34 2039};; 2040.pred.rel "mutex",p11,p12 2041{ .mfi 2042 // store sign of gamma(x) as 32-bit int 2043(p11) st4 [r33] = GR_SignOfGamma 2044 fclass.m p8,p0 = f8,0x1C0 // is arg NaT or NaN? 2045 dep.z GR_Ind = GR_SignExp,3,4 2046} 2047{ .mib 2048 // store sign of gamma(x) as 64-bit int 2049(p12) st8 [r33] = GR_SignOfGamma 2050 and GR_Exp = GR_ExpMask,GR_SignExp 2051(p6) br.ret.spnt b0 // exit for +INF 2052};; 2053{ .mfi 2054 sub GR_PureExp = GR_Exp,GR_ExpBias 2055 fclass.m p9,p0 = f8,0x22 // is arg -INF? 2056 extr.u GR_Ind4T = GR_Sig,55,8 2057} 2058{ .mfb 2059 nop.m 0 2060(p7) fma.s0 FR_tmp = f1,f1,f8 2061(p7) br.cond.sptk lgammaf_core 2062};; 2063{ .mfb 2064 nop.m 0 2065(p8) fms.s.s0 f8 = f8,f1,f8 2066(p8) br.ret.spnt b0 // exit for NaT and NaN 2067};; 2068{ .mfb 2069 nop.m 0 2070(p9) fmerge.s f8 = f1,f8 2071(p9) br.ret.spnt b0 // exit -INF 2072};; 2073// branch for handling negative integers and +/-0 2074//--------------------------------------------------------------------- 2075.align 32 2076lgammaf_singularity: 2077{ .mfi 2078 mov GR_SignOfGamma = 1 // set signgam to 1 2079 fclass.m p6,p0 = f8,0x6 // is x -0? 2080 mov GR_TAG = 109 // negative 2081} 2082{ .mfi 2083 mov GR_ad_SignGam = r33 2084 fma.s1 FR_X = f0,f0,f8 2085 nop.i 0 2086};; 2087{ .mfi 2088 nop.m 0 2089 frcpa.s0 f8,p0 = f1,f0 2090 // set p9 if signgum is 32-bit int 2091 // set p10 if signgum is 64-bit int 2092 cmp.eq p10,p9 = 8,r34 2093} 2094{ .mib 2095 nop.m 0 2096(p6) sub GR_SignOfGamma = r0,GR_SignOfGamma 2097 br.cond.sptk lgammaf_libm_err 2098};; 2099// overflow (x > OVERFLOV_BOUNDARY) 2100//--------------------------------------------------------------------- 2101.align 32 2102lgammaf_overflow: 2103{ .mfi 2104 nop.m 0 2105 nop.f 0 2106 mov r8 = 0x1FFFE 2107};; 2108{ .mfi 2109 setf.exp f9 = r8 2110 fmerge.s FR_X = f8,f8 2111 mov GR_TAG = 108 // overflow 2112};; 2113{ .mfi 2114 mov GR_ad_SignGam = r33 2115 nop.f 0 2116 // set p9 if signgum is 32-bit int 2117 // set p10 if signgum is 64-bit int 2118 cmp.eq p10,p9 = 8,r34 2119} 2120{ .mfi 2121 nop.m 0 2122 fma.s.s0 f8 = f9,f9,f0 // Set I,O and +INF result 2123 nop.i 0 2124};; 2125// gate to __libm_error_support# 2126//--------------------------------------------------------------------- 2127.align 32 2128lgammaf_libm_err: 2129{ .mmi 2130 alloc r32 = ar.pfs,1,4,4,0 2131 mov GR_Parameter_TAG = GR_TAG 2132 nop.i 0 2133};; 2134.pred.rel "mutex",p9,p10 2135{ .mmi 2136 // store sign of gamma(x) as 32-bit int 2137(p9) st4 [GR_ad_SignGam] = GR_SignOfGamma 2138 // store sign of gamma(x) as 64-bit int 2139(p10) st8 [GR_ad_SignGam] = GR_SignOfGamma 2140 nop.i 0 2141};; 2142GLOBAL_LIBM_END(__libm_lgammaf) 2143 2144 2145LOCAL_LIBM_ENTRY(__libm_error_region) 2146.prologue 2147{ .mfi 2148 add GR_Parameter_Y=-32,sp // Parameter 2 value 2149 nop.f 0 2150.save ar.pfs,GR_SAVE_PFS 2151 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs 2152} 2153{ .mfi 2154.fframe 64 2155 add sp=-64,sp // Create new stack 2156 nop.f 0 2157 mov GR_SAVE_GP=gp // Save gp 2158};; 2159{ .mmi 2160 stfs [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack 2161 add GR_Parameter_X = 16,sp // Parameter 1 address 2162.save b0, GR_SAVE_B0 2163 mov GR_SAVE_B0=b0 // Save b0 2164};; 2165.body 2166{ .mib 2167 stfs [GR_Parameter_X] = FR_X // STORE Parameter 1 2168 // on stack 2169 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address 2170 nop.b 0 2171} 2172{ .mib 2173 stfs [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3 2174 // on stack 2175 add GR_Parameter_Y = -16,GR_Parameter_Y 2176 br.call.sptk b0=__libm_error_support# // Call error handling 2177 // function 2178};; 2179{ .mmi 2180 nop.m 0 2181 nop.m 0 2182 add GR_Parameter_RESULT = 48,sp 2183};; 2184{ .mmi 2185 ldfs f8 = [GR_Parameter_RESULT] // Get return result off stack 2186.restore sp 2187 add sp = 64,sp // Restore stack pointer 2188 mov b0 = GR_SAVE_B0 // Restore return address 2189};; 2190{ .mib 2191 mov gp = GR_SAVE_GP // Restore gp 2192 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs 2193 br.ret.sptk b0 // Return 2194};; 2195 2196LOCAL_LIBM_END(__libm_error_region) 2197.type __libm_error_support#,@function 2198.global __libm_error_support# 2199