1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis 22*25c28e83SPiotr Jasiukajtis /* 23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24*25c28e83SPiotr Jasiukajtis */ 25*25c28e83SPiotr Jasiukajtis /* 26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 28*25c28e83SPiotr Jasiukajtis */ 29*25c28e83SPiotr Jasiukajtis 30*25c28e83SPiotr Jasiukajtis #pragma weak tgammaf = __tgammaf 31*25c28e83SPiotr Jasiukajtis 32*25c28e83SPiotr Jasiukajtis /* 33*25c28e83SPiotr Jasiukajtis * True gamma function 34*25c28e83SPiotr Jasiukajtis * 35*25c28e83SPiotr Jasiukajtis * float tgammaf(float x) 36*25c28e83SPiotr Jasiukajtis * 37*25c28e83SPiotr Jasiukajtis * Algorithm: see tgamma.c 38*25c28e83SPiotr Jasiukajtis * 39*25c28e83SPiotr Jasiukajtis * Maximum error observed: 0.87ulp (both positive and negative arguments) 40*25c28e83SPiotr Jasiukajtis */ 41*25c28e83SPiotr Jasiukajtis 42*25c28e83SPiotr Jasiukajtis #include "libm.h" 43*25c28e83SPiotr Jasiukajtis #include "libm_synonyms.h" 44*25c28e83SPiotr Jasiukajtis #include <math.h> 45*25c28e83SPiotr Jasiukajtis #if defined(__SUNPRO_C) 46*25c28e83SPiotr Jasiukajtis #include <sunmath.h> 47*25c28e83SPiotr Jasiukajtis #endif 48*25c28e83SPiotr Jasiukajtis #include <sys/isa_defs.h> 49*25c28e83SPiotr Jasiukajtis 50*25c28e83SPiotr Jasiukajtis #if defined(_BIG_ENDIAN) 51*25c28e83SPiotr Jasiukajtis #define HIWORD 0 52*25c28e83SPiotr Jasiukajtis #define LOWORD 1 53*25c28e83SPiotr Jasiukajtis #else 54*25c28e83SPiotr Jasiukajtis #define HIWORD 1 55*25c28e83SPiotr Jasiukajtis #define LOWORD 0 56*25c28e83SPiotr Jasiukajtis #endif 57*25c28e83SPiotr Jasiukajtis #define __HI(x) ((int *) &x)[HIWORD] 58*25c28e83SPiotr Jasiukajtis #define __LO(x) ((unsigned *) &x)[LOWORD] 59*25c28e83SPiotr Jasiukajtis 60*25c28e83SPiotr Jasiukajtis /* Coefficients for primary intervals GTi() */ 61*25c28e83SPiotr Jasiukajtis static const double cr[] = { 62*25c28e83SPiotr Jasiukajtis /* p1 */ 63*25c28e83SPiotr Jasiukajtis +7.09087253435088360271451613398019280077561279443e-0001, 64*25c28e83SPiotr Jasiukajtis -5.17229560788652108545141978238701790105241761089e-0001, 65*25c28e83SPiotr Jasiukajtis +5.23403394528150789405825222323770647162337764327e-0001, 66*25c28e83SPiotr Jasiukajtis -4.54586308717075010784041566069480411732634814899e-0001, 67*25c28e83SPiotr Jasiukajtis +4.20596490915239085459964590559256913498190955233e-0001, 68*25c28e83SPiotr Jasiukajtis -3.57307589712377520978332185838241458642142185789e-0001, 69*25c28e83SPiotr Jasiukajtis 70*25c28e83SPiotr Jasiukajtis /* p2 */ 71*25c28e83SPiotr Jasiukajtis +4.28486983980295198166056119223984284434264344578e-0001, 72*25c28e83SPiotr Jasiukajtis -1.30704539487709138528680121627899735386650103914e-0001, 73*25c28e83SPiotr Jasiukajtis +1.60856285038051955072861219352655851542955430871e-0001, 74*25c28e83SPiotr Jasiukajtis -9.22285161346010583774458802067371182158937943507e-0002, 75*25c28e83SPiotr Jasiukajtis +7.19240511767225260740890292605070595560626179357e-0002, 76*25c28e83SPiotr Jasiukajtis -4.88158265593355093703112238534484636193260459574e-0002, 77*25c28e83SPiotr Jasiukajtis 78*25c28e83SPiotr Jasiukajtis /* p3 */ 79*25c28e83SPiotr Jasiukajtis +3.82409531118807759081121479786092134814808872880e-0001, 80*25c28e83SPiotr Jasiukajtis +2.65309888180188647956400403013495759365167853426e-0002, 81*25c28e83SPiotr Jasiukajtis +8.06815109775079171923561169415370309376296739835e-0002, 82*25c28e83SPiotr Jasiukajtis -1.54821591666137613928840890835174351674007764799e-0002, 83*25c28e83SPiotr Jasiukajtis +1.76308239242717268530498313416899188157165183405e-0002, 84*25c28e83SPiotr Jasiukajtis 85*25c28e83SPiotr Jasiukajtis /* GZi and TZi */ 86*25c28e83SPiotr Jasiukajtis +0.9382046279096824494097535615803269576988, /* GZ1 */ 87*25c28e83SPiotr Jasiukajtis +0.8856031944108887002788159005825887332080, /* GZ2 */ 88*25c28e83SPiotr Jasiukajtis +0.9367814114636523216188468970808378497426, /* GZ3 */ 89*25c28e83SPiotr Jasiukajtis -0.3517214357852935791015625, /* TZ1 */ 90*25c28e83SPiotr Jasiukajtis +0.280530631542205810546875, /* TZ3 */ 91*25c28e83SPiotr Jasiukajtis }; 92*25c28e83SPiotr Jasiukajtis 93*25c28e83SPiotr Jasiukajtis #define P10 cr[0] 94*25c28e83SPiotr Jasiukajtis #define P11 cr[1] 95*25c28e83SPiotr Jasiukajtis #define P12 cr[2] 96*25c28e83SPiotr Jasiukajtis #define P13 cr[3] 97*25c28e83SPiotr Jasiukajtis #define P14 cr[4] 98*25c28e83SPiotr Jasiukajtis #define P15 cr[5] 99*25c28e83SPiotr Jasiukajtis #define P20 cr[6] 100*25c28e83SPiotr Jasiukajtis #define P21 cr[7] 101*25c28e83SPiotr Jasiukajtis #define P22 cr[8] 102*25c28e83SPiotr Jasiukajtis #define P23 cr[9] 103*25c28e83SPiotr Jasiukajtis #define P24 cr[10] 104*25c28e83SPiotr Jasiukajtis #define P25 cr[11] 105*25c28e83SPiotr Jasiukajtis #define P30 cr[12] 106*25c28e83SPiotr Jasiukajtis #define P31 cr[13] 107*25c28e83SPiotr Jasiukajtis #define P32 cr[14] 108*25c28e83SPiotr Jasiukajtis #define P33 cr[15] 109*25c28e83SPiotr Jasiukajtis #define P34 cr[16] 110*25c28e83SPiotr Jasiukajtis #define GZ1 cr[17] 111*25c28e83SPiotr Jasiukajtis #define GZ2 cr[18] 112*25c28e83SPiotr Jasiukajtis #define GZ3 cr[19] 113*25c28e83SPiotr Jasiukajtis #define TZ1 cr[20] 114*25c28e83SPiotr Jasiukajtis #define TZ3 cr[21] 115*25c28e83SPiotr Jasiukajtis 116*25c28e83SPiotr Jasiukajtis /* compute gamma(y) for y in GT1 = [1.0000, 1.2845] */ 117*25c28e83SPiotr Jasiukajtis static double 118*25c28e83SPiotr Jasiukajtis GT1(double y) { 119*25c28e83SPiotr Jasiukajtis double z, r; 120*25c28e83SPiotr Jasiukajtis 121*25c28e83SPiotr Jasiukajtis z = y * y; 122*25c28e83SPiotr Jasiukajtis r = TZ1 * y + z * ((P10 + y * P11 + z * P12) + (z * y) * (P13 + y * 123*25c28e83SPiotr Jasiukajtis P14 + z * P15)); 124*25c28e83SPiotr Jasiukajtis return (GZ1 + r); 125*25c28e83SPiotr Jasiukajtis } 126*25c28e83SPiotr Jasiukajtis 127*25c28e83SPiotr Jasiukajtis /* compute gamma(y) for y in GT2 = [1.2844, 1.6374] */ 128*25c28e83SPiotr Jasiukajtis static double 129*25c28e83SPiotr Jasiukajtis GT2(double y) { 130*25c28e83SPiotr Jasiukajtis double z; 131*25c28e83SPiotr Jasiukajtis 132*25c28e83SPiotr Jasiukajtis z = y * y; 133*25c28e83SPiotr Jasiukajtis return (GZ2 + z * ((P20 + y * P21 + z * P22) + (z * y) * (P23 + y * 134*25c28e83SPiotr Jasiukajtis P24 + z * P25))); 135*25c28e83SPiotr Jasiukajtis } 136*25c28e83SPiotr Jasiukajtis 137*25c28e83SPiotr Jasiukajtis /* compute gamma(y) for y in GT3 = [1.6373, 2.0000] */ 138*25c28e83SPiotr Jasiukajtis static double 139*25c28e83SPiotr Jasiukajtis GT3(double y) { 140*25c28e83SPiotr Jasiukajtis double z, r; 141*25c28e83SPiotr Jasiukajtis 142*25c28e83SPiotr Jasiukajtis z = y * y; 143*25c28e83SPiotr Jasiukajtis r = TZ3 * y + z * ((P30 + y * P31 + z * P32) + (z * y) * (P33 + y * 144*25c28e83SPiotr Jasiukajtis P34)); 145*25c28e83SPiotr Jasiukajtis return (GZ3 + r); 146*25c28e83SPiotr Jasiukajtis } 147*25c28e83SPiotr Jasiukajtis 148*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 149*25c28e83SPiotr Jasiukajtis static const double c[] = { 150*25c28e83SPiotr Jasiukajtis +1.0, 151*25c28e83SPiotr Jasiukajtis +2.0, 152*25c28e83SPiotr Jasiukajtis +0.5, 153*25c28e83SPiotr Jasiukajtis +1.0e-300, 154*25c28e83SPiotr Jasiukajtis +6.666717231848518054693623697539230e-0001, /* A1=T3[0] */ 155*25c28e83SPiotr Jasiukajtis +8.33333330959694065245736888749042811909994573178e-0002, /* GP[0] */ 156*25c28e83SPiotr Jasiukajtis -2.77765545601667179767706600890361535225507762168e-0003, /* GP[1] */ 157*25c28e83SPiotr Jasiukajtis +7.77830853479775281781085278324621033523037489883e-0004, /* GP[2] */ 158*25c28e83SPiotr Jasiukajtis +4.18938533204672741744150788368695779923320328369e-0001, /* hln2pi */ 159*25c28e83SPiotr Jasiukajtis +2.16608493924982901946e-02, /* ln2_32 */ 160*25c28e83SPiotr Jasiukajtis +4.61662413084468283841e+01, /* invln2_32 */ 161*25c28e83SPiotr Jasiukajtis +5.00004103388988968841156421415669985414073453720e-0001, /* Et1 */ 162*25c28e83SPiotr Jasiukajtis +1.66667656752800761782778277828110208108687545908e-0001, /* Et2 */ 163*25c28e83SPiotr Jasiukajtis }; 164*25c28e83SPiotr Jasiukajtis 165*25c28e83SPiotr Jasiukajtis #define one c[0] 166*25c28e83SPiotr Jasiukajtis #define two c[1] 167*25c28e83SPiotr Jasiukajtis #define half c[2] 168*25c28e83SPiotr Jasiukajtis #define tiny c[3] 169*25c28e83SPiotr Jasiukajtis #define A1 c[4] 170*25c28e83SPiotr Jasiukajtis #define GP0 c[5] 171*25c28e83SPiotr Jasiukajtis #define GP1 c[6] 172*25c28e83SPiotr Jasiukajtis #define GP2 c[7] 173*25c28e83SPiotr Jasiukajtis #define hln2pi c[8] 174*25c28e83SPiotr Jasiukajtis #define ln2_32 c[9] 175*25c28e83SPiotr Jasiukajtis #define invln2_32 c[10] 176*25c28e83SPiotr Jasiukajtis #define Et1 c[11] 177*25c28e83SPiotr Jasiukajtis #define Et2 c[12] 178*25c28e83SPiotr Jasiukajtis 179*25c28e83SPiotr Jasiukajtis /* S[j] = 2**(j/32.) for the final computation of exp(w) */ 180*25c28e83SPiotr Jasiukajtis static const double S[] = { 181*25c28e83SPiotr Jasiukajtis +1.00000000000000000000e+00, /* 3FF0000000000000 */ 182*25c28e83SPiotr Jasiukajtis +1.02189714865411662714e+00, /* 3FF059B0D3158574 */ 183*25c28e83SPiotr Jasiukajtis +1.04427378242741375480e+00, /* 3FF0B5586CF9890F */ 184*25c28e83SPiotr Jasiukajtis +1.06714040067682369717e+00, /* 3FF11301D0125B51 */ 185*25c28e83SPiotr Jasiukajtis +1.09050773266525768967e+00, /* 3FF172B83C7D517B */ 186*25c28e83SPiotr Jasiukajtis +1.11438674259589243221e+00, /* 3FF1D4873168B9AA */ 187*25c28e83SPiotr Jasiukajtis +1.13878863475669156458e+00, /* 3FF2387A6E756238 */ 188*25c28e83SPiotr Jasiukajtis +1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */ 189*25c28e83SPiotr Jasiukajtis +1.18920711500272102690e+00, /* 3FF306FE0A31B715 */ 190*25c28e83SPiotr Jasiukajtis +1.21524735998046895524e+00, /* 3FF371A7373AA9CB */ 191*25c28e83SPiotr Jasiukajtis +1.24185781207348400201e+00, /* 3FF3DEA64C123422 */ 192*25c28e83SPiotr Jasiukajtis +1.26905095719173321989e+00, /* 3FF44E086061892D */ 193*25c28e83SPiotr Jasiukajtis +1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */ 194*25c28e83SPiotr Jasiukajtis +1.32523664315974132322e+00, /* 3FF5342B569D4F82 */ 195*25c28e83SPiotr Jasiukajtis +1.35425554693689265129e+00, /* 3FF5AB07DD485429 */ 196*25c28e83SPiotr Jasiukajtis +1.38390988196383202258e+00, /* 3FF6247EB03A5585 */ 197*25c28e83SPiotr Jasiukajtis +1.41421356237309514547e+00, /* 3FF6A09E667F3BCD */ 198*25c28e83SPiotr Jasiukajtis +1.44518080697704665027e+00, /* 3FF71F75E8EC5F74 */ 199*25c28e83SPiotr Jasiukajtis +1.47682614593949934623e+00, /* 3FF7A11473EB0187 */ 200*25c28e83SPiotr Jasiukajtis +1.50916442759342284141e+00, /* 3FF82589994CCE13 */ 201*25c28e83SPiotr Jasiukajtis +1.54221082540794074411e+00, /* 3FF8ACE5422AA0DB */ 202*25c28e83SPiotr Jasiukajtis +1.57598084510788649659e+00, /* 3FF93737B0CDC5E5 */ 203*25c28e83SPiotr Jasiukajtis +1.61049033194925428347e+00, /* 3FF9C49182A3F090 */ 204*25c28e83SPiotr Jasiukajtis +1.64575547815396494578e+00, /* 3FFA5503B23E255D */ 205*25c28e83SPiotr Jasiukajtis +1.68179283050742900407e+00, /* 3FFAE89F995AD3AD */ 206*25c28e83SPiotr Jasiukajtis +1.71861929812247793414e+00, /* 3FFB7F76F2FB5E47 */ 207*25c28e83SPiotr Jasiukajtis +1.75625216037329945351e+00, /* 3FFC199BDD85529C */ 208*25c28e83SPiotr Jasiukajtis +1.79470907500310716820e+00, /* 3FFCB720DCEF9069 */ 209*25c28e83SPiotr Jasiukajtis +1.83400808640934243066e+00, /* 3FFD5818DCFBA487 */ 210*25c28e83SPiotr Jasiukajtis +1.87416763411029996256e+00, /* 3FFDFC97337B9B5F */ 211*25c28e83SPiotr Jasiukajtis +1.91520656139714740007e+00, /* 3FFEA4AFA2A490DA */ 212*25c28e83SPiotr Jasiukajtis +1.95714412417540017941e+00, /* 3FFF50765B6E4540 */ 213*25c28e83SPiotr Jasiukajtis }; 214*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 215*25c28e83SPiotr Jasiukajtis 216*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 217*25c28e83SPiotr Jasiukajtis /* 218*25c28e83SPiotr Jasiukajtis * return tgammaf(x) in double for 8<x<=35.040096283... using Stirling's formula 219*25c28e83SPiotr Jasiukajtis * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x)) 220*25c28e83SPiotr Jasiukajtis */ 221*25c28e83SPiotr Jasiukajtis /* 222*25c28e83SPiotr Jasiukajtis * compute ss = log(x)-1 223*25c28e83SPiotr Jasiukajtis * 224*25c28e83SPiotr Jasiukajtis * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2, 225*25c28e83SPiotr Jasiukajtis * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and 226*25c28e83SPiotr Jasiukajtis * T1(n-3) = n*log(2)-1, n=3,4,5 227*25c28e83SPiotr Jasiukajtis * T2(j) = log(z[j]), 228*25c28e83SPiotr Jasiukajtis * T3(s) = 2s + A1*s^3 229*25c28e83SPiotr Jasiukajtis * Note 230*25c28e83SPiotr Jasiukajtis * (1) Remez error for T3(s) is bounded by 2**(-35.8) 231*25c28e83SPiotr Jasiukajtis * (see mpremez/work/Log/tgamma_log_2_outr1) 232*25c28e83SPiotr Jasiukajtis */ 233*25c28e83SPiotr Jasiukajtis 234*25c28e83SPiotr Jasiukajtis static const double T1[] = { /* T1[j]=(j+3)*log(2)-1 */ 235*25c28e83SPiotr Jasiukajtis +1.079441541679835928251696364375e+00, 236*25c28e83SPiotr Jasiukajtis +1.772588722239781237668928485833e+00, 237*25c28e83SPiotr Jasiukajtis +2.465735902799726547086160607291e+00, 238*25c28e83SPiotr Jasiukajtis }; 239*25c28e83SPiotr Jasiukajtis 240*25c28e83SPiotr Jasiukajtis static const double T2[] = { /* T2[j]=log(1+j/64+1/128) */ 241*25c28e83SPiotr Jasiukajtis +7.782140442054948947462900061137e-03, 242*25c28e83SPiotr Jasiukajtis +2.316705928153437822879916096229e-02, 243*25c28e83SPiotr Jasiukajtis +3.831886430213659919375532512380e-02, 244*25c28e83SPiotr Jasiukajtis +5.324451451881228286587019378653e-02, 245*25c28e83SPiotr Jasiukajtis +6.795066190850774939456527777263e-02, 246*25c28e83SPiotr Jasiukajtis +8.244366921107459126816006866831e-02, 247*25c28e83SPiotr Jasiukajtis +9.672962645855111229557105648746e-02, 248*25c28e83SPiotr Jasiukajtis +1.108143663402901141948061693232e-01, 249*25c28e83SPiotr Jasiukajtis +1.247034785009572358634065153809e-01, 250*25c28e83SPiotr Jasiukajtis +1.384023228591191356853258736016e-01, 251*25c28e83SPiotr Jasiukajtis +1.519160420258419750718034248969e-01, 252*25c28e83SPiotr Jasiukajtis +1.652495728953071628756114492772e-01, 253*25c28e83SPiotr Jasiukajtis +1.784076574728182971194002415109e-01, 254*25c28e83SPiotr Jasiukajtis +1.913948529996294546092988075613e-01, 255*25c28e83SPiotr Jasiukajtis +2.042155414286908915038203861962e-01, 256*25c28e83SPiotr Jasiukajtis +2.168739383006143596190895257443e-01, 257*25c28e83SPiotr Jasiukajtis +2.293741010648458299914807250461e-01, 258*25c28e83SPiotr Jasiukajtis +2.417199368871451681443075159135e-01, 259*25c28e83SPiotr Jasiukajtis +2.539152099809634441373232979066e-01, 260*25c28e83SPiotr Jasiukajtis +2.659635484971379413391259265375e-01, 261*25c28e83SPiotr Jasiukajtis +2.778684510034563061863500329234e-01, 262*25c28e83SPiotr Jasiukajtis +2.896332925830426768788930555257e-01, 263*25c28e83SPiotr Jasiukajtis +3.012613305781617810128755382338e-01, 264*25c28e83SPiotr Jasiukajtis +3.127557100038968883862465596883e-01, 265*25c28e83SPiotr Jasiukajtis +3.241194686542119760906707604350e-01, 266*25c28e83SPiotr Jasiukajtis +3.353555419211378302571795798142e-01, 267*25c28e83SPiotr Jasiukajtis +3.464667673462085809184621884258e-01, 268*25c28e83SPiotr Jasiukajtis +3.574558889218037742260094901409e-01, 269*25c28e83SPiotr Jasiukajtis +3.683255611587076530482301540504e-01, 270*25c28e83SPiotr Jasiukajtis +3.790783529349694583908533456310e-01, 271*25c28e83SPiotr Jasiukajtis +3.897167511400252133704636040035e-01, 272*25c28e83SPiotr Jasiukajtis +4.002431641270127069293251019951e-01, 273*25c28e83SPiotr Jasiukajtis +4.106599249852683859343062031758e-01, 274*25c28e83SPiotr Jasiukajtis +4.209692946441296361288671615068e-01, 275*25c28e83SPiotr Jasiukajtis +4.311734648183713408591724789556e-01, 276*25c28e83SPiotr Jasiukajtis +4.412745608048752294894964416613e-01, 277*25c28e83SPiotr Jasiukajtis +4.512746441394585851446923830790e-01, 278*25c28e83SPiotr Jasiukajtis +4.611757151221701663679999255979e-01, 279*25c28e83SPiotr Jasiukajtis +4.709797152187910125468978560564e-01, 280*25c28e83SPiotr Jasiukajtis +4.806885293457519076766184554480e-01, 281*25c28e83SPiotr Jasiukajtis +4.903039880451938381503461596457e-01, 282*25c28e83SPiotr Jasiukajtis +4.998278695564493298213314152470e-01, 283*25c28e83SPiotr Jasiukajtis +5.092619017898079468040749192283e-01, 284*25c28e83SPiotr Jasiukajtis +5.186077642080456321529769963648e-01, 285*25c28e83SPiotr Jasiukajtis +5.278670896208423851138922177783e-01, 286*25c28e83SPiotr Jasiukajtis +5.370414658968836545667292441538e-01, 287*25c28e83SPiotr Jasiukajtis +5.461324375981356503823972092312e-01, 288*25c28e83SPiotr Jasiukajtis +5.551415075405015927154803595159e-01, 289*25c28e83SPiotr Jasiukajtis +5.640701382848029660713842900902e-01, 290*25c28e83SPiotr Jasiukajtis +5.729197535617855090927567266263e-01, 291*25c28e83SPiotr Jasiukajtis +5.816917396346224825206107537254e-01, 292*25c28e83SPiotr Jasiukajtis +5.903874466021763746419167081236e-01, 293*25c28e83SPiotr Jasiukajtis +5.990081896460833993816000244617e-01, 294*25c28e83SPiotr Jasiukajtis +6.075552502245417955010851527911e-01, 295*25c28e83SPiotr Jasiukajtis +6.160298772155140196475659281967e-01, 296*25c28e83SPiotr Jasiukajtis +6.244332880118935010425387440547e-01, 297*25c28e83SPiotr Jasiukajtis +6.327666695710378295457864685036e-01, 298*25c28e83SPiotr Jasiukajtis +6.410311794209312910556013344054e-01, 299*25c28e83SPiotr Jasiukajtis +6.492279466251098188908399699053e-01, 300*25c28e83SPiotr Jasiukajtis +6.573580727083600301418900232459e-01, 301*25c28e83SPiotr Jasiukajtis +6.654226325450904489500926100067e-01, 302*25c28e83SPiotr Jasiukajtis +6.734226752121667202979603888010e-01, 303*25c28e83SPiotr Jasiukajtis +6.813592248079030689480715595681e-01, 304*25c28e83SPiotr Jasiukajtis +6.892332812388089803249143378146e-01, 305*25c28e83SPiotr Jasiukajtis }; 306*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 307*25c28e83SPiotr Jasiukajtis 308*25c28e83SPiotr Jasiukajtis static double 309*25c28e83SPiotr Jasiukajtis large_gam(double x) { 310*25c28e83SPiotr Jasiukajtis double ss, zz, z, t1, t2, w, y, u; 311*25c28e83SPiotr Jasiukajtis unsigned lx; 312*25c28e83SPiotr Jasiukajtis int k, ix, j, m; 313*25c28e83SPiotr Jasiukajtis 314*25c28e83SPiotr Jasiukajtis ix = __HI(x); 315*25c28e83SPiotr Jasiukajtis lx = __LO(x); 316*25c28e83SPiotr Jasiukajtis m = (ix >> 20) - 0x3ff; /* exponent of x, range:3-5 */ 317*25c28e83SPiotr Jasiukajtis ix = (ix & 0x000fffff) | 0x3ff00000; /* y = scale x to [1,2] */ 318*25c28e83SPiotr Jasiukajtis __HI(y) = ix; 319*25c28e83SPiotr Jasiukajtis __LO(y) = lx; 320*25c28e83SPiotr Jasiukajtis __HI(z) = (ix & 0xffffc000) | 0x2000; /* z[j]=1+j/64+1/128 */ 321*25c28e83SPiotr Jasiukajtis __LO(z) = 0; 322*25c28e83SPiotr Jasiukajtis j = (ix >> 14) & 0x3f; 323*25c28e83SPiotr Jasiukajtis t1 = y + z; 324*25c28e83SPiotr Jasiukajtis t2 = y - z; 325*25c28e83SPiotr Jasiukajtis u = t2 / t1; 326*25c28e83SPiotr Jasiukajtis ss = T1[m - 3] + T2[j] + u * (two + A1 * (u * u)); 327*25c28e83SPiotr Jasiukajtis /* ss = log(x)-1 */ 328*25c28e83SPiotr Jasiukajtis /* 329*25c28e83SPiotr Jasiukajtis * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2))) 330*25c28e83SPiotr Jasiukajtis * where ss = log(x) - 1 331*25c28e83SPiotr Jasiukajtis */ 332*25c28e83SPiotr Jasiukajtis z = one / x; 333*25c28e83SPiotr Jasiukajtis zz = z * z; 334*25c28e83SPiotr Jasiukajtis w = ((x - half) * ss + hln2pi) + z * (GP0 + zz * GP1 + (zz * zz) * GP2); 335*25c28e83SPiotr Jasiukajtis k = (int) (w * invln2_32 + half); 336*25c28e83SPiotr Jasiukajtis 337*25c28e83SPiotr Jasiukajtis /* compute the exponential of w */ 338*25c28e83SPiotr Jasiukajtis j = k & 0x1f; 339*25c28e83SPiotr Jasiukajtis m = k >> 5; 340*25c28e83SPiotr Jasiukajtis z = w - (double) k *ln2_32; 341*25c28e83SPiotr Jasiukajtis zz = S[j] * (one + z + (z * z) * (Et1 + z * Et2)); 342*25c28e83SPiotr Jasiukajtis __HI(zz) += m << 20; 343*25c28e83SPiotr Jasiukajtis return (zz); 344*25c28e83SPiotr Jasiukajtis } 345*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 346*25c28e83SPiotr Jasiukajtis /* 347*25c28e83SPiotr Jasiukajtis * kpsin(x)= sin(pi*x)/pi 348*25c28e83SPiotr Jasiukajtis * 3 5 7 9 349*25c28e83SPiotr Jasiukajtis * = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x 350*25c28e83SPiotr Jasiukajtis */ 351*25c28e83SPiotr Jasiukajtis static const double ks[] = { 352*25c28e83SPiotr Jasiukajtis -1.64493404985645811354476665052005342839447790544e+0000, 353*25c28e83SPiotr Jasiukajtis +8.11740794458351064092797249069438269367389272270e-0001, 354*25c28e83SPiotr Jasiukajtis -1.90703144603551216933075809162889536878854055202e-0001, 355*25c28e83SPiotr Jasiukajtis +2.55742333994264563281155312271481108635575331201e-0002, 356*25c28e83SPiotr Jasiukajtis }; 357*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 358*25c28e83SPiotr Jasiukajtis 359*25c28e83SPiotr Jasiukajtis static double 360*25c28e83SPiotr Jasiukajtis kpsin(double x) { 361*25c28e83SPiotr Jasiukajtis double z; 362*25c28e83SPiotr Jasiukajtis 363*25c28e83SPiotr Jasiukajtis z = x * x; 364*25c28e83SPiotr Jasiukajtis return (x + (x * z) * ((ks[0] + z * ks[1]) + (z * z) * (ks[2] + z * 365*25c28e83SPiotr Jasiukajtis ks[3]))); 366*25c28e83SPiotr Jasiukajtis } 367*25c28e83SPiotr Jasiukajtis 368*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 369*25c28e83SPiotr Jasiukajtis /* 370*25c28e83SPiotr Jasiukajtis * kpcos(x)= cos(pi*x)/pi 371*25c28e83SPiotr Jasiukajtis * 2 4 6 372*25c28e83SPiotr Jasiukajtis * = kc[0]+kc[1]*x +kc[2]*x +kc[3]*x 373*25c28e83SPiotr Jasiukajtis */ 374*25c28e83SPiotr Jasiukajtis static const double kc[] = { 375*25c28e83SPiotr Jasiukajtis +3.18309886183790671537767526745028724068919291480e-0001, 376*25c28e83SPiotr Jasiukajtis -1.57079581447762568199467875065854538626594937791e+0000, 377*25c28e83SPiotr Jasiukajtis +1.29183528092558692844073004029568674027807393862e+0000, 378*25c28e83SPiotr Jasiukajtis -4.20232949771307685981015914425195471602739075537e-0001, 379*25c28e83SPiotr Jasiukajtis }; 380*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 381*25c28e83SPiotr Jasiukajtis 382*25c28e83SPiotr Jasiukajtis static double 383*25c28e83SPiotr Jasiukajtis kpcos(double x) { 384*25c28e83SPiotr Jasiukajtis double z; 385*25c28e83SPiotr Jasiukajtis 386*25c28e83SPiotr Jasiukajtis z = x * x; 387*25c28e83SPiotr Jasiukajtis return (kc[0] + z * (kc[1] + z * kc[2] + (z * z) * kc[3])); 388*25c28e83SPiotr Jasiukajtis } 389*25c28e83SPiotr Jasiukajtis 390*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 391*25c28e83SPiotr Jasiukajtis static const double 392*25c28e83SPiotr Jasiukajtis t0z1 = 0.134861805732790769689793935774652917006, 393*25c28e83SPiotr Jasiukajtis t0z2 = 0.461632144968362341262659542325721328468, 394*25c28e83SPiotr Jasiukajtis t0z3 = 0.819773101100500601787868704921606996312; 395*25c28e83SPiotr Jasiukajtis /* 1.134861805732790769689793935774652917006 */ 396*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 397*25c28e83SPiotr Jasiukajtis 398*25c28e83SPiotr Jasiukajtis /* 399*25c28e83SPiotr Jasiukajtis * gamma(x+i) for 0 <= x < 1 400*25c28e83SPiotr Jasiukajtis */ 401*25c28e83SPiotr Jasiukajtis static double 402*25c28e83SPiotr Jasiukajtis gam_n(int i, double x) { 403*25c28e83SPiotr Jasiukajtis double rr = 0.0L, yy; 404*25c28e83SPiotr Jasiukajtis double z1, z2; 405*25c28e83SPiotr Jasiukajtis 406*25c28e83SPiotr Jasiukajtis /* compute yy = gamma(x+1) */ 407*25c28e83SPiotr Jasiukajtis if (x > 0.2845) { 408*25c28e83SPiotr Jasiukajtis if (x > 0.6374) 409*25c28e83SPiotr Jasiukajtis yy = GT3(x - t0z3); 410*25c28e83SPiotr Jasiukajtis else 411*25c28e83SPiotr Jasiukajtis yy = GT2(x - t0z2); 412*25c28e83SPiotr Jasiukajtis } else 413*25c28e83SPiotr Jasiukajtis yy = GT1(x - t0z1); 414*25c28e83SPiotr Jasiukajtis 415*25c28e83SPiotr Jasiukajtis /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */ 416*25c28e83SPiotr Jasiukajtis switch (i) { 417*25c28e83SPiotr Jasiukajtis case 0: /* yy/x */ 418*25c28e83SPiotr Jasiukajtis rr = yy / x; 419*25c28e83SPiotr Jasiukajtis break; 420*25c28e83SPiotr Jasiukajtis case 1: /* yy */ 421*25c28e83SPiotr Jasiukajtis rr = yy; 422*25c28e83SPiotr Jasiukajtis break; 423*25c28e83SPiotr Jasiukajtis case 2: /* (x+1)*yy */ 424*25c28e83SPiotr Jasiukajtis rr = (x + one) * yy; 425*25c28e83SPiotr Jasiukajtis break; 426*25c28e83SPiotr Jasiukajtis case 3: /* (x+2)*(x+1)*yy */ 427*25c28e83SPiotr Jasiukajtis rr = (x + one) * (x + two) * yy; 428*25c28e83SPiotr Jasiukajtis break; 429*25c28e83SPiotr Jasiukajtis 430*25c28e83SPiotr Jasiukajtis case 4: /* (x+1)*(x+3)*(x+2)*yy */ 431*25c28e83SPiotr Jasiukajtis rr = (x + one) * (x + two) * ((x + 3.0) * yy); 432*25c28e83SPiotr Jasiukajtis break; 433*25c28e83SPiotr Jasiukajtis case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */ 434*25c28e83SPiotr Jasiukajtis z1 = (x + two) * (x + 3.0) * yy; 435*25c28e83SPiotr Jasiukajtis z2 = (x + one) * (x + 4.0); 436*25c28e83SPiotr Jasiukajtis rr = z1 * z2; 437*25c28e83SPiotr Jasiukajtis break; 438*25c28e83SPiotr Jasiukajtis case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */ 439*25c28e83SPiotr Jasiukajtis z1 = (x + two) * (x + 3.0); 440*25c28e83SPiotr Jasiukajtis z2 = (x + 5.0) * yy; 441*25c28e83SPiotr Jasiukajtis rr = z1 * (z1 - two) * z2; 442*25c28e83SPiotr Jasiukajtis break; 443*25c28e83SPiotr Jasiukajtis case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */ 444*25c28e83SPiotr Jasiukajtis z1 = (x + two) * (x + 3.0); 445*25c28e83SPiotr Jasiukajtis z2 = (x + 5.0) * (x + 6.0) * yy; 446*25c28e83SPiotr Jasiukajtis rr = z1 * (z1 - two) * z2; 447*25c28e83SPiotr Jasiukajtis break; 448*25c28e83SPiotr Jasiukajtis } 449*25c28e83SPiotr Jasiukajtis return (rr); 450*25c28e83SPiotr Jasiukajtis } 451*25c28e83SPiotr Jasiukajtis 452*25c28e83SPiotr Jasiukajtis float 453*25c28e83SPiotr Jasiukajtis tgammaf(float xf) { 454*25c28e83SPiotr Jasiukajtis float zf; 455*25c28e83SPiotr Jasiukajtis double ss, ww; 456*25c28e83SPiotr Jasiukajtis double x, y, z; 457*25c28e83SPiotr Jasiukajtis int i, j, k, ix, hx, xk; 458*25c28e83SPiotr Jasiukajtis 459*25c28e83SPiotr Jasiukajtis hx = *(int *) &xf; 460*25c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff; 461*25c28e83SPiotr Jasiukajtis 462*25c28e83SPiotr Jasiukajtis x = (double) xf; 463*25c28e83SPiotr Jasiukajtis if (ix < 0x33800000) 464*25c28e83SPiotr Jasiukajtis return (1.0F / xf); /* |x| < 2**-24 */ 465*25c28e83SPiotr Jasiukajtis 466*25c28e83SPiotr Jasiukajtis if (ix >= 0x7f800000) 467*25c28e83SPiotr Jasiukajtis return (xf * ((hx < 0)? 0.0F : xf)); /* +-Inf or NaN */ 468*25c28e83SPiotr Jasiukajtis 469*25c28e83SPiotr Jasiukajtis if (hx > 0x420C290F) /* x > 35.040096283... overflow */ 470*25c28e83SPiotr Jasiukajtis return (float)(x / tiny); 471*25c28e83SPiotr Jasiukajtis 472*25c28e83SPiotr Jasiukajtis if (hx >= 0x41000000) /* x >= 8 */ 473*25c28e83SPiotr Jasiukajtis return ((float) large_gam(x)); 474*25c28e83SPiotr Jasiukajtis 475*25c28e83SPiotr Jasiukajtis if (hx > 0) { /* 0 < x < 8 */ 476*25c28e83SPiotr Jasiukajtis i = (int) xf; 477*25c28e83SPiotr Jasiukajtis return ((float) gam_n(i, x - (double) i)); 478*25c28e83SPiotr Jasiukajtis } 479*25c28e83SPiotr Jasiukajtis 480*25c28e83SPiotr Jasiukajtis /* negative x */ 481*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 482*25c28e83SPiotr Jasiukajtis /* 483*25c28e83SPiotr Jasiukajtis * compute xk = 484*25c28e83SPiotr Jasiukajtis * -2 ... x is an even int (-inf is considered even) 485*25c28e83SPiotr Jasiukajtis * -1 ... x is an odd int 486*25c28e83SPiotr Jasiukajtis * +0 ... x is not an int but chopped to an even int 487*25c28e83SPiotr Jasiukajtis * +1 ... x is not an int but chopped to an odd int 488*25c28e83SPiotr Jasiukajtis */ 489*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 490*25c28e83SPiotr Jasiukajtis xk = 0; 491*25c28e83SPiotr Jasiukajtis if (ix >= 0x4b000000) { 492*25c28e83SPiotr Jasiukajtis if (ix > 0x4b000000) 493*25c28e83SPiotr Jasiukajtis xk = -2; 494*25c28e83SPiotr Jasiukajtis else 495*25c28e83SPiotr Jasiukajtis xk = -2 + (ix & 1); 496*25c28e83SPiotr Jasiukajtis } else if (ix >= 0x3f800000) { 497*25c28e83SPiotr Jasiukajtis k = (ix >> 23) - 0x7f; 498*25c28e83SPiotr Jasiukajtis j = ix >> (23 - k); 499*25c28e83SPiotr Jasiukajtis if ((j << (23 - k)) == ix) 500*25c28e83SPiotr Jasiukajtis xk = -2 + (j & 1); 501*25c28e83SPiotr Jasiukajtis else 502*25c28e83SPiotr Jasiukajtis xk = j & 1; 503*25c28e83SPiotr Jasiukajtis } 504*25c28e83SPiotr Jasiukajtis if (xk < 0) { 505*25c28e83SPiotr Jasiukajtis /* 0/0 invalid NaN, ideally gamma(-n)= (-1)**(n+1) * inf */ 506*25c28e83SPiotr Jasiukajtis zf = xf - xf; 507*25c28e83SPiotr Jasiukajtis return (zf / zf); 508*25c28e83SPiotr Jasiukajtis } 509*25c28e83SPiotr Jasiukajtis 510*25c28e83SPiotr Jasiukajtis /* negative underflow thresold */ 511*25c28e83SPiotr Jasiukajtis if (ix > 0x4224000B) { /* x < -(41+11ulp) */ 512*25c28e83SPiotr Jasiukajtis if (xk == 0) 513*25c28e83SPiotr Jasiukajtis z = -tiny; 514*25c28e83SPiotr Jasiukajtis else 515*25c28e83SPiotr Jasiukajtis z = tiny; 516*25c28e83SPiotr Jasiukajtis return ((float)z); 517*25c28e83SPiotr Jasiukajtis } 518*25c28e83SPiotr Jasiukajtis 519*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 520*25c28e83SPiotr Jasiukajtis /* now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */ 521*25c28e83SPiotr Jasiukajtis /* 522*25c28e83SPiotr Jasiukajtis * First compute ss = -sin(pi*y)/pi , so that 523*25c28e83SPiotr Jasiukajtis * gamma(x) = 1/(ss*gamma(1+y)) 524*25c28e83SPiotr Jasiukajtis */ 525*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 526*25c28e83SPiotr Jasiukajtis y = -x; 527*25c28e83SPiotr Jasiukajtis j = (int) y; 528*25c28e83SPiotr Jasiukajtis z = y - (double) j; 529*25c28e83SPiotr Jasiukajtis if (z > 0.3183098861837906715377675) 530*25c28e83SPiotr Jasiukajtis if (z > 0.6816901138162093284622325) 531*25c28e83SPiotr Jasiukajtis ss = kpsin(one - z); 532*25c28e83SPiotr Jasiukajtis else 533*25c28e83SPiotr Jasiukajtis ss = kpcos(0.5 - z); 534*25c28e83SPiotr Jasiukajtis else 535*25c28e83SPiotr Jasiukajtis ss = kpsin(z); 536*25c28e83SPiotr Jasiukajtis if (xk == 0) 537*25c28e83SPiotr Jasiukajtis ss = -ss; 538*25c28e83SPiotr Jasiukajtis 539*25c28e83SPiotr Jasiukajtis /* Then compute ww = gamma(1+y) */ 540*25c28e83SPiotr Jasiukajtis if (j < 7) 541*25c28e83SPiotr Jasiukajtis ww = gam_n(j + 1, z); 542*25c28e83SPiotr Jasiukajtis else 543*25c28e83SPiotr Jasiukajtis ww = large_gam(y + one); 544*25c28e83SPiotr Jasiukajtis 545*25c28e83SPiotr Jasiukajtis /* return 1/(ss*ww) */ 546*25c28e83SPiotr Jasiukajtis return ((float) (one / (ww * ss))); 547*25c28e83SPiotr Jasiukajtis } 548