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 #if defined(ELFOBJ) 31*25c28e83SPiotr Jasiukajtis #pragma weak tgammal = __tgammal 32*25c28e83SPiotr Jasiukajtis #endif 33*25c28e83SPiotr Jasiukajtis 34*25c28e83SPiotr Jasiukajtis #include "libm.h" 35*25c28e83SPiotr Jasiukajtis #include <sys/isa_defs.h> 36*25c28e83SPiotr Jasiukajtis 37*25c28e83SPiotr Jasiukajtis #if defined(_BIG_ENDIAN) 38*25c28e83SPiotr Jasiukajtis #define H0_WORD(x) ((unsigned *) &x)[0] 39*25c28e83SPiotr Jasiukajtis #define H3_WORD(x) ((unsigned *) &x)[3] 40*25c28e83SPiotr Jasiukajtis #define CHOPPED(x) (long double) ((double) (x)) 41*25c28e83SPiotr Jasiukajtis #else 42*25c28e83SPiotr Jasiukajtis #define H0_WORD(x) ((((int *) &x)[2] << 16) | \ 43*25c28e83SPiotr Jasiukajtis (0x0000ffff & (((unsigned *) &x)[1] >> 15))) 44*25c28e83SPiotr Jasiukajtis #define H3_WORD(x) ((unsigned *) &x)[0] 45*25c28e83SPiotr Jasiukajtis #define CHOPPED(x) (long double) ((float) (x)) 46*25c28e83SPiotr Jasiukajtis #endif 47*25c28e83SPiotr Jasiukajtis 48*25c28e83SPiotr Jasiukajtis struct LDouble { 49*25c28e83SPiotr Jasiukajtis long double h, l; 50*25c28e83SPiotr Jasiukajtis }; 51*25c28e83SPiotr Jasiukajtis 52*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 53*25c28e83SPiotr Jasiukajtis /* Primary interval GTi() */ 54*25c28e83SPiotr Jasiukajtis static const long double P1[] = { 55*25c28e83SPiotr Jasiukajtis +0.709086836199777919037185741507610124611513720557L, 56*25c28e83SPiotr Jasiukajtis +4.45754781206489035827915969367354835667391606951e-0001L, 57*25c28e83SPiotr Jasiukajtis +3.21049298735832382311662273882632210062918153852e-0002L, 58*25c28e83SPiotr Jasiukajtis -5.71296796342106617651765245858289197369688864350e-0003L, 59*25c28e83SPiotr Jasiukajtis +6.04666892891998977081619174969855831606965352773e-0003L, 60*25c28e83SPiotr Jasiukajtis +8.99106186996888711939627812174765258822658645168e-0004L, 61*25c28e83SPiotr Jasiukajtis -6.96496846144407741431207008527018441810175568949e-0005L, 62*25c28e83SPiotr Jasiukajtis +1.52597046118984020814225409300131445070213882429e-0005L, 63*25c28e83SPiotr Jasiukajtis +5.68521076168495673844711465407432189190681541547e-0007L, 64*25c28e83SPiotr Jasiukajtis +3.30749673519634895220582062520286565610418952979e-0008L, 65*25c28e83SPiotr Jasiukajtis }; 66*25c28e83SPiotr Jasiukajtis static const long double Q1[] = { 67*25c28e83SPiotr Jasiukajtis +1.0+0000L, 68*25c28e83SPiotr Jasiukajtis +1.35806511721671070408570853537257079579490650668e+0000L, 69*25c28e83SPiotr Jasiukajtis +2.97567810153429553405327140096063086994072952961e-0001L, 70*25c28e83SPiotr Jasiukajtis -1.52956835982588571502954372821681851681118097870e-0001L, 71*25c28e83SPiotr Jasiukajtis -2.88248519561420109768781615289082053597954521218e-0002L, 72*25c28e83SPiotr Jasiukajtis +1.03475311719937405219789948456313936302378395955e-0002L, 73*25c28e83SPiotr Jasiukajtis +4.12310203243891222368965360124391297374822742313e-0004L, 74*25c28e83SPiotr Jasiukajtis -3.12653708152290867248931925120380729518332507388e-0004L, 75*25c28e83SPiotr Jasiukajtis +2.36672170850409745237358105667757760527014332458e-0005L, 76*25c28e83SPiotr Jasiukajtis }; 77*25c28e83SPiotr Jasiukajtis static const long double P2[] = { 78*25c28e83SPiotr Jasiukajtis +0.428486815855585429730209907810650135255270600668084114L, 79*25c28e83SPiotr Jasiukajtis +2.62768479103809762805691743305424077975230551176e-0001L, 80*25c28e83SPiotr Jasiukajtis +3.81187532685392297608310837995193946591425896150e-0002L, 81*25c28e83SPiotr Jasiukajtis +3.00063075891811043820666846129131255948527925381e-0003L, 82*25c28e83SPiotr Jasiukajtis +2.47315407812279164228398470797498649142513408654e-0003L, 83*25c28e83SPiotr Jasiukajtis +3.62838199917848372586173483147214880464782938664e-0004L, 84*25c28e83SPiotr Jasiukajtis +3.43991105975492623982725644046473030098172692423e-0006L, 85*25c28e83SPiotr Jasiukajtis +4.56902151569603272237014240794257659159045432895e-0006L, 86*25c28e83SPiotr Jasiukajtis +2.13734755837595695602045100675540011352948958453e-0007L, 87*25c28e83SPiotr Jasiukajtis +9.74123440547918230781670266967882492234877125358e-0009L, 88*25c28e83SPiotr Jasiukajtis }; 89*25c28e83SPiotr Jasiukajtis static const long double Q2[] = { 90*25c28e83SPiotr Jasiukajtis +1.0L, 91*25c28e83SPiotr Jasiukajtis +9.18284118632506842664645516830761489700556179701e-0001L, 92*25c28e83SPiotr Jasiukajtis -6.41430858837830766045202076965923776189154874947e-0003L, 93*25c28e83SPiotr Jasiukajtis -1.24400885809771073213345747437964149775410921376e-0001L, 94*25c28e83SPiotr Jasiukajtis +4.69803798146251757538856567522481979624746875964e-0003L, 95*25c28e83SPiotr Jasiukajtis +7.18309447069495315914284705109868696262662082731e-0003L, 96*25c28e83SPiotr Jasiukajtis -8.75812626987894695112722600697653425786166399105e-0004L, 97*25c28e83SPiotr Jasiukajtis -1.23539972377769277995959339188431498626674835169e-0004L, 98*25c28e83SPiotr Jasiukajtis +3.10019017590151598732360097849672925448587547746e-0005L, 99*25c28e83SPiotr Jasiukajtis -1.77260223349332617658921874288026777465782364070e-0006L, 100*25c28e83SPiotr Jasiukajtis }; 101*25c28e83SPiotr Jasiukajtis static const long double P3[] = { 102*25c28e83SPiotr Jasiukajtis +0.3824094797345675048502747661075355640070439388902L, 103*25c28e83SPiotr Jasiukajtis +3.42198093076618495415854906335908427159833377774e-0001L, 104*25c28e83SPiotr Jasiukajtis +9.63828189500585568303961406863153237440702754858e-0002L, 105*25c28e83SPiotr Jasiukajtis +8.76069421042696384852462044188520252156846768667e-0003L, 106*25c28e83SPiotr Jasiukajtis +1.86477890389161491224872014149309015261897537488e-0003L, 107*25c28e83SPiotr Jasiukajtis +8.16871354540309895879974742853701311541286944191e-0004L, 108*25c28e83SPiotr Jasiukajtis +6.83783483674600322518695090864659381650125625216e-0005L, 109*25c28e83SPiotr Jasiukajtis -1.10168269719261574708565935172719209272190828456e-0006L, 110*25c28e83SPiotr Jasiukajtis +9.66243228508380420159234853278906717065629721016e-0007L, 111*25c28e83SPiotr Jasiukajtis +2.31858885579177250541163820671121664974334728142e-0008L, 112*25c28e83SPiotr Jasiukajtis }; 113*25c28e83SPiotr Jasiukajtis static const long double Q3[] = { 114*25c28e83SPiotr Jasiukajtis +1.0L, 115*25c28e83SPiotr Jasiukajtis +8.25479821168813634632437430090376252512793067339e-0001L, 116*25c28e83SPiotr Jasiukajtis -1.62251363073937769739639623669295110346015576320e-0002L, 117*25c28e83SPiotr Jasiukajtis -1.10621286905916732758745130629426559691187579852e-0001L, 118*25c28e83SPiotr Jasiukajtis +3.48309693970985612644446415789230015515365291459e-0003L, 119*25c28e83SPiotr Jasiukajtis +6.73553737487488333032431261131289672347043401328e-0003L, 120*25c28e83SPiotr Jasiukajtis -7.63222008393372630162743587811004613050245128051e-0004L, 121*25c28e83SPiotr Jasiukajtis -1.35792670669190631476784768961953711773073251336e-0004L, 122*25c28e83SPiotr Jasiukajtis +3.19610150954223587006220730065608156460205690618e-0005L, 123*25c28e83SPiotr Jasiukajtis -1.82096553862822346610109522015129585693354348322e-0006L, 124*25c28e83SPiotr Jasiukajtis }; 125*25c28e83SPiotr Jasiukajtis 126*25c28e83SPiotr Jasiukajtis static const long double 127*25c28e83SPiotr Jasiukajtis #if defined(__x86) 128*25c28e83SPiotr Jasiukajtis GZ1_h = 0.938204627909682449364570100414084663498215377L, 129*25c28e83SPiotr Jasiukajtis GZ1_l = 4.518346116624229420055327632718530617227944106e-20L, 130*25c28e83SPiotr Jasiukajtis GZ2_h = 0.885603194410888700264725126309883762587560340L, 131*25c28e83SPiotr Jasiukajtis GZ2_l = 1.409077427270497062039119290776508217077297169e-20L, 132*25c28e83SPiotr Jasiukajtis GZ3_h = 0.936781411463652321613537060640553022494714241L, 133*25c28e83SPiotr Jasiukajtis GZ3_l = 5.309836440284827247897772963887219035221996813e-21L, 134*25c28e83SPiotr Jasiukajtis #else 135*25c28e83SPiotr Jasiukajtis GZ1_h = 0.938204627909682449409753561580326910854647031L, 136*25c28e83SPiotr Jasiukajtis GZ1_l = 4.684412162199460089642452580902345976446297037e-35L, 137*25c28e83SPiotr Jasiukajtis GZ2_h = 0.885603194410888700278815900582588658192658794L, 138*25c28e83SPiotr Jasiukajtis GZ2_l = 7.501529273890253789219935569758713534641074860e-35L, 139*25c28e83SPiotr Jasiukajtis GZ3_h = 0.936781411463652321618846897080837818855399840L, 140*25c28e83SPiotr Jasiukajtis GZ3_l = 3.088721217404784363585591914529361687403776917e-35L, 141*25c28e83SPiotr Jasiukajtis #endif 142*25c28e83SPiotr Jasiukajtis TZ1 = -0.3517214357852935791015625L, 143*25c28e83SPiotr Jasiukajtis TZ3 = 0.280530631542205810546875L; 144*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 145*25c28e83SPiotr Jasiukajtis 146*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 147*25c28e83SPiotr Jasiukajtis /* 148*25c28e83SPiotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845] 149*25c28e83SPiotr Jasiukajtis * ...assume yh got 53 or 24(i386) significant bits 150*25c28e83SPiotr Jasiukajtis */ 151*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 152*25c28e83SPiotr Jasiukajtis static struct LDouble 153*25c28e83SPiotr Jasiukajtis GT1(long double yh, long double yl) { 154*25c28e83SPiotr Jasiukajtis long double t3, t4, y; 155*25c28e83SPiotr Jasiukajtis int i; 156*25c28e83SPiotr Jasiukajtis struct LDouble r; 157*25c28e83SPiotr Jasiukajtis 158*25c28e83SPiotr Jasiukajtis y = yh + yl; 159*25c28e83SPiotr Jasiukajtis for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) { 160*25c28e83SPiotr Jasiukajtis t4 = t4 * y + Q1[i]; 161*25c28e83SPiotr Jasiukajtis t3 = t3 * y + P1[i]; 162*25c28e83SPiotr Jasiukajtis } 163*25c28e83SPiotr Jasiukajtis t3 = (y * y) * t3 / t4; 164*25c28e83SPiotr Jasiukajtis t3 += (TZ1 * yl + GZ1_l); 165*25c28e83SPiotr Jasiukajtis t4 = TZ1 * yh; 166*25c28e83SPiotr Jasiukajtis r.h = CHOPPED((t4 + GZ1_h + t3)); 167*25c28e83SPiotr Jasiukajtis t3 += (t4 - (r.h - GZ1_h)); 168*25c28e83SPiotr Jasiukajtis r.l = t3; 169*25c28e83SPiotr Jasiukajtis return (r); 170*25c28e83SPiotr Jasiukajtis } 171*25c28e83SPiotr Jasiukajtis 172*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 173*25c28e83SPiotr Jasiukajtis /* 174*25c28e83SPiotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374] 175*25c28e83SPiotr Jasiukajtis * ...assume yh got 53 significant bits 176*25c28e83SPiotr Jasiukajtis */ 177*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 178*25c28e83SPiotr Jasiukajtis static struct LDouble 179*25c28e83SPiotr Jasiukajtis GT2(long double yh, long double yl) { 180*25c28e83SPiotr Jasiukajtis long double t3, t4, y; 181*25c28e83SPiotr Jasiukajtis int i; 182*25c28e83SPiotr Jasiukajtis struct LDouble r; 183*25c28e83SPiotr Jasiukajtis 184*25c28e83SPiotr Jasiukajtis y = yh + yl; 185*25c28e83SPiotr Jasiukajtis for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) { 186*25c28e83SPiotr Jasiukajtis t4 = t4 * y + Q2[i]; 187*25c28e83SPiotr Jasiukajtis t3 = t3 * y + P2[i]; 188*25c28e83SPiotr Jasiukajtis } 189*25c28e83SPiotr Jasiukajtis t3 = GZ2_l + (y * y) * t3 / t4; 190*25c28e83SPiotr Jasiukajtis r.h = CHOPPED((GZ2_h + t3)); 191*25c28e83SPiotr Jasiukajtis r.l = t3 - (r.h - GZ2_h); 192*25c28e83SPiotr Jasiukajtis return (r); 193*25c28e83SPiotr Jasiukajtis } 194*25c28e83SPiotr Jasiukajtis 195*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 196*25c28e83SPiotr Jasiukajtis /* 197*25c28e83SPiotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000] 198*25c28e83SPiotr Jasiukajtis * ...assume yh got 53 significant bits 199*25c28e83SPiotr Jasiukajtis */ 200*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 201*25c28e83SPiotr Jasiukajtis static struct LDouble 202*25c28e83SPiotr Jasiukajtis GT3(long double yh, long double yl) { 203*25c28e83SPiotr Jasiukajtis long double t3, t4, y; 204*25c28e83SPiotr Jasiukajtis int i; 205*25c28e83SPiotr Jasiukajtis struct LDouble r; 206*25c28e83SPiotr Jasiukajtis 207*25c28e83SPiotr Jasiukajtis y = yh + yl; 208*25c28e83SPiotr Jasiukajtis for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) { 209*25c28e83SPiotr Jasiukajtis t4 = t4 * y + Q3[i]; 210*25c28e83SPiotr Jasiukajtis t3 = t3 * y + P3[i]; 211*25c28e83SPiotr Jasiukajtis } 212*25c28e83SPiotr Jasiukajtis t3 = (y * y) * t3 / t4; 213*25c28e83SPiotr Jasiukajtis t3 += (TZ3 * yl + GZ3_l); 214*25c28e83SPiotr Jasiukajtis t4 = TZ3 * yh; 215*25c28e83SPiotr Jasiukajtis r.h = CHOPPED((t4 + GZ3_h + t3)); 216*25c28e83SPiotr Jasiukajtis t3 += (t4 - (r.h - GZ3_h)); 217*25c28e83SPiotr Jasiukajtis r.l = t3; 218*25c28e83SPiotr Jasiukajtis return (r); 219*25c28e83SPiotr Jasiukajtis } 220*25c28e83SPiotr Jasiukajtis 221*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 222*25c28e83SPiotr Jasiukajtis /* Hex value of GP[0] shoule be 3FB55555 55555555 */ 223*25c28e83SPiotr Jasiukajtis static const long double GP[] = { 224*25c28e83SPiotr Jasiukajtis +0.083333333333333333333333333333333172839171301L, 225*25c28e83SPiotr Jasiukajtis -2.77777777777777777777777777492501211999399424104e-0003L, 226*25c28e83SPiotr Jasiukajtis +7.93650793650793650793635650541638236350020883243e-0004L, 227*25c28e83SPiotr Jasiukajtis -5.95238095238095238057299772679324503339241961704e-0004L, 228*25c28e83SPiotr Jasiukajtis +8.41750841750841696138422987977683524926142600321e-0004L, 229*25c28e83SPiotr Jasiukajtis -1.91752691752686682825032547823699662178842123308e-0003L, 230*25c28e83SPiotr Jasiukajtis +6.41025641022403480921891559356473451161279359322e-0003L, 231*25c28e83SPiotr Jasiukajtis -2.95506535798414019189819587455577003732808185071e-0002L, 232*25c28e83SPiotr Jasiukajtis +1.79644367229970031486079180060923073476568732136e-0001L, 233*25c28e83SPiotr Jasiukajtis -1.39243086487274662174562872567057200255649290646e+0000L, 234*25c28e83SPiotr Jasiukajtis +1.34025874044417962188677816477842265259608269775e+0001L, 235*25c28e83SPiotr Jasiukajtis -1.56803713480127469414495545399982508700748274318e+0002L, 236*25c28e83SPiotr Jasiukajtis +2.18739841656201561694927630335099313968924493891e+0003L, 237*25c28e83SPiotr Jasiukajtis -3.55249848644100338419187038090925410976237921269e+0004L, 238*25c28e83SPiotr Jasiukajtis +6.43464880437835286216768959439484376449179576452e+0005L, 239*25c28e83SPiotr Jasiukajtis -1.20459154385577014992600342782821389605893904624e+0007L, 240*25c28e83SPiotr Jasiukajtis +2.09263249637351298563934942349749718491071093210e+0008L, 241*25c28e83SPiotr Jasiukajtis -2.96247483183169219343745316433899599834685703457e+0009L, 242*25c28e83SPiotr Jasiukajtis +2.88984933605896033154727626086506756972327292981e+0010L, 243*25c28e83SPiotr Jasiukajtis -1.40960434146030007732838382416230610302678063984e+0011L, /* 19 */ 244*25c28e83SPiotr Jasiukajtis }; 245*25c28e83SPiotr Jasiukajtis 246*25c28e83SPiotr Jasiukajtis static const long double T3[] = { 247*25c28e83SPiotr Jasiukajtis +0.666666666666666666666666666666666634567834260213L, /* T3[0] */ 248*25c28e83SPiotr Jasiukajtis +0.400000000000000000000000000040853636176634934140L, /* T3[1] */ 249*25c28e83SPiotr Jasiukajtis +0.285714285714285714285696975252753987869020263448L, /* T3[2] */ 250*25c28e83SPiotr Jasiukajtis +0.222222222222222225593221101192317258554772129875L, /* T3[3] */ 251*25c28e83SPiotr Jasiukajtis +0.181818181817850192105847183461778186703779262916L, /* T3[4] */ 252*25c28e83SPiotr Jasiukajtis +0.153846169861348633757101285952333369222567014596L, /* T3[5] */ 253*25c28e83SPiotr Jasiukajtis +0.133033462889260193922261296772841229985047571265L, /* T3[6] */ 254*25c28e83SPiotr Jasiukajtis }; 255*25c28e83SPiotr Jasiukajtis 256*25c28e83SPiotr Jasiukajtis static const long double c[] = { 257*25c28e83SPiotr Jasiukajtis 0.0L, 258*25c28e83SPiotr Jasiukajtis 1.0L, 259*25c28e83SPiotr Jasiukajtis 2.0L, 260*25c28e83SPiotr Jasiukajtis 0.5L, 261*25c28e83SPiotr Jasiukajtis 1.0e-4930L, /* tiny */ 262*25c28e83SPiotr Jasiukajtis 4.18937683105468750000e-01L, /* hln2pim1_h */ 263*25c28e83SPiotr Jasiukajtis 8.50099203991780329736405617639861397473637783412817152e-07L, /* hln2pim1_l */ 264*25c28e83SPiotr Jasiukajtis 0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */ 265*25c28e83SPiotr Jasiukajtis 2.16608493865351192653179168701171875e-02L, /* ln2_32hi */ 266*25c28e83SPiotr Jasiukajtis 5.96317165397058692545083025235937919875797669127130e-12L, /* ln2_32lo */ 267*25c28e83SPiotr Jasiukajtis 46.16624130844682903551758979206054839765267053289554989233L, /* invln2_32 */ 268*25c28e83SPiotr Jasiukajtis #if defined(__x86) 269*25c28e83SPiotr Jasiukajtis 1.7555483429044629170023839037639845628291e+03L, /* overflow */ 270*25c28e83SPiotr Jasiukajtis #else 271*25c28e83SPiotr Jasiukajtis 1.7555483429044629170038892160702032034177e+03L, /* overflow */ 272*25c28e83SPiotr Jasiukajtis #endif 273*25c28e83SPiotr Jasiukajtis }; 274*25c28e83SPiotr Jasiukajtis 275*25c28e83SPiotr Jasiukajtis #define zero c[0] 276*25c28e83SPiotr Jasiukajtis #define one c[1] 277*25c28e83SPiotr Jasiukajtis #define two c[2] 278*25c28e83SPiotr Jasiukajtis #define half c[3] 279*25c28e83SPiotr Jasiukajtis #define tiny c[4] 280*25c28e83SPiotr Jasiukajtis #define hln2pim1_h c[5] 281*25c28e83SPiotr Jasiukajtis #define hln2pim1_l c[6] 282*25c28e83SPiotr Jasiukajtis #define hln2pim1 c[7] 283*25c28e83SPiotr Jasiukajtis #define ln2_32hi c[8] 284*25c28e83SPiotr Jasiukajtis #define ln2_32lo c[9] 285*25c28e83SPiotr Jasiukajtis #define invln2_32 c[10] 286*25c28e83SPiotr Jasiukajtis #define overflow c[11] 287*25c28e83SPiotr Jasiukajtis 288*25c28e83SPiotr Jasiukajtis /* 289*25c28e83SPiotr Jasiukajtis * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64 290*25c28e83SPiotr Jasiukajtis */ 291*25c28e83SPiotr Jasiukajtis static const long double Et[] = { 292*25c28e83SPiotr Jasiukajtis +5.0000000000000000000e-1L, 293*25c28e83SPiotr Jasiukajtis +1.66666666666666666666666666666828835166292152466e-0001L, 294*25c28e83SPiotr Jasiukajtis +4.16666666666666666666666666666693398646592712189e-0002L, 295*25c28e83SPiotr Jasiukajtis +8.33333333333333333333331748774512601775591115951e-0003L, 296*25c28e83SPiotr Jasiukajtis +1.38888888888888888888888845356011511394764753997e-0003L, 297*25c28e83SPiotr Jasiukajtis +1.98412698412698413237140350092993252684198882102e-0004L, 298*25c28e83SPiotr Jasiukajtis +2.48015873015873016080222025357442659895814371694e-0005L, 299*25c28e83SPiotr Jasiukajtis +2.75573192239028921114572986441972140933432317798e-0006L, 300*25c28e83SPiotr Jasiukajtis +2.75573192239448470555548102895526369739856219317e-0007L, 301*25c28e83SPiotr Jasiukajtis +2.50521677867683935940853997995937600214167232477e-0008L, 302*25c28e83SPiotr Jasiukajtis +2.08767928899010367374984448513685566514152147362e-0009L, 303*25c28e83SPiotr Jasiukajtis }; 304*25c28e83SPiotr Jasiukajtis 305*25c28e83SPiotr Jasiukajtis /* 306*25c28e83SPiotr Jasiukajtis * long double precision coefficients for computing log(x)-1 in tgamma. 307*25c28e83SPiotr Jasiukajtis * See "algorithm" for details 308*25c28e83SPiotr Jasiukajtis * 309*25c28e83SPiotr Jasiukajtis * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2, 310*25c28e83SPiotr Jasiukajtis * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and 311*25c28e83SPiotr Jasiukajtis * T1(n) = T1[2n,2n+1] = n*log(2)-1, 312*25c28e83SPiotr Jasiukajtis * T2(j) = T2[2j,2j+1] = log(z[j]), 313*25c28e83SPiotr Jasiukajtis * T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15 314*25c28e83SPiotr Jasiukajtis * Note 315*25c28e83SPiotr Jasiukajtis * (1) the leading entries are truncated to 24 binary point. 316*25c28e83SPiotr Jasiukajtis * (2) Remez error for T3(s) is bounded by 2**(-136.54) 317*25c28e83SPiotr Jasiukajtis */ 318*25c28e83SPiotr Jasiukajtis static const long double T1[] = { 319*25c28e83SPiotr Jasiukajtis -1.000000000000000000000000000000000000000000e+00L, 320*25c28e83SPiotr Jasiukajtis +0.000000000000000000000000000000000000000000e+00L, 321*25c28e83SPiotr Jasiukajtis -3.068528175354003906250000000000000000000000e-01L, 322*25c28e83SPiotr Jasiukajtis -1.904654299957767878541823431924500011926579e-09L, 323*25c28e83SPiotr Jasiukajtis +3.862943053245544433593750000000000000000000e-01L, 324*25c28e83SPiotr Jasiukajtis +5.579533617547508924291635313615100141107647e-08L, 325*25c28e83SPiotr Jasiukajtis +1.079441487789154052734375000000000000000000e+00L, 326*25c28e83SPiotr Jasiukajtis +5.389068187551732136437452970422650211661470e-08L, 327*25c28e83SPiotr Jasiukajtis +1.772588670253753662109375000000000000000000e+00L, 328*25c28e83SPiotr Jasiukajtis +5.198602757555955348583270627230200282215294e-08L, 329*25c28e83SPiotr Jasiukajtis +2.465735852718353271484375000000000000000000e+00L, 330*25c28e83SPiotr Jasiukajtis +5.008137327560178560729088284037750352769117e-08L, 331*25c28e83SPiotr Jasiukajtis +3.158883035182952880859375000000000000000000e+00L, 332*25c28e83SPiotr Jasiukajtis +4.817671897564401772874905940845299849351090e-08L, 333*25c28e83SPiotr Jasiukajtis +3.852030217647552490234375000000000000000000e+00L, 334*25c28e83SPiotr Jasiukajtis +4.627206467568624985020723597652849919904913e-08L, 335*25c28e83SPiotr Jasiukajtis +4.545177400112152099609375000000000000000000e+00L, 336*25c28e83SPiotr Jasiukajtis +4.436741037572848197166541254460399990458737e-08L, 337*25c28e83SPiotr Jasiukajtis +5.238324582576751708984375000000000000000000e+00L, 338*25c28e83SPiotr Jasiukajtis +4.246275607577071409312358911267950061012560e-08L, 339*25c28e83SPiotr Jasiukajtis +5.931471765041351318359375000000000000000000e+00L, 340*25c28e83SPiotr Jasiukajtis +4.055810177581294621458176568075500131566384e-08L, 341*25c28e83SPiotr Jasiukajtis }; 342*25c28e83SPiotr Jasiukajtis 343*25c28e83SPiotr Jasiukajtis /* 344*25c28e83SPiotr Jasiukajtis * T2[2i,2i+1] = log(1+i/64+1/128) 345*25c28e83SPiotr Jasiukajtis */ 346*25c28e83SPiotr Jasiukajtis static const long double T2[] = { 347*25c28e83SPiotr Jasiukajtis +7.7821016311645507812500000000000000000000e-03L, 348*25c28e83SPiotr Jasiukajtis +3.8810890398166212900061136763678127453570e-08L, 349*25c28e83SPiotr Jasiukajtis +2.3167014122009277343750000000000000000000e-02L, 350*25c28e83SPiotr Jasiukajtis +4.5159525100885049160962289916579411752759e-08L, 351*25c28e83SPiotr Jasiukajtis +3.8318812847137451171875000000000000000000e-02L, 352*25c28e83SPiotr Jasiukajtis +5.1454999148021880325123797290345960518164e-08L, 353*25c28e83SPiotr Jasiukajtis +5.3244471549987792968750000000000000000000e-02L, 354*25c28e83SPiotr Jasiukajtis +4.2968824489897120193786528776939573415076e-08L, 355*25c28e83SPiotr Jasiukajtis +6.7950606346130371093750000000000000000000e-02L, 356*25c28e83SPiotr Jasiukajtis +5.5562377378300815277772629414034632394030e-08L, 357*25c28e83SPiotr Jasiukajtis +8.2443654537200927734375000000000000000000e-02L, 358*25c28e83SPiotr Jasiukajtis +1.4673873663533785068668307805914095366600e-08L, 359*25c28e83SPiotr Jasiukajtis +9.6729576587677001953125000000000000000000e-02L, 360*25c28e83SPiotr Jasiukajtis +4.9870874110342446056487463437015041543346e-08L, 361*25c28e83SPiotr Jasiukajtis +1.1081433296203613281250000000000000000000e-01L, 362*25c28e83SPiotr Jasiukajtis +3.3378253981382306169323211928098474801099e-08L, 363*25c28e83SPiotr Jasiukajtis +1.2470346689224243164062500000000000000000e-01L, 364*25c28e83SPiotr Jasiukajtis +1.1608714804222781515380863268491613205318e-08L, 365*25c28e83SPiotr Jasiukajtis +1.3840228319168090820312500000000000000000e-01L, 366*25c28e83SPiotr Jasiukajtis +3.9667438227482200873601649187393160823607e-08L, 367*25c28e83SPiotr Jasiukajtis +1.5191602706909179687500000000000000000000e-01L, 368*25c28e83SPiotr Jasiukajtis +1.4956750178196803424896884511327584958252e-08L, 369*25c28e83SPiotr Jasiukajtis +1.6524952650070190429687500000000000000000e-01L, 370*25c28e83SPiotr Jasiukajtis +4.6394605258578736449277240313729237989366e-08L, 371*25c28e83SPiotr Jasiukajtis +1.7840760946273803710937500000000000000000e-01L, 372*25c28e83SPiotr Jasiukajtis +4.8010080260010025241510941968354682199540e-08L, 373*25c28e83SPiotr Jasiukajtis +1.9139480590820312500000000000000000000000e-01L, 374*25c28e83SPiotr Jasiukajtis +4.7091426329609298807561308873447039132856e-08L, 375*25c28e83SPiotr Jasiukajtis +2.0421552658081054687500000000000000000000e-01L, 376*25c28e83SPiotr Jasiukajtis +1.4847880344628820386196239272213742113867e-08L, 377*25c28e83SPiotr Jasiukajtis +2.1687388420104980468750000000000000000000e-01L, 378*25c28e83SPiotr Jasiukajtis +5.4099564554931589525744347498478964801484e-08L, 379*25c28e83SPiotr Jasiukajtis +2.2937405109405517578125000000000000000000e-01L, 380*25c28e83SPiotr Jasiukajtis +4.9970790654210230725046139871550961365282e-08L, 381*25c28e83SPiotr Jasiukajtis +2.4171990156173706054687500000000000000000e-01L, 382*25c28e83SPiotr Jasiukajtis +3.5325408107597432515913513900103385655073e-08L, 383*25c28e83SPiotr Jasiukajtis +2.5391519069671630859375000000000000000000e-01L, 384*25c28e83SPiotr Jasiukajtis +1.9284247135543573297906606667466299224747e-08L, 385*25c28e83SPiotr Jasiukajtis +2.6596349477767944335937500000000000000000e-01L, 386*25c28e83SPiotr Jasiukajtis +5.3719458497979750926537543389268821141517e-08L, 387*25c28e83SPiotr Jasiukajtis +2.7786844968795776367187500000000000000000e-01L, 388*25c28e83SPiotr Jasiukajtis +1.3154985425144750329234012330820349974537e-09L, 389*25c28e83SPiotr Jasiukajtis +2.8963327407836914062500000000000000000000e-01L, 390*25c28e83SPiotr Jasiukajtis +1.8504673536253893055525668970003860369760e-08L, 391*25c28e83SPiotr Jasiukajtis +3.0126130580902099609375000000000000000000e-01L, 392*25c28e83SPiotr Jasiukajtis +2.4769140784919125538233755492657352680723e-08L, 393*25c28e83SPiotr Jasiukajtis +3.1275570392608642578125000000000000000000e-01L, 394*25c28e83SPiotr Jasiukajtis +6.0778104626049965596883190321597861455475e-09L, 395*25c28e83SPiotr Jasiukajtis +3.2411944866180419921875000000000000000000e-01L, 396*25c28e83SPiotr Jasiukajtis +1.9992407776871920760434987352182336158873e-08L, 397*25c28e83SPiotr Jasiukajtis +3.3535552024841308593750000000000000000000e-01L, 398*25c28e83SPiotr Jasiukajtis +2.1672724744319679579814166199074433006807e-08L, 399*25c28e83SPiotr Jasiukajtis +3.4646672010421752929687500000000000000000e-01L, 400*25c28e83SPiotr Jasiukajtis +4.7241991051621587188425772950711830538414e-08L, 401*25c28e83SPiotr Jasiukajtis +3.5745584964752197265625000000000000000000e-01L, 402*25c28e83SPiotr Jasiukajtis +3.9274281801569759490140904474434669956562e-08L, 403*25c28e83SPiotr Jasiukajtis +3.6832553148269653320312500000000000000000e-01L, 404*25c28e83SPiotr Jasiukajtis +2.9676011119845105154050398826897178765758e-08L, 405*25c28e83SPiotr Jasiukajtis +3.7907832860946655273437500000000000000000e-01L, 406*25c28e83SPiotr Jasiukajtis +2.4325502905656478345631019858881408009210e-08L, 407*25c28e83SPiotr Jasiukajtis +3.8971674442291259765625000000000000000000e-01L, 408*25c28e83SPiotr Jasiukajtis +6.7171126157142136040035208670510556529487e-09L, 409*25c28e83SPiotr Jasiukajtis +4.0024316310882568359375000000000000000000e-01L, 410*25c28e83SPiotr Jasiukajtis +1.0181870233355751019951311700799406124957e-09L, 411*25c28e83SPiotr Jasiukajtis +4.1065990924835205078125000000000000000000e-01L, 412*25c28e83SPiotr Jasiukajtis +1.5736916335153056203175822787661567534220e-08L, 413*25c28e83SPiotr Jasiukajtis +4.2096924781799316406250000000000000000000e-01L, 414*25c28e83SPiotr Jasiukajtis +4.6826136472066367161506795972449857268707e-08L, 415*25c28e83SPiotr Jasiukajtis +4.3117344379425048828125000000000000000000e-01L, 416*25c28e83SPiotr Jasiukajtis +2.1024120852577922478955594998480144051225e-08L, 417*25c28e83SPiotr Jasiukajtis +4.4127452373504638671875000000000000000000e-01L, 418*25c28e83SPiotr Jasiukajtis +3.7069828842770746441661301225362605528786e-08L, 419*25c28e83SPiotr Jasiukajtis +4.5127463340759277343750000000000000000000e-01L, 420*25c28e83SPiotr Jasiukajtis +1.0731865811707192383079012478685922879010e-08L, 421*25c28e83SPiotr Jasiukajtis +4.6117568016052246093750000000000000000000e-01L, 422*25c28e83SPiotr Jasiukajtis +3.4961647705430499925597855358603099030515e-08L, 423*25c28e83SPiotr Jasiukajtis +4.7097969055175781250000000000000000000000e-01L, 424*25c28e83SPiotr Jasiukajtis +2.4667033200046897856056359251373510964634e-08L, 425*25c28e83SPiotr Jasiukajtis +4.8068851232528686523437500000000000000000e-01L, 426*25c28e83SPiotr Jasiukajtis +1.7020465042442243455448011551208861216878e-08L, 427*25c28e83SPiotr Jasiukajtis +4.9030393362045288085937500000000000000000e-01L, 428*25c28e83SPiotr Jasiukajtis +5.4424740957290971159645746860530583309571e-08L, 429*25c28e83SPiotr Jasiukajtis +4.9982786178588867187500000000000000000000e-01L, 430*25c28e83SPiotr Jasiukajtis +7.7705606579463314152470441415126573566105e-09L, 431*25c28e83SPiotr Jasiukajtis +5.0926184654235839843750000000000000000000e-01L, 432*25c28e83SPiotr Jasiukajtis +5.5247449548366574919228323824878565745713e-08L, 433*25c28e83SPiotr Jasiukajtis +5.1860773563385009765625000000000000000000e-01L, 434*25c28e83SPiotr Jasiukajtis +2.8574195534496726996364798698556235730848e-08L, 435*25c28e83SPiotr Jasiukajtis +5.2786707878112792968750000000000000000000e-01L, 436*25c28e83SPiotr Jasiukajtis +1.0839714455426392217778300963558522088193e-08L, 437*25c28e83SPiotr Jasiukajtis +5.3704142570495605468750000000000000000000e-01L, 438*25c28e83SPiotr Jasiukajtis +4.0191927599879229244153832299023744345999e-08L, 439*25c28e83SPiotr Jasiukajtis +5.4613238573074340820312500000000000000000e-01L, 440*25c28e83SPiotr Jasiukajtis +5.1867392242179272209231209163864971792889e-08L, 441*25c28e83SPiotr Jasiukajtis +5.5514144897460937500000000000000000000000e-01L, 442*25c28e83SPiotr Jasiukajtis +5.8565892217715480359515904050170125743178e-08L, 443*25c28e83SPiotr Jasiukajtis +5.6407010555267333984375000000000000000000e-01L, 444*25c28e83SPiotr Jasiukajtis +3.2732129626227634290090190711817681692354e-08L, 445*25c28e83SPiotr Jasiukajtis +5.7291972637176513671875000000000000000000e-01L, 446*25c28e83SPiotr Jasiukajtis +2.7190020372374006726626261068626400393936e-08L, 447*25c28e83SPiotr Jasiukajtis +5.8169168233871459960937500000000000000000e-01L, 448*25c28e83SPiotr Jasiukajtis +5.7295907882911235753725372340709967597394e-08L, 449*25c28e83SPiotr Jasiukajtis +5.9038740396499633789062500000000000000000e-01L, 450*25c28e83SPiotr Jasiukajtis +4.2637180036751291708123598757577783615014e-08L, 451*25c28e83SPiotr Jasiukajtis +5.9900814294815063476562500000000000000000e-01L, 452*25c28e83SPiotr Jasiukajtis +4.6697932764615975024461651502060474048774e-08L, 453*25c28e83SPiotr Jasiukajtis +6.0755521059036254882812500000000000000000e-01L, 454*25c28e83SPiotr Jasiukajtis +3.9634179246672960152791125371893149820625e-08L, 455*25c28e83SPiotr Jasiukajtis +6.1602985858917236328125000000000000000000e-01L, 456*25c28e83SPiotr Jasiukajtis +1.8626341656366315928196700650292529688219e-08L, 457*25c28e83SPiotr Jasiukajtis +6.2443327903747558593750000000000000000000e-01L, 458*25c28e83SPiotr Jasiukajtis +8.9744179151050387440546731199093039879228e-09L, 459*25c28e83SPiotr Jasiukajtis +6.3276666402816772460937500000000000000000e-01L, 460*25c28e83SPiotr Jasiukajtis +5.5428701049364114685035797584887586099726e-09L, 461*25c28e83SPiotr Jasiukajtis +6.4103114604949951171875000000000000000000e-01L, 462*25c28e83SPiotr Jasiukajtis +3.3371431779336851334405392546708949047361e-08L, 463*25c28e83SPiotr Jasiukajtis +6.4922791719436645507812500000000000000000e-01L, 464*25c28e83SPiotr Jasiukajtis +2.9430743363812714969905311122271269100885e-08L, 465*25c28e83SPiotr Jasiukajtis +6.5735805034637451171875000000000000000000e-01L, 466*25c28e83SPiotr Jasiukajtis +2.2361985518423140023245936165514147093250e-08L, 467*25c28e83SPiotr Jasiukajtis +6.6542261838912963867187500000000000000000e-01L, 468*25c28e83SPiotr Jasiukajtis +1.4155960810278217610006660181148303091649e-08L, 469*25c28e83SPiotr Jasiukajtis +6.7342263460159301757812500000000000000000e-01L, 470*25c28e83SPiotr Jasiukajtis +4.0610573702719835388801017264750843477878e-08L, 471*25c28e83SPiotr Jasiukajtis +6.8135917186737060546875000000000000000000e-01L, 472*25c28e83SPiotr Jasiukajtis +5.2940532463479321559568089441735584156689e-08L, 473*25c28e83SPiotr Jasiukajtis +6.8923324346542358398437500000000000000000e-01L, 474*25c28e83SPiotr Jasiukajtis +3.7773385396340539337814603903232796216537e-08L, 475*25c28e83SPiotr Jasiukajtis }; 476*25c28e83SPiotr Jasiukajtis 477*25c28e83SPiotr Jasiukajtis /* 478*25c28e83SPiotr Jasiukajtis * S[j],S_trail[j] = 2**(j/32.) for the final computation of exp(t+w) 479*25c28e83SPiotr Jasiukajtis */ 480*25c28e83SPiotr Jasiukajtis static const long double S[] = { 481*25c28e83SPiotr Jasiukajtis #if defined(__x86) 482*25c28e83SPiotr Jasiukajtis +1.0000000000000000000000000e+00L, 483*25c28e83SPiotr Jasiukajtis +1.0218971486541166782081522e+00L, 484*25c28e83SPiotr Jasiukajtis +1.0442737824274138402382006e+00L, 485*25c28e83SPiotr Jasiukajtis +1.0671404006768236181297224e+00L, 486*25c28e83SPiotr Jasiukajtis +1.0905077326652576591003302e+00L, 487*25c28e83SPiotr Jasiukajtis +1.1143867425958925362894369e+00L, 488*25c28e83SPiotr Jasiukajtis +1.1387886347566916536971221e+00L, 489*25c28e83SPiotr Jasiukajtis +1.1637248587775775137938619e+00L, 490*25c28e83SPiotr Jasiukajtis +1.1892071150027210666875674e+00L, 491*25c28e83SPiotr Jasiukajtis +1.2152473599804688780476325e+00L, 492*25c28e83SPiotr Jasiukajtis +1.2418578120734840485256747e+00L, 493*25c28e83SPiotr Jasiukajtis +1.2690509571917332224885722e+00L, 494*25c28e83SPiotr Jasiukajtis +1.2968395546510096659215822e+00L, 495*25c28e83SPiotr Jasiukajtis +1.3252366431597412945939118e+00L, 496*25c28e83SPiotr Jasiukajtis +1.3542555469368927282668852e+00L, 497*25c28e83SPiotr Jasiukajtis +1.3839098819638319548151403e+00L, 498*25c28e83SPiotr Jasiukajtis +1.4142135623730950487637881e+00L, 499*25c28e83SPiotr Jasiukajtis +1.4451808069770466200253470e+00L, 500*25c28e83SPiotr Jasiukajtis +1.4768261459394993113155431e+00L, 501*25c28e83SPiotr Jasiukajtis +1.5091644275934227397133885e+00L, 502*25c28e83SPiotr Jasiukajtis +1.5422108254079408235859630e+00L, 503*25c28e83SPiotr Jasiukajtis +1.5759808451078864864006862e+00L, 504*25c28e83SPiotr Jasiukajtis +1.6104903319492543080837174e+00L, 505*25c28e83SPiotr Jasiukajtis +1.6457554781539648445110730e+00L, 506*25c28e83SPiotr Jasiukajtis +1.6817928305074290860378350e+00L, 507*25c28e83SPiotr Jasiukajtis +1.7186192981224779156032914e+00L, 508*25c28e83SPiotr Jasiukajtis +1.7562521603732994831094730e+00L, 509*25c28e83SPiotr Jasiukajtis +1.7947090750031071864148413e+00L, 510*25c28e83SPiotr Jasiukajtis +1.8340080864093424633989166e+00L, 511*25c28e83SPiotr Jasiukajtis +1.8741676341102999013002103e+00L, 512*25c28e83SPiotr Jasiukajtis +1.9152065613971472938202589e+00L, 513*25c28e83SPiotr Jasiukajtis +1.9571441241754002689657438e+00L, 514*25c28e83SPiotr Jasiukajtis #else 515*25c28e83SPiotr Jasiukajtis +1.00000000000000000000000000000000000e+00L, 516*25c28e83SPiotr Jasiukajtis +1.02189714865411667823448013478329942e+00L, 517*25c28e83SPiotr Jasiukajtis +1.04427378242741384032196647873992910e+00L, 518*25c28e83SPiotr Jasiukajtis +1.06714040067682361816952112099280918e+00L, 519*25c28e83SPiotr Jasiukajtis +1.09050773266525765920701065576070789e+00L, 520*25c28e83SPiotr Jasiukajtis +1.11438674259589253630881295691960313e+00L, 521*25c28e83SPiotr Jasiukajtis +1.13878863475669165370383028384151134e+00L, 522*25c28e83SPiotr Jasiukajtis +1.16372485877757751381357359909218536e+00L, 523*25c28e83SPiotr Jasiukajtis +1.18920711500272106671749997056047593e+00L, 524*25c28e83SPiotr Jasiukajtis +1.21524735998046887811652025133879836e+00L, 525*25c28e83SPiotr Jasiukajtis +1.24185781207348404859367746872659561e+00L, 526*25c28e83SPiotr Jasiukajtis +1.26905095719173322255441908103233805e+00L, 527*25c28e83SPiotr Jasiukajtis +1.29683955465100966593375411779245118e+00L, 528*25c28e83SPiotr Jasiukajtis +1.32523664315974129462953709549872168e+00L, 529*25c28e83SPiotr Jasiukajtis +1.35425554693689272829801474014070273e+00L, 530*25c28e83SPiotr Jasiukajtis +1.38390988196383195487265952726519287e+00L, 531*25c28e83SPiotr Jasiukajtis +1.41421356237309504880168872420969798e+00L, 532*25c28e83SPiotr Jasiukajtis +1.44518080697704662003700624147167095e+00L, 533*25c28e83SPiotr Jasiukajtis +1.47682614593949931138690748037404985e+00L, 534*25c28e83SPiotr Jasiukajtis +1.50916442759342273976601955103319352e+00L, 535*25c28e83SPiotr Jasiukajtis +1.54221082540794082361229186209073479e+00L, 536*25c28e83SPiotr Jasiukajtis +1.57598084510788648645527016018190504e+00L, 537*25c28e83SPiotr Jasiukajtis +1.61049033194925430817952066735740067e+00L, 538*25c28e83SPiotr Jasiukajtis +1.64575547815396484451875672472582254e+00L, 539*25c28e83SPiotr Jasiukajtis +1.68179283050742908606225095246642969e+00L, 540*25c28e83SPiotr Jasiukajtis +1.71861929812247791562934437645631244e+00L, 541*25c28e83SPiotr Jasiukajtis +1.75625216037329948311216061937531314e+00L, 542*25c28e83SPiotr Jasiukajtis +1.79470907500310718642770324212778174e+00L, 543*25c28e83SPiotr Jasiukajtis +1.83400808640934246348708318958828892e+00L, 544*25c28e83SPiotr Jasiukajtis +1.87416763411029990132999894995444645e+00L, 545*25c28e83SPiotr Jasiukajtis +1.91520656139714729387261127029583086e+00L, 546*25c28e83SPiotr Jasiukajtis +1.95714412417540026901832225162687149e+00L, 547*25c28e83SPiotr Jasiukajtis #endif 548*25c28e83SPiotr Jasiukajtis }; 549*25c28e83SPiotr Jasiukajtis static const long double S_trail[] = { 550*25c28e83SPiotr Jasiukajtis #if defined(__x86) 551*25c28e83SPiotr Jasiukajtis +0.0000000000000000000000000e+00L, 552*25c28e83SPiotr Jasiukajtis +2.6327965667180882569382524e-20L, 553*25c28e83SPiotr Jasiukajtis +8.3765863521895191129661899e-20L, 554*25c28e83SPiotr Jasiukajtis +3.9798705777454504249209575e-20L, 555*25c28e83SPiotr Jasiukajtis +1.0668046596651558640993042e-19L, 556*25c28e83SPiotr Jasiukajtis +1.9376009847285360448117114e-20L, 557*25c28e83SPiotr Jasiukajtis +6.7081819456112953751277576e-21L, 558*25c28e83SPiotr Jasiukajtis +1.9711680502629186462729727e-20L, 559*25c28e83SPiotr Jasiukajtis +2.9932584438449523689104569e-20L, 560*25c28e83SPiotr Jasiukajtis +6.8887754153039109411061914e-20L, 561*25c28e83SPiotr Jasiukajtis +6.8002718741225378942847820e-20L, 562*25c28e83SPiotr Jasiukajtis +6.5846917376975403439742349e-20L, 563*25c28e83SPiotr Jasiukajtis +1.2171958727511372194876001e-20L, 564*25c28e83SPiotr Jasiukajtis +3.5625253228704087115438260e-20L, 565*25c28e83SPiotr Jasiukajtis +3.1129551559077560956309179e-20L, 566*25c28e83SPiotr Jasiukajtis +5.7519192396164779846216492e-20L, 567*25c28e83SPiotr Jasiukajtis +3.7900651177865141593101239e-20L, 568*25c28e83SPiotr Jasiukajtis +1.1659262405698741798080115e-20L, 569*25c28e83SPiotr Jasiukajtis +7.1364385105284695967172478e-20L, 570*25c28e83SPiotr Jasiukajtis +5.2631003710812203588788949e-20L, 571*25c28e83SPiotr Jasiukajtis +2.6328853788732632868460580e-20L, 572*25c28e83SPiotr Jasiukajtis +5.4583950085438242788190141e-20L, 573*25c28e83SPiotr Jasiukajtis +9.5803254376938269960718656e-20L, 574*25c28e83SPiotr Jasiukajtis +7.6837733983874245823512279e-21L, 575*25c28e83SPiotr Jasiukajtis +2.4415965910835093824202087e-20L, 576*25c28e83SPiotr Jasiukajtis +2.6052966871016580981769728e-20L, 577*25c28e83SPiotr Jasiukajtis +2.6876456344632553875309579e-21L, 578*25c28e83SPiotr Jasiukajtis +1.2861930155613700201703279e-20L, 579*25c28e83SPiotr Jasiukajtis +8.8166633394037485606572294e-20L, 580*25c28e83SPiotr Jasiukajtis +2.9788615389580190940837037e-20L, 581*25c28e83SPiotr Jasiukajtis +5.2352341619805098677422139e-20L, 582*25c28e83SPiotr Jasiukajtis +5.2578463064010463732242363e-20L, 583*25c28e83SPiotr Jasiukajtis #else 584*25c28e83SPiotr Jasiukajtis +0.00000000000000000000000000000000000e+00L, 585*25c28e83SPiotr Jasiukajtis +1.80506787420330954745573333054573786e-35L, 586*25c28e83SPiotr Jasiukajtis -9.37452029228042742195756741973083214e-35L, 587*25c28e83SPiotr Jasiukajtis -1.59696844729275877071290963023149997e-35L, 588*25c28e83SPiotr Jasiukajtis +9.11249341012502297851168610167248666e-35L, 589*25c28e83SPiotr Jasiukajtis -6.50422820697854828723037477525938871e-35L, 590*25c28e83SPiotr Jasiukajtis -8.14846884452585113732569176748815532e-35L, 591*25c28e83SPiotr Jasiukajtis -5.06621457672180031337233074514290335e-35L, 592*25c28e83SPiotr Jasiukajtis -1.35983097468881697374987563824591912e-35L, 593*25c28e83SPiotr Jasiukajtis +9.49742763556319647030771056643324660e-35L, 594*25c28e83SPiotr Jasiukajtis -3.28317052317699860161506596533391526e-36L, 595*25c28e83SPiotr Jasiukajtis -5.01723570938719041029018653045842895e-35L, 596*25c28e83SPiotr Jasiukajtis -2.39147479768910917162283430160264014e-35L, 597*25c28e83SPiotr Jasiukajtis -8.35057135763390881529889073794408385e-36L, 598*25c28e83SPiotr Jasiukajtis +7.03675688907326504242173719067187644e-35L, 599*25c28e83SPiotr Jasiukajtis -5.18248485306464645753689301856695619e-35L, 600*25c28e83SPiotr Jasiukajtis +9.42224254862183206569211673639406488e-35L, 601*25c28e83SPiotr Jasiukajtis -3.96750082539886230916730613021641828e-35L, 602*25c28e83SPiotr Jasiukajtis +7.14352899156330061452327361509276724e-35L, 603*25c28e83SPiotr Jasiukajtis +1.15987125286798512424651783410044433e-35L, 604*25c28e83SPiotr Jasiukajtis +4.69693347835811549530973921320187447e-35L, 605*25c28e83SPiotr Jasiukajtis -3.38651317599500471079924198499981917e-35L, 606*25c28e83SPiotr Jasiukajtis -8.58731877429824706886865593510387445e-35L, 607*25c28e83SPiotr Jasiukajtis -9.60595154874935050318549936224606909e-35L, 608*25c28e83SPiotr Jasiukajtis +9.60973393212801278450755869714178581e-35L, 609*25c28e83SPiotr Jasiukajtis +6.37839792144002843924476144978084855e-35L, 610*25c28e83SPiotr Jasiukajtis +7.79243078569586424945646112516927770e-35L, 611*25c28e83SPiotr Jasiukajtis +7.36133776758845652413193083663393220e-35L, 612*25c28e83SPiotr Jasiukajtis -6.47299514791334723003521457561217053e-35L, 613*25c28e83SPiotr Jasiukajtis +8.58747441795369869427879806229522962e-35L, 614*25c28e83SPiotr Jasiukajtis +2.37181542282517483569165122830269098e-35L, 615*25c28e83SPiotr Jasiukajtis -3.02689168209611877300459737342190031e-37L, 616*25c28e83SPiotr Jasiukajtis #endif 617*25c28e83SPiotr Jasiukajtis }; 618*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 619*25c28e83SPiotr Jasiukajtis 620*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 621*25c28e83SPiotr Jasiukajtis /* 622*25c28e83SPiotr Jasiukajtis * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula 623*25c28e83SPiotr Jasiukajtis * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x)) 624*25c28e83SPiotr Jasiukajtis * = L1 + L2 + L3, 625*25c28e83SPiotr Jasiukajtis */ 626*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 627*25c28e83SPiotr Jasiukajtis static struct LDouble 628*25c28e83SPiotr Jasiukajtis large_gam(long double x, int *m) { 629*25c28e83SPiotr Jasiukajtis long double z, t1, t2, t3, z2, t5, w, y, u, r, v; 630*25c28e83SPiotr Jasiukajtis long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L; 631*25c28e83SPiotr Jasiukajtis int n2, j2, k, ix, j, i; 632*25c28e83SPiotr Jasiukajtis struct LDouble zz; 633*25c28e83SPiotr Jasiukajtis long double u2, ss_h, ss_l, r_h, w_h, w_l, t4; 634*25c28e83SPiotr Jasiukajtis 635*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 636*25c28e83SPiotr Jasiukajtis /* 637*25c28e83SPiotr Jasiukajtis * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details) 638*25c28e83SPiotr Jasiukajtis * 639*25c28e83SPiotr Jasiukajtis * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2, 640*25c28e83SPiotr Jasiukajtis * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and 641*25c28e83SPiotr Jasiukajtis * T1(n) = T1[2n,2n+1] = n*log(2)-1, 642*25c28e83SPiotr Jasiukajtis * T2(j) = T2[2j,2j+1] = log(z[j]), 643*25c28e83SPiotr Jasiukajtis * T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15 644*25c28e83SPiotr Jasiukajtis * Note 645*25c28e83SPiotr Jasiukajtis * (1) the leading entries are truncated to 24 binary point. 646*25c28e83SPiotr Jasiukajtis * (2) Remez error for T3(s) is bounded by 2**(-72.4) 647*25c28e83SPiotr Jasiukajtis * 2**(-24) 648*25c28e83SPiotr Jasiukajtis * _________V___________________ 649*25c28e83SPiotr Jasiukajtis * T1(n): |_________|___________________| 650*25c28e83SPiotr Jasiukajtis * _______ ______________________ 651*25c28e83SPiotr Jasiukajtis * T2(j): |_______|______________________| 652*25c28e83SPiotr Jasiukajtis * ____ _______________________ 653*25c28e83SPiotr Jasiukajtis * 2s: |____|_______________________| 654*25c28e83SPiotr Jasiukajtis * __________________________ 655*25c28e83SPiotr Jasiukajtis * + T3(s)-2s: |__________________________| 656*25c28e83SPiotr Jasiukajtis * ------------------------------------------- 657*25c28e83SPiotr Jasiukajtis * [leading] + [Trailing] 658*25c28e83SPiotr Jasiukajtis */ 659*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 660*25c28e83SPiotr Jasiukajtis ix = H0_WORD(x); 661*25c28e83SPiotr Jasiukajtis n2 = (ix >> 16) - 0x3fff; /* exponent of x, range:3-10 */ 662*25c28e83SPiotr Jasiukajtis y = scalbnl(x, -n2); /* y = scale x to [1,2] */ 663*25c28e83SPiotr Jasiukajtis n2 += n2; /* 2n */ 664*25c28e83SPiotr Jasiukajtis j = (ix >> 10) & 0x3f; /* j */ 665*25c28e83SPiotr Jasiukajtis z = 1.0078125L + (long double) j * 0.015625L; /* z[j]=1+j/64+1/128 */ 666*25c28e83SPiotr Jasiukajtis j2 = j + j; 667*25c28e83SPiotr Jasiukajtis t1 = y + z; 668*25c28e83SPiotr Jasiukajtis t2 = y - z; 669*25c28e83SPiotr Jasiukajtis r = one / t1; 670*25c28e83SPiotr Jasiukajtis u = r * t2; /* u = (y-z)/(y+z) */ 671*25c28e83SPiotr Jasiukajtis t1 = CHOPPED(t1); 672*25c28e83SPiotr Jasiukajtis t4 = T2[j2 + 1] + T1[n2 + 1]; 673*25c28e83SPiotr Jasiukajtis z2 = u * u; 674*25c28e83SPiotr Jasiukajtis k = H0_WORD(u) & 0x7fffffff; 675*25c28e83SPiotr Jasiukajtis t3 = T2[j2] + T1[n2]; 676*25c28e83SPiotr Jasiukajtis for (t5 = T3[6], i = 5; i >= 0; i--) 677*25c28e83SPiotr Jasiukajtis t5 = z2 * t5 + T3[i]; 678*25c28e83SPiotr Jasiukajtis if ((k >> 16) < 0x3fec) { /* |u|<2**-19 */ 679*25c28e83SPiotr Jasiukajtis t2 = t4 + u * (two + z2 * t5); 680*25c28e83SPiotr Jasiukajtis } else { 681*25c28e83SPiotr Jasiukajtis t5 = t4 + (u * z2) * t5; 682*25c28e83SPiotr Jasiukajtis u2 = u + u; 683*25c28e83SPiotr Jasiukajtis v = (long double) ((int) (u2 * t24)) * p24; 684*25c28e83SPiotr Jasiukajtis t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z))); 685*25c28e83SPiotr Jasiukajtis t3 += v; 686*25c28e83SPiotr Jasiukajtis } 687*25c28e83SPiotr Jasiukajtis ss_h = CHOPPED((t2 + t3)); 688*25c28e83SPiotr Jasiukajtis ss_l = t2 - (ss_h - t3); 689*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 690*25c28e83SPiotr Jasiukajtis /* 691*25c28e83SPiotr Jasiukajtis * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2))) 692*25c28e83SPiotr Jasiukajtis * where ss = log(x) - 1 in already in extra precision 693*25c28e83SPiotr Jasiukajtis */ 694*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 695*25c28e83SPiotr Jasiukajtis z = one / x; 696*25c28e83SPiotr Jasiukajtis r = x - half; 697*25c28e83SPiotr Jasiukajtis r_h = CHOPPED((r)); 698*25c28e83SPiotr Jasiukajtis w_h = r_h * ss_h + hln2pim1_h; 699*25c28e83SPiotr Jasiukajtis z2 = z * z; 700*25c28e83SPiotr Jasiukajtis w = (r - r_h) * ss_h + r * ss_l; 701*25c28e83SPiotr Jasiukajtis t1 = GP[19]; 702*25c28e83SPiotr Jasiukajtis for (i = 18; i > 0; i--) 703*25c28e83SPiotr Jasiukajtis t1 = z2 * t1 + GP[i]; 704*25c28e83SPiotr Jasiukajtis w += hln2pim1_l; 705*25c28e83SPiotr Jasiukajtis w_l = z * (GP[0] + z2 * t1) + w; 706*25c28e83SPiotr Jasiukajtis k = (int) ((w_h + w_l) * invln2_32 + half); 707*25c28e83SPiotr Jasiukajtis 708*25c28e83SPiotr Jasiukajtis /* compute the exponential of w_h+w_l */ 709*25c28e83SPiotr Jasiukajtis 710*25c28e83SPiotr Jasiukajtis j = k & 0x1f; 711*25c28e83SPiotr Jasiukajtis *m = k >> 5; 712*25c28e83SPiotr Jasiukajtis t3 = (long double) k; 713*25c28e83SPiotr Jasiukajtis 714*25c28e83SPiotr Jasiukajtis /* perform w - k*ln2_32 (represent as w_h - w_l) */ 715*25c28e83SPiotr Jasiukajtis t1 = w_h - t3 * ln2_32hi; 716*25c28e83SPiotr Jasiukajtis t2 = t3 * ln2_32lo; 717*25c28e83SPiotr Jasiukajtis w = t2 - w_l; 718*25c28e83SPiotr Jasiukajtis w_h = t1 - w; 719*25c28e83SPiotr Jasiukajtis w_l = w - (t1 - w_h); 720*25c28e83SPiotr Jasiukajtis 721*25c28e83SPiotr Jasiukajtis /* compute exp(w_h-w_l) */ 722*25c28e83SPiotr Jasiukajtis z = w_h - w_l; 723*25c28e83SPiotr Jasiukajtis for (t1 = Et[10], i = 9; i >= 0; i--) 724*25c28e83SPiotr Jasiukajtis t1 = z * t1 + Et[i]; 725*25c28e83SPiotr Jasiukajtis t3 = w_h - (w_l - (z * z) * t1); /* t3 = expm1(z) */ 726*25c28e83SPiotr Jasiukajtis zz.l = S_trail[j] * (one + t3) + S[j] * t3; 727*25c28e83SPiotr Jasiukajtis zz.h = S[j]; 728*25c28e83SPiotr Jasiukajtis return (zz); 729*25c28e83SPiotr Jasiukajtis } 730*25c28e83SPiotr Jasiukajtis 731*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 732*25c28e83SPiotr Jasiukajtis /* 733*25c28e83SPiotr Jasiukajtis * kpsin(x)= sin(pi*x)/pi 734*25c28e83SPiotr Jasiukajtis * 3 5 7 9 11 27 735*25c28e83SPiotr Jasiukajtis * = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x + ... + ks[12]*x 736*25c28e83SPiotr Jasiukajtis */ 737*25c28e83SPiotr Jasiukajtis static const long double ks[] = { 738*25c28e83SPiotr Jasiukajtis -1.64493406684822643647241516664602518705158902870e+0000L, 739*25c28e83SPiotr Jasiukajtis +8.11742425283353643637002772405874238094995726160e-0001L, 740*25c28e83SPiotr Jasiukajtis -1.90751824122084213696472111835337366232282723933e-0001L, 741*25c28e83SPiotr Jasiukajtis +2.61478478176548005046532613563241288115395517084e-0002L, 742*25c28e83SPiotr Jasiukajtis -2.34608103545582363750893072647117829448016479971e-0003L, 743*25c28e83SPiotr Jasiukajtis +1.48428793031071003684606647212534027556262040158e-0004L, 744*25c28e83SPiotr Jasiukajtis -6.97587366165638046518462722252768122615952898698e-0006L, 745*25c28e83SPiotr Jasiukajtis +2.53121740413702536928659271747187500934840057929e-0007L, 746*25c28e83SPiotr Jasiukajtis -7.30471182221385990397683641695766121301933621956e-0009L, 747*25c28e83SPiotr Jasiukajtis +1.71653847451163495739958249695549313987973589884e-0010L, 748*25c28e83SPiotr Jasiukajtis -3.34813314714560776122245796929054813458341420565e-0012L, 749*25c28e83SPiotr Jasiukajtis +5.50724992262622033449487808306969135431411753047e-0014L, 750*25c28e83SPiotr Jasiukajtis -7.67678132753577998601234393215802221104236979928e-0016L, 751*25c28e83SPiotr Jasiukajtis }; 752*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 753*25c28e83SPiotr Jasiukajtis 754*25c28e83SPiotr Jasiukajtis /* 755*25c28e83SPiotr Jasiukajtis * assume x is not tiny and positive 756*25c28e83SPiotr Jasiukajtis */ 757*25c28e83SPiotr Jasiukajtis static struct LDouble 758*25c28e83SPiotr Jasiukajtis kpsin(long double x) { 759*25c28e83SPiotr Jasiukajtis long double z, t1, t2; 760*25c28e83SPiotr Jasiukajtis struct LDouble xx; 761*25c28e83SPiotr Jasiukajtis int i; 762*25c28e83SPiotr Jasiukajtis 763*25c28e83SPiotr Jasiukajtis z = x * x; 764*25c28e83SPiotr Jasiukajtis xx.h = x; 765*25c28e83SPiotr Jasiukajtis for (t2 = ks[12], i = 11; i > 0; i--) 766*25c28e83SPiotr Jasiukajtis t2 = z * t2 + ks[i]; 767*25c28e83SPiotr Jasiukajtis t1 = z * x; 768*25c28e83SPiotr Jasiukajtis t2 *= z * t1; 769*25c28e83SPiotr Jasiukajtis xx.l = t1 * ks[0] + t2; 770*25c28e83SPiotr Jasiukajtis return (xx); 771*25c28e83SPiotr Jasiukajtis } 772*25c28e83SPiotr Jasiukajtis 773*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 774*25c28e83SPiotr Jasiukajtis /* 775*25c28e83SPiotr Jasiukajtis * kpcos(x)= cos(pi*x)/pi 776*25c28e83SPiotr Jasiukajtis * 2 4 6 8 10 12 777*25c28e83SPiotr Jasiukajtis * = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x +kc[5]*x 778*25c28e83SPiotr Jasiukajtis * 779*25c28e83SPiotr Jasiukajtis * 2 4 6 8 10 22 780*25c28e83SPiotr Jasiukajtis * = 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +...+kc[9]*x 781*25c28e83SPiotr Jasiukajtis * 782*25c28e83SPiotr Jasiukajtis * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l) 783*25c28e83SPiotr Jasiukajtis * = npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x 784*25c28e83SPiotr Jasiukajtis * = npi_2_h*x_f*x_f + npi_2_h*(x*x-x_f*x_f) + npi_2_l*x*x 785*25c28e83SPiotr Jasiukajtis * = npi_2_h*x_f*x_f + npi_2_h*(x+x_f)*(x-x_f) + npi_2_l*x*x 786*25c28e83SPiotr Jasiukajtis * Here x_f = (long double) (float)x 787*25c28e83SPiotr Jasiukajtis * Note that pi/2(in hex) = 788*25c28e83SPiotr Jasiukajtis * 1.921FB54442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 789*25c28e83SPiotr Jasiukajtis * npi_2_h = -pi/2 chopped to 25 bits = -1.921FB50000000000000000000000000 = 790*25c28e83SPiotr Jasiukajtis * -1.570796310901641845703125000000000 and 791*25c28e83SPiotr Jasiukajtis * npi_2_l = 792*25c28e83SPiotr Jasiukajtis * -0.0000004442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 = 793*25c28e83SPiotr Jasiukajtis * -.0000000158932547735281966916397514420985846996875529104874722961539 = 794*25c28e83SPiotr Jasiukajtis * -1.5893254773528196691639751442098584699687552910487472296153e-8 795*25c28e83SPiotr Jasiukajtis * 1/pi(in hex) = 796*25c28e83SPiotr Jasiukajtis * .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B 797*25c28e83SPiotr Jasiukajtis * will be splitted into: 798*25c28e83SPiotr Jasiukajtis * one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000... and 799*25c28e83SPiotr Jasiukajtis * one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B 800*25c28e83SPiotr Jasiukajtis */ 801*25c28e83SPiotr Jasiukajtis 802*25c28e83SPiotr Jasiukajtis static const long double 803*25c28e83SPiotr Jasiukajtis #if defined(__x86) 804*25c28e83SPiotr Jasiukajtis one_pi_h = 0.3183098861481994390487670898437500L, /* 31 bits */ 805*25c28e83SPiotr Jasiukajtis one_pi_l = 3.559123248900043690127872406891929148e-11L, 806*25c28e83SPiotr Jasiukajtis #else 807*25c28e83SPiotr Jasiukajtis one_pi_h = 0.31830988618379052468299050815403461456298828125L, 808*25c28e83SPiotr Jasiukajtis one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L, 809*25c28e83SPiotr Jasiukajtis #endif 810*25c28e83SPiotr Jasiukajtis npi_2_h = -1.570796310901641845703125000000000L, 811*25c28e83SPiotr Jasiukajtis npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L; 812*25c28e83SPiotr Jasiukajtis 813*25c28e83SPiotr Jasiukajtis static const long double kc[] = { 814*25c28e83SPiotr Jasiukajtis +1.29192819501249250731151312779548918765320728489e+0000L, 815*25c28e83SPiotr Jasiukajtis -4.25027339979557573976029596929319207009444090366e-0001L, 816*25c28e83SPiotr Jasiukajtis +7.49080661650990096109672954618317623888421628613e-0002L, 817*25c28e83SPiotr Jasiukajtis -8.21458866111282287985539464173976555436050215120e-0003L, 818*25c28e83SPiotr Jasiukajtis +6.14202578809529228503205255165761204750211603402e-0004L, 819*25c28e83SPiotr Jasiukajtis -3.33073432691149607007217330302595267179545908740e-0005L, 820*25c28e83SPiotr Jasiukajtis +1.36970959047832085796809745461530865597993680204e-0006L, 821*25c28e83SPiotr Jasiukajtis -4.41780774262583514450246512727201806217271097336e-0008L, 822*25c28e83SPiotr Jasiukajtis +1.14741409212381858820016567664488123478660705759e-0009L, 823*25c28e83SPiotr Jasiukajtis -2.44261236114707374558437500654381006300502749632e-0011L, 824*25c28e83SPiotr Jasiukajtis }; 825*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 826*25c28e83SPiotr Jasiukajtis 827*25c28e83SPiotr Jasiukajtis /* 828*25c28e83SPiotr Jasiukajtis * assume x is not tiny and positive 829*25c28e83SPiotr Jasiukajtis */ 830*25c28e83SPiotr Jasiukajtis static struct LDouble 831*25c28e83SPiotr Jasiukajtis kpcos(long double x) { 832*25c28e83SPiotr Jasiukajtis long double z, t1, t2, t3, t4, x4, x8; 833*25c28e83SPiotr Jasiukajtis int i; 834*25c28e83SPiotr Jasiukajtis struct LDouble xx; 835*25c28e83SPiotr Jasiukajtis 836*25c28e83SPiotr Jasiukajtis z = x * x; 837*25c28e83SPiotr Jasiukajtis xx.h = one_pi_h; 838*25c28e83SPiotr Jasiukajtis t1 = (long double) ((float) x); 839*25c28e83SPiotr Jasiukajtis x4 = z * z; 840*25c28e83SPiotr Jasiukajtis t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1); 841*25c28e83SPiotr Jasiukajtis for (i = 8, t3 = kc[9]; i >= 0; i--) 842*25c28e83SPiotr Jasiukajtis t3 = z * t3 + kc[i]; 843*25c28e83SPiotr Jasiukajtis t3 = one_pi_l + x4 * t3; 844*25c28e83SPiotr Jasiukajtis t4 = t1 * t1 * npi_2_h; 845*25c28e83SPiotr Jasiukajtis x8 = t2 + t3; 846*25c28e83SPiotr Jasiukajtis xx.l = x8 + t4; 847*25c28e83SPiotr Jasiukajtis return (xx); 848*25c28e83SPiotr Jasiukajtis } 849*25c28e83SPiotr Jasiukajtis 850*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 851*25c28e83SPiotr Jasiukajtis static const long double 852*25c28e83SPiotr Jasiukajtis /* 0.13486180573279076968979393577465291700642511139552429398233 */ 853*25c28e83SPiotr Jasiukajtis #if defined(__x86) 854*25c28e83SPiotr Jasiukajtis t0z1 = 0.1348618057327907696779385054997035808810L, 855*25c28e83SPiotr Jasiukajtis t0z1_l = 1.1855430274949336125392717150257379614654e-20L, 856*25c28e83SPiotr Jasiukajtis #else 857*25c28e83SPiotr Jasiukajtis t0z1 = 0.1348618057327907696897939357746529168654L, 858*25c28e83SPiotr Jasiukajtis t0z1_l = 1.4102088588676879418739164486159514674310e-37L, 859*25c28e83SPiotr Jasiukajtis #endif 860*25c28e83SPiotr Jasiukajtis /* 0.46163214496836234126265954232572132846819620400644635129599 */ 861*25c28e83SPiotr Jasiukajtis #if defined(__x86) 862*25c28e83SPiotr Jasiukajtis t0z2 = 0.4616321449683623412538115843295472018326L, 863*25c28e83SPiotr Jasiukajtis t0z2_l = 8.84795799617412663558532305039261747030640e-21L, 864*25c28e83SPiotr Jasiukajtis #else 865*25c28e83SPiotr Jasiukajtis t0z2 = 0.46163214496836234126265954232572132343318L, 866*25c28e83SPiotr Jasiukajtis t0z2_l = 5.03501162329616380465302666480916271611101e-36L, 867*25c28e83SPiotr Jasiukajtis #endif 868*25c28e83SPiotr Jasiukajtis /* 0.81977310110050060178786870492160699631174407846245179119586 */ 869*25c28e83SPiotr Jasiukajtis #if defined(__x86) 870*25c28e83SPiotr Jasiukajtis t0z3 = 0.81977310110050060178773362329351925836817L, 871*25c28e83SPiotr Jasiukajtis t0z3_l = 1.350816280877379435658077052534574556256230e-22L 872*25c28e83SPiotr Jasiukajtis #else 873*25c28e83SPiotr Jasiukajtis t0z3 = 0.8197731011005006017878687049216069516957449L, 874*25c28e83SPiotr Jasiukajtis t0z3_l = 4.461599916947014419045492615933551648857380e-35L 875*25c28e83SPiotr Jasiukajtis #endif 876*25c28e83SPiotr Jasiukajtis ; 877*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 878*25c28e83SPiotr Jasiukajtis 879*25c28e83SPiotr Jasiukajtis /* 880*25c28e83SPiotr Jasiukajtis * gamma(x+i) for 0 <= x < 1 881*25c28e83SPiotr Jasiukajtis */ 882*25c28e83SPiotr Jasiukajtis static struct LDouble 883*25c28e83SPiotr Jasiukajtis gam_n(int i, long double x) { 884*25c28e83SPiotr Jasiukajtis struct LDouble rr = {0.0L, 0.0L}, yy; 885*25c28e83SPiotr Jasiukajtis long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl; 886*25c28e83SPiotr Jasiukajtis 887*25c28e83SPiotr Jasiukajtis /* compute yy = gamma(x+1) */ 888*25c28e83SPiotr Jasiukajtis if (x > 0.2845L) { 889*25c28e83SPiotr Jasiukajtis if (x > 0.6374L) { 890*25c28e83SPiotr Jasiukajtis r1 = x - t0z3; 891*25c28e83SPiotr Jasiukajtis r2 = CHOPPED((r1 - t0z3_l)); 892*25c28e83SPiotr Jasiukajtis t2 = r1 - r2; 893*25c28e83SPiotr Jasiukajtis yy = GT3(r2, t2 - t0z3_l); 894*25c28e83SPiotr Jasiukajtis } else { 895*25c28e83SPiotr Jasiukajtis r1 = x - t0z2; 896*25c28e83SPiotr Jasiukajtis r2 = CHOPPED((r1 - t0z2_l)); 897*25c28e83SPiotr Jasiukajtis t2 = r1 - r2; 898*25c28e83SPiotr Jasiukajtis yy = GT2(r2, t2 - t0z2_l); 899*25c28e83SPiotr Jasiukajtis } 900*25c28e83SPiotr Jasiukajtis } else { 901*25c28e83SPiotr Jasiukajtis r1 = x - t0z1; 902*25c28e83SPiotr Jasiukajtis r2 = CHOPPED((r1 - t0z1_l)); 903*25c28e83SPiotr Jasiukajtis t2 = r1 - r2; 904*25c28e83SPiotr Jasiukajtis yy = GT1(r2, t2 - t0z1_l); 905*25c28e83SPiotr Jasiukajtis } 906*25c28e83SPiotr Jasiukajtis /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */ 907*25c28e83SPiotr Jasiukajtis switch (i) { 908*25c28e83SPiotr Jasiukajtis case 0: /* yy/x */ 909*25c28e83SPiotr Jasiukajtis r1 = one / x; 910*25c28e83SPiotr Jasiukajtis xh = CHOPPED((x)); /* x is not tiny */ 911*25c28e83SPiotr Jasiukajtis rr.h = CHOPPED(((yy.h + yy.l) * r1)); 912*25c28e83SPiotr Jasiukajtis rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) - 913*25c28e83SPiotr Jasiukajtis r1 * yy.l); 914*25c28e83SPiotr Jasiukajtis break; 915*25c28e83SPiotr Jasiukajtis case 1: /* yy */ 916*25c28e83SPiotr Jasiukajtis rr.h = yy.h; 917*25c28e83SPiotr Jasiukajtis rr.l = yy.l; 918*25c28e83SPiotr Jasiukajtis break; 919*25c28e83SPiotr Jasiukajtis case 2: /* (x+1)*yy */ 920*25c28e83SPiotr Jasiukajtis z = x + one; /* may not be exact */ 921*25c28e83SPiotr Jasiukajtis zh = CHOPPED((z)); 922*25c28e83SPiotr Jasiukajtis rr.h = zh * yy.h; 923*25c28e83SPiotr Jasiukajtis rr.l = z * yy.l + (x - (zh - one)) * yy.h; 924*25c28e83SPiotr Jasiukajtis break; 925*25c28e83SPiotr Jasiukajtis case 3: /* (x+2)*(x+1)*yy */ 926*25c28e83SPiotr Jasiukajtis z1 = x + one; 927*25c28e83SPiotr Jasiukajtis z2 = x + 2.0L; 928*25c28e83SPiotr Jasiukajtis z = z1 * z2; 929*25c28e83SPiotr Jasiukajtis xh = CHOPPED((z)); 930*25c28e83SPiotr Jasiukajtis zh = CHOPPED((z1)); 931*25c28e83SPiotr Jasiukajtis xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one)); 932*25c28e83SPiotr Jasiukajtis 933*25c28e83SPiotr Jasiukajtis rr.h = xh * yy.h; 934*25c28e83SPiotr Jasiukajtis rr.l = z * yy.l + xl * yy.h; 935*25c28e83SPiotr Jasiukajtis break; 936*25c28e83SPiotr Jasiukajtis 937*25c28e83SPiotr Jasiukajtis case 4: /* (x+1)*(x+3)*(x+2)*yy */ 938*25c28e83SPiotr Jasiukajtis z1 = x + 2.0L; 939*25c28e83SPiotr Jasiukajtis z2 = (x + one) * (x + 3.0L); 940*25c28e83SPiotr Jasiukajtis zh = CHOPPED(z1); 941*25c28e83SPiotr Jasiukajtis zl = x - (zh - 2.0L); 942*25c28e83SPiotr Jasiukajtis xh = CHOPPED(z2); 943*25c28e83SPiotr Jasiukajtis xl = zl * (zh + z1) - (xh - (zh * zh - one)); 944*25c28e83SPiotr Jasiukajtis 945*25c28e83SPiotr Jasiukajtis /* wh+wl=(x+2)*yy */ 946*25c28e83SPiotr Jasiukajtis wh = CHOPPED((z1 * (yy.h + yy.l))); 947*25c28e83SPiotr Jasiukajtis wl = (zl * yy.h + z1 * yy.l) - (wh - zh * yy.h); 948*25c28e83SPiotr Jasiukajtis 949*25c28e83SPiotr Jasiukajtis rr.h = xh * wh; 950*25c28e83SPiotr Jasiukajtis rr.l = z2 * wl + xl * wh; 951*25c28e83SPiotr Jasiukajtis 952*25c28e83SPiotr Jasiukajtis break; 953*25c28e83SPiotr Jasiukajtis case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */ 954*25c28e83SPiotr Jasiukajtis z1 = x + 2.0L; 955*25c28e83SPiotr Jasiukajtis z2 = x + 3.0L; 956*25c28e83SPiotr Jasiukajtis z = z1 * z2; 957*25c28e83SPiotr Jasiukajtis zh = CHOPPED((z1)); 958*25c28e83SPiotr Jasiukajtis yh = CHOPPED((z)); 959*25c28e83SPiotr Jasiukajtis yl = (x - (zh - 2.0L)) * (z2 + zh) - (yh - zh * (zh + one)); 960*25c28e83SPiotr Jasiukajtis z2 = z - 2.0L; 961*25c28e83SPiotr Jasiukajtis z *= z2; 962*25c28e83SPiotr Jasiukajtis xh = CHOPPED((z)); 963*25c28e83SPiotr Jasiukajtis xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L)); 964*25c28e83SPiotr Jasiukajtis rr.h = xh * yy.h; 965*25c28e83SPiotr Jasiukajtis rr.l = z * yy.l + xl * yy.h; 966*25c28e83SPiotr Jasiukajtis break; 967*25c28e83SPiotr Jasiukajtis case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */ 968*25c28e83SPiotr Jasiukajtis z1 = x + 2.0L; 969*25c28e83SPiotr Jasiukajtis z2 = x + 3.0L; 970*25c28e83SPiotr Jasiukajtis z = z1 * z2; 971*25c28e83SPiotr Jasiukajtis zh = CHOPPED((z1)); 972*25c28e83SPiotr Jasiukajtis yh = CHOPPED((z)); 973*25c28e83SPiotr Jasiukajtis z1 = x - (zh - 2.0L); 974*25c28e83SPiotr Jasiukajtis yl = z1 * (z2 + zh) - (yh - zh * (zh + one)); 975*25c28e83SPiotr Jasiukajtis z2 = z - 2.0L; 976*25c28e83SPiotr Jasiukajtis x5 = x + 5.0L; 977*25c28e83SPiotr Jasiukajtis z *= z2; 978*25c28e83SPiotr Jasiukajtis xh = CHOPPED(z); 979*25c28e83SPiotr Jasiukajtis zh += 3.0; 980*25c28e83SPiotr Jasiukajtis xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L)); 981*25c28e83SPiotr Jasiukajtis /* xh+xl=(x+1)*...*(x+4) */ 982*25c28e83SPiotr Jasiukajtis /* wh+wl=(x+5)*yy */ 983*25c28e83SPiotr Jasiukajtis wh = CHOPPED((x5 * (yy.h + yy.l))); 984*25c28e83SPiotr Jasiukajtis wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h); 985*25c28e83SPiotr Jasiukajtis rr.h = wh * xh; 986*25c28e83SPiotr Jasiukajtis rr.l = z * wl + xl * wh; 987*25c28e83SPiotr Jasiukajtis break; 988*25c28e83SPiotr Jasiukajtis case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */ 989*25c28e83SPiotr Jasiukajtis z1 = x + 3.0L; 990*25c28e83SPiotr Jasiukajtis z2 = x + 4.0L; 991*25c28e83SPiotr Jasiukajtis z = z2 * z1; 992*25c28e83SPiotr Jasiukajtis zh = CHOPPED((z1)); 993*25c28e83SPiotr Jasiukajtis yh = CHOPPED((z)); /* yh+yl = (x+3)(x+4) */ 994*25c28e83SPiotr Jasiukajtis yl = (x - (zh - 3.0L)) * (z2 + zh) - (yh - (zh * (zh + one))); 995*25c28e83SPiotr Jasiukajtis z1 = x + 6.0L; 996*25c28e83SPiotr Jasiukajtis z2 = z - 2.0L; /* z2 = (x+2)*(x+5) */ 997*25c28e83SPiotr Jasiukajtis z *= z2; 998*25c28e83SPiotr Jasiukajtis xh = CHOPPED((z)); 999*25c28e83SPiotr Jasiukajtis xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L)); 1000*25c28e83SPiotr Jasiukajtis /* xh+xl=(x+2)*...*(x+5) */ 1001*25c28e83SPiotr Jasiukajtis /* wh+wl=(x+1)(x+6)*yy */ 1002*25c28e83SPiotr Jasiukajtis z2 -= 4.0L; /* z2 = (x+1)(x+6) */ 1003*25c28e83SPiotr Jasiukajtis wh = CHOPPED((z2 * (yy.h + yy.l))); 1004*25c28e83SPiotr Jasiukajtis wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0L) * yy.h); 1005*25c28e83SPiotr Jasiukajtis rr.h = wh * xh; 1006*25c28e83SPiotr Jasiukajtis rr.l = z * wl + xl * wh; 1007*25c28e83SPiotr Jasiukajtis } 1008*25c28e83SPiotr Jasiukajtis return (rr); 1009*25c28e83SPiotr Jasiukajtis } 1010*25c28e83SPiotr Jasiukajtis 1011*25c28e83SPiotr Jasiukajtis long double 1012*25c28e83SPiotr Jasiukajtis tgammal(long double x) { 1013*25c28e83SPiotr Jasiukajtis struct LDouble ss, ww; 1014*25c28e83SPiotr Jasiukajtis long double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5; 1015*25c28e83SPiotr Jasiukajtis int i, j, m, ix, hx, xk; 1016*25c28e83SPiotr Jasiukajtis unsigned lx; 1017*25c28e83SPiotr Jasiukajtis 1018*25c28e83SPiotr Jasiukajtis hx = H0_WORD(x); 1019*25c28e83SPiotr Jasiukajtis lx = H3_WORD(x); 1020*25c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff; 1021*25c28e83SPiotr Jasiukajtis y = x; 1022*25c28e83SPiotr Jasiukajtis if (ix < 0x3f8e0000) { /* x < 2**-113 */ 1023*25c28e83SPiotr Jasiukajtis return (one / x); 1024*25c28e83SPiotr Jasiukajtis } 1025*25c28e83SPiotr Jasiukajtis if (ix >= 0x7fff0000) 1026*25c28e83SPiotr Jasiukajtis return (x * ((hx < 0)? zero : x)); /* Inf or NaN */ 1027*25c28e83SPiotr Jasiukajtis if (x > overflow) /* overflow threshold */ 1028*25c28e83SPiotr Jasiukajtis return (x * 1.0e4932L); 1029*25c28e83SPiotr Jasiukajtis if (hx >= 0x40020000) { /* x >= 8 */ 1030*25c28e83SPiotr Jasiukajtis ww = large_gam(x, &m); 1031*25c28e83SPiotr Jasiukajtis w = ww.h + ww.l; 1032*25c28e83SPiotr Jasiukajtis return (scalbnl(w, m)); 1033*25c28e83SPiotr Jasiukajtis } 1034*25c28e83SPiotr Jasiukajtis 1035*25c28e83SPiotr Jasiukajtis if (hx > 0) { /* 0 < x < 8 */ 1036*25c28e83SPiotr Jasiukajtis i = (int) x; 1037*25c28e83SPiotr Jasiukajtis ww = gam_n(i, x - (long double) i); 1038*25c28e83SPiotr Jasiukajtis return (ww.h + ww.l); 1039*25c28e83SPiotr Jasiukajtis } 1040*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 1041*25c28e83SPiotr Jasiukajtis /* negative x */ 1042*25c28e83SPiotr Jasiukajtis /* 1043*25c28e83SPiotr Jasiukajtis * compute xk = 1044*25c28e83SPiotr Jasiukajtis * -2 ... x is an even int (-inf is considered an even #) 1045*25c28e83SPiotr Jasiukajtis * -1 ... x is an odd int 1046*25c28e83SPiotr Jasiukajtis * +0 ... x is not an int but chopped to an even int 1047*25c28e83SPiotr Jasiukajtis * +1 ... x is not an int but chopped to an odd int 1048*25c28e83SPiotr Jasiukajtis */ 1049*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 1050*25c28e83SPiotr Jasiukajtis xk = 0; 1051*25c28e83SPiotr Jasiukajtis #if defined(__x86) 1052*25c28e83SPiotr Jasiukajtis if (ix >= 0x403e0000) { /* x >= 2**63 } */ 1053*25c28e83SPiotr Jasiukajtis if (ix >= 0x403f0000) 1054*25c28e83SPiotr Jasiukajtis xk = -2; 1055*25c28e83SPiotr Jasiukajtis else 1056*25c28e83SPiotr Jasiukajtis xk = -2 + (lx & 1); 1057*25c28e83SPiotr Jasiukajtis #else 1058*25c28e83SPiotr Jasiukajtis if (ix >= 0x406f0000) { /* x >= 2**112 */ 1059*25c28e83SPiotr Jasiukajtis if (ix >= 0x40700000) 1060*25c28e83SPiotr Jasiukajtis xk = -2; 1061*25c28e83SPiotr Jasiukajtis else 1062*25c28e83SPiotr Jasiukajtis xk = -2 + (lx & 1); 1063*25c28e83SPiotr Jasiukajtis #endif 1064*25c28e83SPiotr Jasiukajtis } else if (ix >= 0x3fff0000) { 1065*25c28e83SPiotr Jasiukajtis w = -x; 1066*25c28e83SPiotr Jasiukajtis t1 = floorl(w); 1067*25c28e83SPiotr Jasiukajtis t2 = t1 * half; 1068*25c28e83SPiotr Jasiukajtis t3 = floorl(t2); 1069*25c28e83SPiotr Jasiukajtis if (t1 == w) { 1070*25c28e83SPiotr Jasiukajtis if (t2 == t3) 1071*25c28e83SPiotr Jasiukajtis xk = -2; 1072*25c28e83SPiotr Jasiukajtis else 1073*25c28e83SPiotr Jasiukajtis xk = -1; 1074*25c28e83SPiotr Jasiukajtis } else { 1075*25c28e83SPiotr Jasiukajtis if (t2 == t3) 1076*25c28e83SPiotr Jasiukajtis xk = 0; 1077*25c28e83SPiotr Jasiukajtis else 1078*25c28e83SPiotr Jasiukajtis xk = 1; 1079*25c28e83SPiotr Jasiukajtis } 1080*25c28e83SPiotr Jasiukajtis } 1081*25c28e83SPiotr Jasiukajtis 1082*25c28e83SPiotr Jasiukajtis if (xk < 0) { 1083*25c28e83SPiotr Jasiukajtis /* return NaN. Ideally gamma(-n)= (-1)**(n+1) * inf */ 1084*25c28e83SPiotr Jasiukajtis return (x - x) / (x - x); 1085*25c28e83SPiotr Jasiukajtis } 1086*25c28e83SPiotr Jasiukajtis 1087*25c28e83SPiotr Jasiukajtis /* 1088*25c28e83SPiotr Jasiukajtis * negative underflow thresold -(1774+9ulp) 1089*25c28e83SPiotr Jasiukajtis */ 1090*25c28e83SPiotr Jasiukajtis if (x < -1774.0000000000000000000000000000017749370L) { 1091*25c28e83SPiotr Jasiukajtis z = tiny / x; 1092*25c28e83SPiotr Jasiukajtis if (xk == 1) 1093*25c28e83SPiotr Jasiukajtis z = -z; 1094*25c28e83SPiotr Jasiukajtis return (z * tiny); 1095*25c28e83SPiotr Jasiukajtis } 1096*25c28e83SPiotr Jasiukajtis 1097*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 1098*25c28e83SPiotr Jasiukajtis /* 1099*25c28e83SPiotr Jasiukajtis * now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x 1100*25c28e83SPiotr Jasiukajtis */ 1101*25c28e83SPiotr Jasiukajtis /* 1102*25c28e83SPiotr Jasiukajtis * First compute ss = -sin(pi*y)/pi so that 1103*25c28e83SPiotr Jasiukajtis * gamma(x) = 1/(ss*gamma(1+y)) 1104*25c28e83SPiotr Jasiukajtis */ 1105*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 1106*25c28e83SPiotr Jasiukajtis y = -x; 1107*25c28e83SPiotr Jasiukajtis j = (int) y; 1108*25c28e83SPiotr Jasiukajtis z = y - (long double) j; 1109*25c28e83SPiotr Jasiukajtis if (z > 0.3183098861837906715377675L) 1110*25c28e83SPiotr Jasiukajtis if (z > 0.6816901138162093284622325L) 1111*25c28e83SPiotr Jasiukajtis ss = kpsin(one - z); 1112*25c28e83SPiotr Jasiukajtis else 1113*25c28e83SPiotr Jasiukajtis ss = kpcos(0.5L - z); 1114*25c28e83SPiotr Jasiukajtis else 1115*25c28e83SPiotr Jasiukajtis ss = kpsin(z); 1116*25c28e83SPiotr Jasiukajtis if (xk == 0) { 1117*25c28e83SPiotr Jasiukajtis ss.h = -ss.h; 1118*25c28e83SPiotr Jasiukajtis ss.l = -ss.l; 1119*25c28e83SPiotr Jasiukajtis } 1120*25c28e83SPiotr Jasiukajtis 1121*25c28e83SPiotr Jasiukajtis /* Then compute ww = gamma(1+y), note that result scale to 2**m */ 1122*25c28e83SPiotr Jasiukajtis m = 0; 1123*25c28e83SPiotr Jasiukajtis if (j < 7) { 1124*25c28e83SPiotr Jasiukajtis ww = gam_n(j + 1, z); 1125*25c28e83SPiotr Jasiukajtis } else { 1126*25c28e83SPiotr Jasiukajtis w = y + one; 1127*25c28e83SPiotr Jasiukajtis if ((lx & 1) == 0) { /* y+1 exact (note that y<184) */ 1128*25c28e83SPiotr Jasiukajtis ww = large_gam(w, &m); 1129*25c28e83SPiotr Jasiukajtis } else { 1130*25c28e83SPiotr Jasiukajtis t = w - one; 1131*25c28e83SPiotr Jasiukajtis if (t == y) { /* y+one exact */ 1132*25c28e83SPiotr Jasiukajtis ww = large_gam(w, &m); 1133*25c28e83SPiotr Jasiukajtis } else { /* use y*gamma(y) */ 1134*25c28e83SPiotr Jasiukajtis if (j == 7) 1135*25c28e83SPiotr Jasiukajtis ww = gam_n(j, z); 1136*25c28e83SPiotr Jasiukajtis else 1137*25c28e83SPiotr Jasiukajtis ww = large_gam(y, &m); 1138*25c28e83SPiotr Jasiukajtis t4 = ww.h + ww.l; 1139*25c28e83SPiotr Jasiukajtis t1 = CHOPPED((y)); 1140*25c28e83SPiotr Jasiukajtis t2 = CHOPPED((t4)); 1141*25c28e83SPiotr Jasiukajtis /* t4 will not be too large */ 1142*25c28e83SPiotr Jasiukajtis ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2; 1143*25c28e83SPiotr Jasiukajtis ww.h = t1 * t2; 1144*25c28e83SPiotr Jasiukajtis } 1145*25c28e83SPiotr Jasiukajtis } 1146*25c28e83SPiotr Jasiukajtis } 1147*25c28e83SPiotr Jasiukajtis 1148*25c28e83SPiotr Jasiukajtis /* compute 1/(ss*ww) */ 1149*25c28e83SPiotr Jasiukajtis t3 = ss.h + ss.l; 1150*25c28e83SPiotr Jasiukajtis t4 = ww.h + ww.l; 1151*25c28e83SPiotr Jasiukajtis t1 = CHOPPED((t3)); 1152*25c28e83SPiotr Jasiukajtis t2 = CHOPPED((t4)); 1153*25c28e83SPiotr Jasiukajtis z1 = ss.l - (t1 - ss.h); /* (t1,z1) = ss */ 1154*25c28e83SPiotr Jasiukajtis z2 = ww.l - (t2 - ww.h); /* (t2,z2) = ww */ 1155*25c28e83SPiotr Jasiukajtis t3 = t3 * t4; /* t3 = ss*ww */ 1156*25c28e83SPiotr Jasiukajtis z3 = one / t3; /* z3 = 1/(ss*ww) */ 1157*25c28e83SPiotr Jasiukajtis t5 = t1 * t2; 1158*25c28e83SPiotr Jasiukajtis z5 = z1 * t4 + t1 * z2; /* (t5,z5) = ss*ww */ 1159*25c28e83SPiotr Jasiukajtis t1 = CHOPPED((t3)); /* (t1,z1) = ss*ww */ 1160*25c28e83SPiotr Jasiukajtis z1 = z5 - (t1 - t5); 1161*25c28e83SPiotr Jasiukajtis t2 = CHOPPED((z3)); /* leading 1/(ss*ww) */ 1162*25c28e83SPiotr Jasiukajtis z2 = z3 * (t2 * z1 - (one - t2 * t1)); 1163*25c28e83SPiotr Jasiukajtis z = t2 - z2; 1164*25c28e83SPiotr Jasiukajtis 1165*25c28e83SPiotr Jasiukajtis return (scalbnl(z, -m)); 1166*25c28e83SPiotr Jasiukajtis } 1167