Print this page
11210 libm should be cstyle(1ONBLD) clean


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */

  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak __tgammal = tgammal
  31 
  32 #include "libm.h"
  33 #include <sys/isa_defs.h>
  34 
  35 #if defined(_BIG_ENDIAN)
  36 #define H0_WORD(x)      ((unsigned *) &x)[0]
  37 #define H3_WORD(x)      ((unsigned *) &x)[3]
  38 #define CHOPPED(x)      (long double) ((double) (x))
  39 #else
  40 #define H0_WORD(x)      ((((int *) &x)[2] << 16) | \
  41                         (0x0000ffff & (((unsigned *) &x)[1] >> 15)))
  42 #define H3_WORD(x)      ((unsigned *) &x)[0]
  43 #define CHOPPED(x)      (long double) ((float) (x))
  44 #endif
  45 
  46 struct LDouble {
  47         long double h, l;
  48 };
  49 
  50 /* INDENT OFF */
  51 /* Primary interval GTi() */

  52 static const long double P1[] = {
  53         +0.709086836199777919037185741507610124611513720557L,
  54         +4.45754781206489035827915969367354835667391606951e-0001L,
  55         +3.21049298735832382311662273882632210062918153852e-0002L,
  56         -5.71296796342106617651765245858289197369688864350e-0003L,
  57         +6.04666892891998977081619174969855831606965352773e-0003L,
  58         +8.99106186996888711939627812174765258822658645168e-0004L,
  59         -6.96496846144407741431207008527018441810175568949e-0005L,
  60         +1.52597046118984020814225409300131445070213882429e-0005L,
  61         +5.68521076168495673844711465407432189190681541547e-0007L,
  62         +3.30749673519634895220582062520286565610418952979e-0008L,
  63 };

  64 static const long double Q1[] = {
  65         +1.0+0000L,
  66         +1.35806511721671070408570853537257079579490650668e+0000L,
  67         +2.97567810153429553405327140096063086994072952961e-0001L,
  68         -1.52956835982588571502954372821681851681118097870e-0001L,
  69         -2.88248519561420109768781615289082053597954521218e-0002L,
  70         +1.03475311719937405219789948456313936302378395955e-0002L,
  71         +4.12310203243891222368965360124391297374822742313e-0004L,
  72         -3.12653708152290867248931925120380729518332507388e-0004L,
  73         +2.36672170850409745237358105667757760527014332458e-0005L,
  74 };

  75 static const long double P2[] = {
  76         +0.428486815855585429730209907810650135255270600668084114L,
  77         +2.62768479103809762805691743305424077975230551176e-0001L,
  78         +3.81187532685392297608310837995193946591425896150e-0002L,
  79         +3.00063075891811043820666846129131255948527925381e-0003L,
  80         +2.47315407812279164228398470797498649142513408654e-0003L,
  81         +3.62838199917848372586173483147214880464782938664e-0004L,
  82         +3.43991105975492623982725644046473030098172692423e-0006L,
  83         +4.56902151569603272237014240794257659159045432895e-0006L,
  84         +2.13734755837595695602045100675540011352948958453e-0007L,
  85         +9.74123440547918230781670266967882492234877125358e-0009L,
  86 };

  87 static const long double Q2[] = {
  88         +1.0L,
  89         +9.18284118632506842664645516830761489700556179701e-0001L,
  90         -6.41430858837830766045202076965923776189154874947e-0003L,
  91         -1.24400885809771073213345747437964149775410921376e-0001L,
  92         +4.69803798146251757538856567522481979624746875964e-0003L,
  93         +7.18309447069495315914284705109868696262662082731e-0003L,
  94         -8.75812626987894695112722600697653425786166399105e-0004L,
  95         -1.23539972377769277995959339188431498626674835169e-0004L,
  96         +3.10019017590151598732360097849672925448587547746e-0005L,
  97         -1.77260223349332617658921874288026777465782364070e-0006L,
  98 };

  99 static const long double P3[] = {
 100         +0.3824094797345675048502747661075355640070439388902L,
 101         +3.42198093076618495415854906335908427159833377774e-0001L,
 102         +9.63828189500585568303961406863153237440702754858e-0002L,
 103         +8.76069421042696384852462044188520252156846768667e-0003L,
 104         +1.86477890389161491224872014149309015261897537488e-0003L,
 105         +8.16871354540309895879974742853701311541286944191e-0004L,
 106         +6.83783483674600322518695090864659381650125625216e-0005L,
 107         -1.10168269719261574708565935172719209272190828456e-0006L,
 108         +9.66243228508380420159234853278906717065629721016e-0007L,
 109         +2.31858885579177250541163820671121664974334728142e-0008L,
 110 };

 111 static const long double Q3[] = {
 112         +1.0L,
 113         +8.25479821168813634632437430090376252512793067339e-0001L,
 114         -1.62251363073937769739639623669295110346015576320e-0002L,
 115         -1.10621286905916732758745130629426559691187579852e-0001L,
 116         +3.48309693970985612644446415789230015515365291459e-0003L,
 117         +6.73553737487488333032431261131289672347043401328e-0003L,
 118         -7.63222008393372630162743587811004613050245128051e-0004L,
 119         -1.35792670669190631476784768961953711773073251336e-0004L,
 120         +3.19610150954223587006220730065608156460205690618e-0005L,
 121         -1.82096553862822346610109522015129585693354348322e-0006L,
 122 };
 123 
 124 static const long double
 125 #if defined(__x86)
 126 GZ1_h   =  0.938204627909682449364570100414084663498215377L,
 127 GZ1_l   =  4.518346116624229420055327632718530617227944106e-20L,
 128 GZ2_h   =  0.885603194410888700264725126309883762587560340L,
 129 GZ2_l   =  1.409077427270497062039119290776508217077297169e-20L,
 130 GZ3_h   =  0.936781411463652321613537060640553022494714241L,
 131 GZ3_l   =  5.309836440284827247897772963887219035221996813e-21L,
 132 #else
 133 GZ1_h   =  0.938204627909682449409753561580326910854647031L,
 134 GZ1_l   =  4.684412162199460089642452580902345976446297037e-35L,
 135 GZ2_h   =  0.885603194410888700278815900582588658192658794L,
 136 GZ2_l   =  7.501529273890253789219935569758713534641074860e-35L,
 137 GZ3_h   =  0.936781411463652321618846897080837818855399840L,
 138 GZ3_l   =  3.088721217404784363585591914529361687403776917e-35L,
 139 #endif
 140 TZ1     = -0.3517214357852935791015625L,
 141 TZ3     =  0.280530631542205810546875L;
 142 /* INDENT ON */
 143 
 144 /* INDENT OFF */
 145 /*
 146  * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
 147  * ...assume yh got 53 or 24(i386) significant bits
 148  */
 149 /* INDENT ON */
 150 static struct LDouble
 151 GT1(long double yh, long double yl) {

 152         long double t3, t4, y;
 153         int i;
 154         struct LDouble r;
 155 
 156         y = yh + yl;

 157         for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) {
 158                 t4 = t4 * y + Q1[i];
 159                 t3 = t3 * y + P1[i];
 160         }

 161         t3 = (y * y) * t3 / t4;
 162         t3 += (TZ1 * yl + GZ1_l);
 163         t4 = TZ1 * yh;
 164         r.h = CHOPPED((t4 + GZ1_h + t3));
 165         t3 += (t4 - (r.h - GZ1_h));
 166         r.l = t3;
 167         return (r);
 168 }
 169 
 170 /* INDENT OFF */
 171 /*
 172  * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
 173  * ...assume yh got 53 significant bits
 174  */
 175 /* INDENT ON */
 176 static struct LDouble
 177 GT2(long double yh, long double yl) {

 178         long double t3, t4, y;
 179         int i;
 180         struct LDouble r;
 181 
 182         y = yh + yl;

 183         for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) {
 184                 t4 = t4 * y + Q2[i];
 185                 t3 = t3 * y + P2[i];
 186         }

 187         t3 = GZ2_l + (y * y) * t3 / t4;
 188         r.h = CHOPPED((GZ2_h + t3));
 189         r.l = t3 - (r.h - GZ2_h);
 190         return (r);
 191 }
 192 
 193 /* INDENT OFF */
 194 /*
 195  * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
 196  * ...assume yh got 53 significant bits
 197  */
 198 /* INDENT ON */
 199 static struct LDouble
 200 GT3(long double yh, long double yl) {

 201         long double t3, t4, y;
 202         int i;
 203         struct LDouble r;
 204 
 205         y = yh + yl;

 206         for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) {
 207                 t4 = t4 * y + Q3[i];
 208                 t3 = t3 * y + P3[i];
 209         }

 210         t3 = (y * y) * t3 / t4;
 211         t3 += (TZ3 * yl + GZ3_l);
 212         t4 = TZ3 * yh;
 213         r.h = CHOPPED((t4 + GZ3_h + t3));
 214         t3 += (t4 - (r.h - GZ3_h));
 215         r.l = t3;
 216         return (r);
 217 }
 218 
 219 /* INDENT OFF */
 220 /* Hex value of GP[0] shoule be 3FB55555 55555555 */

 221 static const long double GP[] = {
 222         +0.083333333333333333333333333333333172839171301L,
 223         -2.77777777777777777777777777492501211999399424104e-0003L,
 224         +7.93650793650793650793635650541638236350020883243e-0004L,
 225         -5.95238095238095238057299772679324503339241961704e-0004L,
 226         +8.41750841750841696138422987977683524926142600321e-0004L,
 227         -1.91752691752686682825032547823699662178842123308e-0003L,
 228         +6.41025641022403480921891559356473451161279359322e-0003L,
 229         -2.95506535798414019189819587455577003732808185071e-0002L,
 230         +1.79644367229970031486079180060923073476568732136e-0001L,
 231         -1.39243086487274662174562872567057200255649290646e+0000L,
 232         +1.34025874044417962188677816477842265259608269775e+0001L,
 233         -1.56803713480127469414495545399982508700748274318e+0002L,
 234         +2.18739841656201561694927630335099313968924493891e+0003L,
 235         -3.55249848644100338419187038090925410976237921269e+0004L,
 236         +6.43464880437835286216768959439484376449179576452e+0005L,
 237         -1.20459154385577014992600342782821389605893904624e+0007L,
 238         +2.09263249637351298563934942349749718491071093210e+0008L,
 239         -2.96247483183169219343745316433899599834685703457e+0009L,
 240         +2.88984933605896033154727626086506756972327292981e+0010L,
 241         -1.40960434146030007732838382416230610302678063984e+0011L,      /* 19 */
 242 };
 243 
 244 static const long double T3[] = {
 245         +0.666666666666666666666666666666666634567834260213L,   /* T3[0] */
 246         +0.400000000000000000000000000040853636176634934140L,   /* T3[1] */
 247         +0.285714285714285714285696975252753987869020263448L,   /* T3[2] */
 248         +0.222222222222222225593221101192317258554772129875L,   /* T3[3] */
 249         +0.181818181817850192105847183461778186703779262916L,   /* T3[4] */
 250         +0.153846169861348633757101285952333369222567014596L,   /* T3[5] */
 251         +0.133033462889260193922261296772841229985047571265L,   /* T3[6] */
 252 };
 253 
 254 static const long double c[] = {
 255 0.0L,
 256 1.0L,
 257 2.0L,
 258 0.5L,
 259 1.0e-4930L,                                                     /* tiny */
 260 4.18937683105468750000e-01L,                                    /* hln2pim1_h */
 261 8.50099203991780329736405617639861397473637783412817152e-07L,   /* hln2pim1_l */
 262 0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
 263 2.16608493865351192653179168701171875e-02L,                     /* ln2_32hi */
 264 5.96317165397058692545083025235937919875797669127130e-12L,      /* ln2_32lo */
 265 46.16624130844682903551758979206054839765267053289554989233L,   /* invln2_32 */
 266 #if defined(__x86)
 267 1.7555483429044629170023839037639845628291e+03L,                /* overflow */
 268 #else
 269 1.7555483429044629170038892160702032034177e+03L,                /* overflow */
 270 #endif
 271 };

 272 
 273 #define zero            c[0]
 274 #define one             c[1]
 275 #define two             c[2]
 276 #define half            c[3]
 277 #define tiny            c[4]
 278 #define hln2pim1_h      c[5]
 279 #define hln2pim1_l      c[6]
 280 #define hln2pim1        c[7]
 281 #define ln2_32hi        c[8]
 282 #define ln2_32lo        c[9]
 283 #define invln2_32       c[10]
 284 #define overflow        c[11]
 285 
 286 /*
 287  * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64
 288  */
 289 static const long double Et[] = {
 290         +5.0000000000000000000e-1L,
 291         +1.66666666666666666666666666666828835166292152466e-0001L,


 297         +2.75573192239028921114572986441972140933432317798e-0006L,
 298         +2.75573192239448470555548102895526369739856219317e-0007L,
 299         +2.50521677867683935940853997995937600214167232477e-0008L,
 300         +2.08767928899010367374984448513685566514152147362e-0009L,
 301 };
 302 
 303 /*
 304  * long double precision coefficients for computing log(x)-1 in tgamma.
 305  *  See "algorithm" for details
 306  *
 307  *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 308  *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 309  *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 310  *       T2(j) = T2[2j,2j+1] = log(z[j]),
 311  *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15
 312  *  Note
 313  *  (1) the leading entries are truncated to 24 binary point.
 314  *  (2) Remez error for T3(s) is bounded by 2**(-136.54)
 315  */
 316 static const long double T1[] = {
 317 -1.000000000000000000000000000000000000000000e+00L,
 318         +0.000000000000000000000000000000000000000000e+00L,
 319 -3.068528175354003906250000000000000000000000e-01L,
 320 -1.904654299957767878541823431924500011926579e-09L,
 321         +3.862943053245544433593750000000000000000000e-01L,
 322         +5.579533617547508924291635313615100141107647e-08L,
 323         +1.079441487789154052734375000000000000000000e+00L,
 324         +5.389068187551732136437452970422650211661470e-08L,
 325         +1.772588670253753662109375000000000000000000e+00L,
 326         +5.198602757555955348583270627230200282215294e-08L,
 327         +2.465735852718353271484375000000000000000000e+00L,
 328         +5.008137327560178560729088284037750352769117e-08L,
 329         +3.158883035182952880859375000000000000000000e+00L,
 330         +4.817671897564401772874905940845299849351090e-08L,
 331         +3.852030217647552490234375000000000000000000e+00L,
 332         +4.627206467568624985020723597652849919904913e-08L,
 333         +4.545177400112152099609375000000000000000000e+00L,
 334         +4.436741037572848197166541254460399990458737e-08L,
 335         +5.238324582576751708984375000000000000000000e+00L,
 336         +4.246275607577071409312358911267950061012560e-08L,
 337         +5.931471765041351318359375000000000000000000e+00L,
 338         +4.055810177581294621458176568075500131566384e-08L,
 339 };
 340 


 527         +1.35425554693689272829801474014070273e+00L,
 528         +1.38390988196383195487265952726519287e+00L,
 529         +1.41421356237309504880168872420969798e+00L,
 530         +1.44518080697704662003700624147167095e+00L,
 531         +1.47682614593949931138690748037404985e+00L,
 532         +1.50916442759342273976601955103319352e+00L,
 533         +1.54221082540794082361229186209073479e+00L,
 534         +1.57598084510788648645527016018190504e+00L,
 535         +1.61049033194925430817952066735740067e+00L,
 536         +1.64575547815396484451875672472582254e+00L,
 537         +1.68179283050742908606225095246642969e+00L,
 538         +1.71861929812247791562934437645631244e+00L,
 539         +1.75625216037329948311216061937531314e+00L,
 540         +1.79470907500310718642770324212778174e+00L,
 541         +1.83400808640934246348708318958828892e+00L,
 542         +1.87416763411029990132999894995444645e+00L,
 543         +1.91520656139714729387261127029583086e+00L,
 544         +1.95714412417540026901832225162687149e+00L,
 545 #endif
 546 };

 547 static const long double S_trail[] = {
 548 #if defined(__x86)
 549         +0.0000000000000000000000000e+00L,
 550         +2.6327965667180882569382524e-20L,
 551         +8.3765863521895191129661899e-20L,
 552         +3.9798705777454504249209575e-20L,
 553         +1.0668046596651558640993042e-19L,
 554         +1.9376009847285360448117114e-20L,
 555         +6.7081819456112953751277576e-21L,
 556         +1.9711680502629186462729727e-20L,
 557         +2.9932584438449523689104569e-20L,
 558         +6.8887754153039109411061914e-20L,
 559         +6.8002718741225378942847820e-20L,
 560         +6.5846917376975403439742349e-20L,
 561         +1.2171958727511372194876001e-20L,
 562         +3.5625253228704087115438260e-20L,
 563         +3.1129551559077560956309179e-20L,
 564         +5.7519192396164779846216492e-20L,
 565         +3.7900651177865141593101239e-20L,
 566         +1.1659262405698741798080115e-20L,
 567         +7.1364385105284695967172478e-20L,
 568         +5.2631003710812203588788949e-20L,
 569         +2.6328853788732632868460580e-20L,
 570         +5.4583950085438242788190141e-20L,
 571         +9.5803254376938269960718656e-20L,
 572         +7.6837733983874245823512279e-21L,
 573         +2.4415965910835093824202087e-20L,
 574         +2.6052966871016580981769728e-20L,
 575         +2.6876456344632553875309579e-21L,
 576         +1.2861930155613700201703279e-20L,
 577         +8.8166633394037485606572294e-20L,
 578         +2.9788615389580190940837037e-20L,
 579         +5.2352341619805098677422139e-20L,
 580         +5.2578463064010463732242363e-20L,
 581 #else
 582         +0.00000000000000000000000000000000000e+00L,
 583         +1.80506787420330954745573333054573786e-35L,
 584 -9.37452029228042742195756741973083214e-35L,
 585 -1.59696844729275877071290963023149997e-35L,
 586         +9.11249341012502297851168610167248666e-35L,
 587 -6.50422820697854828723037477525938871e-35L,
 588 -8.14846884452585113732569176748815532e-35L,
 589 -5.06621457672180031337233074514290335e-35L,
 590 -1.35983097468881697374987563824591912e-35L,
 591         +9.49742763556319647030771056643324660e-35L,
 592 -3.28317052317699860161506596533391526e-36L,
 593 -5.01723570938719041029018653045842895e-35L,
 594 -2.39147479768910917162283430160264014e-35L,
 595 -8.35057135763390881529889073794408385e-36L,
 596         +7.03675688907326504242173719067187644e-35L,
 597 -5.18248485306464645753689301856695619e-35L,
 598         +9.42224254862183206569211673639406488e-35L,
 599 -3.96750082539886230916730613021641828e-35L,
 600         +7.14352899156330061452327361509276724e-35L,
 601         +1.15987125286798512424651783410044433e-35L,
 602         +4.69693347835811549530973921320187447e-35L,
 603 -3.38651317599500471079924198499981917e-35L,
 604 -8.58731877429824706886865593510387445e-35L,
 605 -9.60595154874935050318549936224606909e-35L,
 606         +9.60973393212801278450755869714178581e-35L,
 607         +6.37839792144002843924476144978084855e-35L,
 608         +7.79243078569586424945646112516927770e-35L,
 609         +7.36133776758845652413193083663393220e-35L,
 610 -6.47299514791334723003521457561217053e-35L,
 611         +8.58747441795369869427879806229522962e-35L,
 612         +2.37181542282517483569165122830269098e-35L,
 613 -3.02689168209611877300459737342190031e-37L,
 614 #endif
 615 };
 616 /* INDENT ON */
 617 
 618 /* INDENT OFF */
 619 /*
 620  * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
 621  *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
 622  *                = L1 + L2 + L3,
 623  */
 624 /* INDENT ON */
 625 static struct LDouble
 626 large_gam(long double x, int *m) {

 627         long double z, t1, t2, t3, z2, t5, w, y, u, r, v;
 628         long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L;
 629         int n2, j2, k, ix, j, i;
 630         struct LDouble zz;
 631         long double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
 632 
 633 /* INDENT OFF */
 634 /*
 635  * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
 636  *
 637  *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 638  *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 639  *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 640  *       T2(j) = T2[2j,2j+1] = log(z[j]),
 641  *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15
 642  *  Note
 643  *  (1) the leading entries are truncated to 24 binary point.
 644  *  (2) Remez error for T3(s) is bounded by 2**(-72.4)
 645  *                                   2**(-24)
 646  *                           _________V___________________
 647  *               T1(n):     |_________|___________________|
 648  *                             _______ ______________________
 649  *               T2(j):       |_______|______________________|
 650  *                                ____ _______________________
 651  *               2s:             |____|_______________________|
 652  *                                    __________________________
 653  *          +    T3(s)-2s:           |__________________________|
 654  *                       -------------------------------------------
 655  *                          [leading] + [Trailing]
 656  */
 657         /* INDENT ON */
 658         ix = H0_WORD(x);
 659         n2 = (ix >> 16) - 0x3fff; /* exponent of x, range:3-10 */
 660         y = scalbnl(x, -n2);    /* y = scale x to [1,2] */
 661         n2 += n2;               /* 2n */
 662         j = (ix >> 10) & 0x3f;        /* j */
 663         z = 1.0078125L + (long double) j * 0.015625L;   /* z[j]=1+j/64+1/128 */
 664         j2 = j + j;
 665         t1 = y + z;
 666         t2 = y - z;
 667         r = one / t1;
 668         u = r * t2;             /* u = (y-z)/(y+z) */
 669         t1 = CHOPPED(t1);
 670         t4 = T2[j2 + 1] + T1[n2 + 1];
 671         z2 = u * u;
 672         k = H0_WORD(u) & 0x7fffffff;
 673         t3 = T2[j2] + T1[n2];

 674         for (t5 = T3[6], i = 5; i >= 0; i--)
 675                 t5 = z2 * t5 + T3[i];

 676         if ((k >> 16) < 0x3fec) {      /* |u|<2**-19 */
 677                 t2 = t4 + u * (two + z2 * t5);
 678         } else {
 679                 t5 = t4 + (u * z2) * t5;
 680                 u2 = u + u;
 681                 v = (long double) ((int) (u2 * t24)) * p24;
 682                 t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
 683                 t3 += v;
 684         }

 685         ss_h = CHOPPED((t2 + t3));
 686         ss_l = t2 - (ss_h - t3);
 687 /* INDENT OFF */
 688 /*
 689  * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
 690  * where ss = log(x) - 1 in already in extra precision
 691  */
 692         /* INDENT ON */
 693         z = one / x;
 694         r = x - half;
 695         r_h = CHOPPED((r));
 696         w_h = r_h * ss_h + hln2pim1_h;
 697         z2 = z * z;
 698         w = (r - r_h) * ss_h + r * ss_l;
 699         t1 = GP[19];

 700         for (i = 18; i > 0; i--)
 701                 t1 = z2 * t1 + GP[i];

 702         w += hln2pim1_l;
 703         w_l = z * (GP[0] + z2 * t1) + w;
 704         k = (int) ((w_h + w_l) * invln2_32 + half);
 705 
 706         /* compute the exponential of w_h+w_l */
 707 
 708         j = k & 0x1f;
 709         *m = k >> 5;
 710         t3 = (long double) k;
 711 
 712         /* perform w - k*ln2_32 (represent as w_h - w_l) */
 713         t1 = w_h - t3 * ln2_32hi;
 714         t2 = t3 * ln2_32lo;
 715         w = t2 - w_l;
 716         w_h = t1 - w;
 717         w_l = w - (t1 - w_h);
 718 
 719         /* compute exp(w_h-w_l) */
 720         z = w_h - w_l;

 721         for (t1 = Et[10], i = 9; i >= 0; i--)
 722                 t1 = z * t1 + Et[i];

 723         t3 = w_h - (w_l - (z * z) * t1);        /* t3 = expm1(z) */
 724         zz.l = S_trail[j] * (one + t3) + S[j] * t3;
 725         zz.h = S[j];
 726         return (zz);
 727 }
 728 
 729 /* INDENT OFF */
 730 /*
 731  * kpsin(x)= sin(pi*x)/pi
 732  *                 3        5        7        9        11                27
 733  *      = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x  + ... + ks[12]*x
 734  */
 735 static const long double ks[] = {
 736         -1.64493406684822643647241516664602518705158902870e+0000L,
 737         +8.11742425283353643637002772405874238094995726160e-0001L,
 738         -1.90751824122084213696472111835337366232282723933e-0001L,
 739         +2.61478478176548005046532613563241288115395517084e-0002L,
 740         -2.34608103545582363750893072647117829448016479971e-0003L,
 741         +1.48428793031071003684606647212534027556262040158e-0004L,
 742         -6.97587366165638046518462722252768122615952898698e-0006L,
 743         +2.53121740413702536928659271747187500934840057929e-0007L,
 744         -7.30471182221385990397683641695766121301933621956e-0009L,
 745         +1.71653847451163495739958249695549313987973589884e-0010L,
 746         -3.34813314714560776122245796929054813458341420565e-0012L,
 747         +5.50724992262622033449487808306969135431411753047e-0014L,
 748         -7.67678132753577998601234393215802221104236979928e-0016L,
 749 };
 750 /* INDENT ON */
 751 
 752 /*
 753  * assume x is not tiny and positive
 754  */
 755 static struct LDouble
 756 kpsin(long double x) {

 757         long double z, t1, t2;
 758         struct LDouble xx;
 759         int i;
 760 
 761         z = x * x;
 762         xx.h = x;

 763         for (t2 = ks[12], i = 11; i > 0; i--)
 764                 t2 = z * t2 + ks[i];

 765         t1 = z * x;
 766         t2 *= z * t1;
 767         xx.l = t1 * ks[0] + t2;
 768         return (xx);
 769 }
 770 
 771 /* INDENT OFF */
 772 /*
 773  * kpcos(x)= cos(pi*x)/pi
 774  *                     2        4        6        8        10        12
 775  *      = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x  +kc[5]*x
 776  *
 777  *                     2        4        6        8        10            22
 778  *      = 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x  +...+kc[9]*x
 779  *
 780  * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l)
 781  *         =  npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x
 782  *         =  npi_2_h*x_f*x_f + npi_2_h*(x*x-x_f*x_f) + npi_2_l*x*x
 783  *         =  npi_2_h*x_f*x_f + npi_2_h*(x+x_f)*(x-x_f) + npi_2_l*x*x
 784  * Here x_f = (long double) (float)x
 785  * Note that pi/2(in hex) =
 786  *  1.921FB54442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29
 787  * npi_2_h = -pi/2 chopped to 25 bits = -1.921FB50000000000000000000000000 =
 788  *  -1.570796310901641845703125000000000 and
 789  * npi_2_l =
 790  *  -0.0000004442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 =
 791  *  -.0000000158932547735281966916397514420985846996875529104874722961539 =
 792  *  -1.5893254773528196691639751442098584699687552910487472296153e-8
 793  * 1/pi(in hex) =
 794  *  .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 795  * will be splitted into:
 796  *  one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000...  and
 797  *  one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 798  */
 799 
 800 static const long double
 801 #if defined(__x86)
 802 one_pi_h = 0.3183098861481994390487670898437500L,       /* 31 bits */
 803 one_pi_l = 3.559123248900043690127872406891929148e-11L,
 804 #else
 805 one_pi_h = 0.31830988618379052468299050815403461456298828125L,
 806 one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
 807 #endif
 808 npi_2_h = -1.570796310901641845703125000000000L,
 809 npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;
 810 
 811 static const long double kc[] = {
 812         +1.29192819501249250731151312779548918765320728489e+0000L,
 813         -4.25027339979557573976029596929319207009444090366e-0001L,
 814         +7.49080661650990096109672954618317623888421628613e-0002L,
 815         -8.21458866111282287985539464173976555436050215120e-0003L,
 816         +6.14202578809529228503205255165761204750211603402e-0004L,
 817         -3.33073432691149607007217330302595267179545908740e-0005L,
 818         +1.36970959047832085796809745461530865597993680204e-0006L,
 819         -4.41780774262583514450246512727201806217271097336e-0008L,
 820         +1.14741409212381858820016567664488123478660705759e-0009L,
 821         -2.44261236114707374558437500654381006300502749632e-0011L,
 822 };
 823 /* INDENT ON */
 824 
 825 /*
 826  * assume x is not tiny and positive
 827  */
 828 static struct LDouble
 829 kpcos(long double x) {

 830         long double z, t1, t2, t3, t4, x4, x8;
 831         int i;
 832         struct LDouble xx;
 833 
 834         z = x * x;
 835         xx.h = one_pi_h;
 836         t1 = (long double) ((float) x);
 837         x4 = z * z;
 838         t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);

 839         for (i = 8, t3 = kc[9]; i >= 0; i--)
 840                 t3 = z * t3 + kc[i];

 841         t3 = one_pi_l + x4 * t3;
 842         t4 = t1 * t1 * npi_2_h;
 843         x8 = t2 + t3;
 844         xx.l = x8 + t4;
 845         return (xx);
 846 }
 847 
 848 /* INDENT OFF */
 849 static const long double
 850         /* 0.13486180573279076968979393577465291700642511139552429398233 */
 851 #if defined(__x86)
 852 t0z1   =  0.1348618057327907696779385054997035808810L,
 853 t0z1_l =  1.1855430274949336125392717150257379614654e-20L,
 854 #else
 855 t0z1   =  0.1348618057327907696897939357746529168654L,
 856 t0z1_l =  1.4102088588676879418739164486159514674310e-37L,
 857 #endif
 858         /* 0.46163214496836234126265954232572132846819620400644635129599 */
 859 #if defined(__x86)
 860 t0z2   =  0.4616321449683623412538115843295472018326L,
 861 t0z2_l =  8.84795799617412663558532305039261747030640e-21L,
 862 #else
 863 t0z2   =  0.46163214496836234126265954232572132343318L,
 864 t0z2_l =  5.03501162329616380465302666480916271611101e-36L,
 865 #endif
 866         /* 0.81977310110050060178786870492160699631174407846245179119586 */
 867 #if defined(__x86)
 868 t0z3   =  0.81977310110050060178773362329351925836817L,
 869 t0z3_l =  1.350816280877379435658077052534574556256230e-22L
 870 #else
 871 t0z3   =  0.8197731011005006017878687049216069516957449L,
 872 t0z3_l =  4.461599916947014419045492615933551648857380e-35L
 873 #endif
 874 ;
 875 /* INDENT ON */
 876 
 877 /*
 878  * gamma(x+i) for 0 <= x < 1
 879  */
 880 static struct LDouble
 881 gam_n(int i, long double x) {
 882         struct LDouble rr = {0.0L, 0.0L}, yy;

 883         long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
 884 
 885         /* compute yy = gamma(x+1) */
 886         if (x > 0.2845L) {
 887                 if (x > 0.6374L) {
 888                         r1 = x - t0z3;
 889                         r2 = CHOPPED((r1 - t0z3_l));
 890                         t2 = r1 - r2;
 891                         yy = GT3(r2, t2 - t0z3_l);
 892                 } else {
 893                         r1 = x - t0z2;
 894                         r2 = CHOPPED((r1 - t0z2_l));
 895                         t2 = r1 - r2;
 896                         yy = GT2(r2, t2 - t0z2_l);
 897                 }
 898         } else {
 899                 r1 = x - t0z1;
 900                 r2 = CHOPPED((r1 - t0z1_l));
 901                 t2 = r1 - r2;
 902                 yy = GT1(r2, t2 - t0z1_l);
 903         }

 904         /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
 905         switch (i) {
 906         case 0:         /* yy/x */
 907                 r1 = one / x;
 908                 xh = CHOPPED((x));      /* x is not tiny */
 909                 rr.h = CHOPPED(((yy.h + yy.l) * r1));
 910                 rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) -
 911                         r1 * yy.l);
 912                 break;
 913         case 1:         /* yy */
 914                 rr.h = yy.h;
 915                 rr.l = yy.l;
 916                 break;
 917         case 2:         /* (x+1)*yy */
 918                 z = x + one;    /* may not be exact */
 919                 zh = CHOPPED((z));
 920                 rr.h = zh * yy.h;
 921                 rr.l = z * yy.l + (x - (zh - one)) * yy.h;
 922                 break;
 923         case 3:         /* (x+2)*(x+1)*yy */
 924                 z1 = x + one;
 925                 z2 = x + 2.0L;
 926                 z = z1 * z2;
 927                 xh = CHOPPED((z));
 928                 zh = CHOPPED((z1));
 929                 xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
 930 
 931                 rr.h = xh * yy.h;


 959                 z *= z2;
 960                 xh = CHOPPED((z));
 961                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 962                 rr.h = xh * yy.h;
 963                 rr.l = z * yy.l + xl * yy.h;
 964                 break;
 965         case 6:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
 966                 z1 = x + 2.0L;
 967                 z2 = x + 3.0L;
 968                 z = z1 * z2;
 969                 zh = CHOPPED((z1));
 970                 yh = CHOPPED((z));
 971                 z1 = x - (zh - 2.0L);
 972                 yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
 973                 z2 = z - 2.0L;
 974                 x5 = x + 5.0L;
 975                 z *= z2;
 976                 xh = CHOPPED(z);
 977                 zh += 3.0;
 978                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 979                                                 /* xh+xl=(x+1)*...*(x+4) */
 980                 /* wh+wl=(x+5)*yy */



 981                 wh = CHOPPED((x5 * (yy.h + yy.l)));
 982                 wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
 983                 rr.h = wh * xh;
 984                 rr.l = z * wl + xl * wh;
 985                 break;
 986         case 7:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
 987                 z1 = x + 3.0L;
 988                 z2 = x + 4.0L;
 989                 z = z2 * z1;
 990                 zh = CHOPPED((z1));
 991                 yh = CHOPPED((z));      /* yh+yl = (x+3)(x+4) */
 992                 yl = (x - (zh - 3.0L)) * (z2 + zh) - (yh - (zh * (zh + one)));
 993                 z1 = x + 6.0L;
 994                 z2 = z - 2.0L;  /* z2 = (x+2)*(x+5) */
 995                 z *= z2;
 996                 xh = CHOPPED((z));
 997                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 998                                                 /* xh+xl=(x+2)*...*(x+5) */
 999                 /* wh+wl=(x+1)(x+6)*yy */



1000                 z2 -= 4.0L;     /* z2 = (x+1)(x+6) */
1001                 wh = CHOPPED((z2 * (yy.h + yy.l)));
1002                 wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0L) * yy.h);
1003                 rr.h = wh * xh;
1004                 rr.l = z * wl + xl * wh;
1005         }

1006         return (rr);
1007 }
1008 
1009 long double
1010 tgammal(long double x) {

1011         struct LDouble ss, ww;
1012         long double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
1013         int i, j, m, ix, hx, xk;
1014         unsigned lx;
1015 
1016         hx = H0_WORD(x);
1017         lx = H3_WORD(x);
1018         ix = hx & 0x7fffffff;
1019         y = x;
1020         if (ix < 0x3f8e0000) {       /* x < 2**-113 */

1021                 return (one / x);
1022         }
1023         if (ix >= 0x7fff0000)
1024                 return (x * ((hx < 0)? zero : x));   /* Inf or NaN */

1025         if (x > overflow)    /* overflow threshold */
1026                 return (x * 1.0e4932L);

1027         if (hx >= 0x40020000) {      /* x >= 8 */
1028                 ww = large_gam(x, &m);
1029                 w = ww.h + ww.l;
1030                 return (scalbnl(w, m));
1031         }
1032 
1033         if (hx > 0) {                /* 0 < x < 8 */
1034                 i = (int) x;
1035                 ww = gam_n(i, x - (long double) i);
1036                 return (ww.h + ww.l);
1037         }
1038         /* INDENT OFF */
1039         /* negative x */



1040         /*
1041          * compute xk =
1042          *      -2 ... x is an even int (-inf is considered an even #)
1043          *      -1 ... x is an odd int
1044          *      +0 ... x is not an int but chopped to an even int
1045          *      +1 ... x is not an int but chopped to an odd int
1046          */
1047         /* INDENT ON */
1048         xk = 0;
1049 #if defined(__x86)
1050         if (ix >= 0x403e0000) {      /* x >= 2**63 } */
1051                 if (ix >= 0x403f0000)
1052                         xk = -2;
1053                 else
1054                         xk = -2 + (lx & 1);
1055 #else
1056         if (ix >= 0x406f0000) {      /* x >= 2**112 */
1057                 if (ix >= 0x40700000)
1058                         xk = -2;
1059                 else
1060                         xk = -2 + (lx & 1);
1061 #endif
1062         } else if (ix >= 0x3fff0000) {
1063                 w = -x;
1064                 t1 = floorl(w);
1065                 t2 = t1 * half;
1066                 t3 = floorl(t2);

1067                 if (t1 == w) {
1068                         if (t2 == t3)
1069                                 xk = -2;
1070                         else
1071                                 xk = -1;
1072                 } else {
1073                         if (t2 == t3)
1074                                 xk = 0;
1075                         else
1076                                 xk = 1;
1077                 }
1078         }
1079 
1080         if (xk < 0) {
1081                 /* return NaN. Ideally gamma(-n)= (-1)**(n+1) * inf */
1082                 return (x - x) / (x - x);
1083         }
1084 
1085         /*
1086          * negative underflow thresold -(1774+9ulp)
1087          */
1088         if (x < -1774.0000000000000000000000000000017749370L) {
1089                 z = tiny / x;

1090                 if (xk == 1)
1091                         z = -z;

1092                 return (z * tiny);
1093         }
1094 
1095         /* INDENT OFF */
1096         /*
1097          * now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
1098          */

1099         /*
1100          * First compute ss = -sin(pi*y)/pi so that
1101          * gamma(x) = 1/(ss*gamma(1+y))
1102          */
1103         /* INDENT ON */
1104         y = -x;
1105         j = (int) y;
1106         z = y - (long double) j;
1107         if (z > 0.3183098861837906715377675L)

1108                 if (z > 0.6816901138162093284622325L)
1109                         ss = kpsin(one - z);
1110                 else
1111                         ss = kpcos(0.5L - z);
1112         else
1113                 ss = kpsin(z);


1114         if (xk == 0) {
1115                 ss.h = -ss.h;
1116                 ss.l = -ss.l;
1117         }
1118 
1119         /* Then compute ww = gamma(1+y), note that result scale to 2**m */
1120         m = 0;

1121         if (j < 7) {
1122                 ww = gam_n(j + 1, z);
1123         } else {
1124                 w = y + one;

1125                 if ((lx & 1) == 0) {        /* y+1 exact (note that y<184) */
1126                         ww = large_gam(w, &m);
1127                 } else {
1128                         t = w - one;

1129                         if (t == y) {   /* y+one exact */
1130                                 ww = large_gam(w, &m);
1131                         } else {        /* use y*gamma(y) */
1132                                 if (j == 7)
1133                                         ww = gam_n(j, z);
1134                                 else
1135                                         ww = large_gam(y, &m);

1136                                 t4 = ww.h + ww.l;
1137                                 t1 = CHOPPED((y));
1138                                 t2 = CHOPPED((t4));
1139                                                 /* t4 will not be too large */
1140                                 ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
1141                                 ww.h = t1 * t2;
1142                         }
1143                 }
1144         }
1145 
1146         /* compute 1/(ss*ww) */
1147         t3 = ss.h + ss.l;
1148         t4 = ww.h + ww.l;
1149         t1 = CHOPPED((t3));
1150         t2 = CHOPPED((t4));
1151         z1 = ss.l - (t1 - ss.h);        /* (t1,z1) = ss */
1152         z2 = ww.l - (t2 - ww.h);        /* (t2,z2) = ww */
1153         t3 = t3 * t4;                   /* t3 = ss*ww */
1154         z3 = one / t3;                  /* z3 = 1/(ss*ww) */
1155         t5 = t1 * t2;


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __tgammal = tgammal
  32 
  33 #include "libm.h"
  34 #include <sys/isa_defs.h>
  35 
  36 #if defined(_BIG_ENDIAN)
  37 #define H0_WORD(x)              ((unsigned *)&x)[0]
  38 #define H3_WORD(x)              ((unsigned *)&x)[3]
  39 #define CHOPPED(x)              (long double)((double)(x))
  40 #else
  41 #define H0_WORD(x)              ((((int *)&x)[2] << 16) | (0x0000ffff & \
  42         (((unsigned *)&x)[1] >> 15)))
  43 #define H3_WORD(x)              ((unsigned *)&x)[0]
  44 #define CHOPPED(x)              (long double)((float)(x))
  45 #endif
  46 
  47 struct LDouble {
  48         long double h, l;
  49 };
  50 
  51 /*
  52  * Primary interval GTi()
  53  */
  54 static const long double P1[] = {
  55         +0.709086836199777919037185741507610124611513720557L,
  56         +4.45754781206489035827915969367354835667391606951e-0001L,
  57         +3.21049298735832382311662273882632210062918153852e-0002L,
  58         -5.71296796342106617651765245858289197369688864350e-0003L,
  59         +6.04666892891998977081619174969855831606965352773e-0003L,
  60         +8.99106186996888711939627812174765258822658645168e-0004L,
  61         -6.96496846144407741431207008527018441810175568949e-0005L,
  62         +1.52597046118984020814225409300131445070213882429e-0005L,
  63         +5.68521076168495673844711465407432189190681541547e-0007L,
  64         +3.30749673519634895220582062520286565610418952979e-0008L,
  65 };
  66 
  67 static const long double Q1[] = {
  68         +1.0 + 0000L,
  69         +1.35806511721671070408570853537257079579490650668e+0000L,
  70         +2.97567810153429553405327140096063086994072952961e-0001L,
  71         -1.52956835982588571502954372821681851681118097870e-0001L,
  72         -2.88248519561420109768781615289082053597954521218e-0002L,
  73         +1.03475311719937405219789948456313936302378395955e-0002L,
  74         +4.12310203243891222368965360124391297374822742313e-0004L,
  75         -3.12653708152290867248931925120380729518332507388e-0004L,
  76         +2.36672170850409745237358105667757760527014332458e-0005L,
  77 };
  78 
  79 static const long double P2[] = {
  80         +0.428486815855585429730209907810650135255270600668084114L,
  81         +2.62768479103809762805691743305424077975230551176e-0001L,
  82         +3.81187532685392297608310837995193946591425896150e-0002L,
  83         +3.00063075891811043820666846129131255948527925381e-0003L,
  84         +2.47315407812279164228398470797498649142513408654e-0003L,
  85         +3.62838199917848372586173483147214880464782938664e-0004L,
  86         +3.43991105975492623982725644046473030098172692423e-0006L,
  87         +4.56902151569603272237014240794257659159045432895e-0006L,
  88         +2.13734755837595695602045100675540011352948958453e-0007L,
  89         +9.74123440547918230781670266967882492234877125358e-0009L,
  90 };
  91 
  92 static const long double Q2[] = {
  93         +1.0L,
  94         +9.18284118632506842664645516830761489700556179701e-0001L,
  95         -6.41430858837830766045202076965923776189154874947e-0003L,
  96         -1.24400885809771073213345747437964149775410921376e-0001L,
  97         +4.69803798146251757538856567522481979624746875964e-0003L,
  98         +7.18309447069495315914284705109868696262662082731e-0003L,
  99         -8.75812626987894695112722600697653425786166399105e-0004L,
 100         -1.23539972377769277995959339188431498626674835169e-0004L,
 101         +3.10019017590151598732360097849672925448587547746e-0005L,
 102         -1.77260223349332617658921874288026777465782364070e-0006L,
 103 };
 104 
 105 static const long double P3[] = {
 106         +0.3824094797345675048502747661075355640070439388902L,
 107         +3.42198093076618495415854906335908427159833377774e-0001L,
 108         +9.63828189500585568303961406863153237440702754858e-0002L,
 109         +8.76069421042696384852462044188520252156846768667e-0003L,
 110         +1.86477890389161491224872014149309015261897537488e-0003L,
 111         +8.16871354540309895879974742853701311541286944191e-0004L,
 112         +6.83783483674600322518695090864659381650125625216e-0005L,
 113         -1.10168269719261574708565935172719209272190828456e-0006L,
 114         +9.66243228508380420159234853278906717065629721016e-0007L,
 115         +2.31858885579177250541163820671121664974334728142e-0008L,
 116 };
 117 
 118 static const long double Q3[] = {
 119         +1.0L, +8.25479821168813634632437430090376252512793067339e-0001L,

 120         -1.62251363073937769739639623669295110346015576320e-0002L,
 121         -1.10621286905916732758745130629426559691187579852e-0001L,
 122         +3.48309693970985612644446415789230015515365291459e-0003L,
 123         +6.73553737487488333032431261131289672347043401328e-0003L,
 124         -7.63222008393372630162743587811004613050245128051e-0004L,
 125         -1.35792670669190631476784768961953711773073251336e-0004L,
 126         +3.19610150954223587006220730065608156460205690618e-0005L,
 127         -1.82096553862822346610109522015129585693354348322e-0006L,
 128 };
 129 
 130 static const long double
 131 #if defined(__x86)
 132         GZ1_h = 0.938204627909682449364570100414084663498215377L,
 133         GZ1_l = 4.518346116624229420055327632718530617227944106e-20L,
 134         GZ2_h = 0.885603194410888700264725126309883762587560340L,
 135         GZ2_l = 1.409077427270497062039119290776508217077297169e-20L,
 136         GZ3_h = 0.936781411463652321613537060640553022494714241L,
 137         GZ3_l = 5.309836440284827247897772963887219035221996813e-21L,
 138 #else
 139         GZ1_h = 0.938204627909682449409753561580326910854647031L,
 140         GZ1_l = 4.684412162199460089642452580902345976446297037e-35L,
 141         GZ2_h = 0.885603194410888700278815900582588658192658794L,
 142         GZ2_l = 7.501529273890253789219935569758713534641074860e-35L,
 143         GZ3_h = 0.936781411463652321618846897080837818855399840L,
 144         GZ3_l = 3.088721217404784363585591914529361687403776917e-35L,
 145 #endif
 146         TZ1 = -0.3517214357852935791015625L,
 147         TZ3 = 0.280530631542205810546875L;
 148 
 149 

 150 /*
 151  * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
 152  * ...assume yh got 53 or 24(i386) significant bits
 153  */

 154 static struct LDouble
 155 GT1(long double yh, long double yl)
 156 {
 157         long double t3, t4, y;
 158         int i;
 159         struct LDouble r;
 160 
 161         y = yh + yl;
 162 
 163         for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) {
 164                 t4 = t4 * y + Q1[i];
 165                 t3 = t3 * y + P1[i];
 166         }
 167 
 168         t3 = (y * y) * t3 / t4;
 169         t3 += (TZ1 * yl + GZ1_l);
 170         t4 = TZ1 * yh;
 171         r.h = CHOPPED((t4 + GZ1_h + t3));
 172         t3 += (t4 - (r.h - GZ1_h));
 173         r.l = t3;
 174         return (r);
 175 }
 176 
 177 
 178 /*
 179  * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
 180  * ...assume yh got 53 significant bits
 181  */

 182 static struct LDouble
 183 GT2(long double yh, long double yl)
 184 {
 185         long double t3, t4, y;
 186         int i;
 187         struct LDouble r;
 188 
 189         y = yh + yl;
 190 
 191         for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) {
 192                 t4 = t4 * y + Q2[i];
 193                 t3 = t3 * y + P2[i];
 194         }
 195 
 196         t3 = GZ2_l + (y * y) * t3 / t4;
 197         r.h = CHOPPED((GZ2_h + t3));
 198         r.l = t3 - (r.h - GZ2_h);
 199         return (r);
 200 }
 201 
 202 
 203 /*
 204  * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
 205  * ...assume yh got 53 significant bits
 206  */

 207 static struct LDouble
 208 GT3(long double yh, long double yl)
 209 {
 210         long double t3, t4, y;
 211         int i;
 212         struct LDouble r;
 213 
 214         y = yh + yl;
 215 
 216         for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) {
 217                 t4 = t4 * y + Q3[i];
 218                 t3 = t3 * y + P3[i];
 219         }
 220 
 221         t3 = (y * y) * t3 / t4;
 222         t3 += (TZ3 * yl + GZ3_l);
 223         t4 = TZ3 * yh;
 224         r.h = CHOPPED((t4 + GZ3_h + t3));
 225         t3 += (t4 - (r.h - GZ3_h));
 226         r.l = t3;
 227         return (r);
 228 }
 229 
 230 /*
 231  * Hex value of GP[0] shoule be 3FB55555 55555555
 232  */
 233 static const long double GP[] = {
 234         +0.083333333333333333333333333333333172839171301L,
 235         -2.77777777777777777777777777492501211999399424104e-0003L,
 236         +7.93650793650793650793635650541638236350020883243e-0004L,
 237         -5.95238095238095238057299772679324503339241961704e-0004L,
 238         +8.41750841750841696138422987977683524926142600321e-0004L,
 239         -1.91752691752686682825032547823699662178842123308e-0003L,
 240         +6.41025641022403480921891559356473451161279359322e-0003L,
 241         -2.95506535798414019189819587455577003732808185071e-0002L,
 242         +1.79644367229970031486079180060923073476568732136e-0001L,
 243         -1.39243086487274662174562872567057200255649290646e+0000L,
 244         +1.34025874044417962188677816477842265259608269775e+0001L,
 245         -1.56803713480127469414495545399982508700748274318e+0002L,
 246         +2.18739841656201561694927630335099313968924493891e+0003L,
 247         -3.55249848644100338419187038090925410976237921269e+0004L,
 248         +6.43464880437835286216768959439484376449179576452e+0005L,
 249         -1.20459154385577014992600342782821389605893904624e+0007L,
 250         +2.09263249637351298563934942349749718491071093210e+0008L,
 251         -2.96247483183169219343745316433899599834685703457e+0009L,
 252         +2.88984933605896033154727626086506756972327292981e+0010L,
 253         -1.40960434146030007732838382416230610302678063984e+0011L, /* 19 */
 254 };
 255 
 256 static const long double T3[] = {
 257         +0.666666666666666666666666666666666634567834260213L, /* T3[0] */
 258         +0.400000000000000000000000000040853636176634934140L, /* T3[1] */
 259         +0.285714285714285714285696975252753987869020263448L, /* T3[2] */
 260         +0.222222222222222225593221101192317258554772129875L, /* T3[3] */
 261         +0.181818181817850192105847183461778186703779262916L, /* T3[4] */
 262         +0.153846169861348633757101285952333369222567014596L, /* T3[5] */
 263         +0.133033462889260193922261296772841229985047571265L, /* T3[6] */
 264 };
 265 /* BEGIN CSTYLED */
 266 static const long double c[] = {
 267         0.0L,
 268         1.0L,
 269         2.0L,
 270         0.5L,
 271         1.0e-4930L,                  /* tiny */
 272         4.18937683105468750000e-01L, /* hln2pim1_h */
 273         8.50099203991780329736405617639861397473637783412817152e-07L, /* hln2pim1_l */
 274         0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
 275         2.16608493865351192653179168701171875e-02L, /* ln2_32hi */
 276         5.96317165397058692545083025235937919875797669127130e-12L, /* ln2_32lo */
 277         46.16624130844682903551758979206054839765267053289554989233L, /* invln2_32 */
 278 #if defined(__x86)
 279         1.7555483429044629170023839037639845628291e+03L, /* overflow */
 280 #else
 281         1.7555483429044629170038892160702032034177e+03L, /* overflow */
 282 #endif
 283 };
 284 /* END CSTYLED */
 285 
 286 #define zero                    c[0]
 287 #define one                     c[1]
 288 #define two                     c[2]
 289 #define half                    c[3]
 290 #define tiny                    c[4]
 291 #define hln2pim1_h              c[5]
 292 #define hln2pim1_l              c[6]
 293 #define hln2pim1                c[7]
 294 #define ln2_32hi                c[8]
 295 #define ln2_32lo                c[9]
 296 #define invln2_32               c[10]
 297 #define overflow                c[11]
 298 
 299 /*
 300  * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64
 301  */
 302 static const long double Et[] = {
 303         +5.0000000000000000000e-1L,
 304         +1.66666666666666666666666666666828835166292152466e-0001L,


 310         +2.75573192239028921114572986441972140933432317798e-0006L,
 311         +2.75573192239448470555548102895526369739856219317e-0007L,
 312         +2.50521677867683935940853997995937600214167232477e-0008L,
 313         +2.08767928899010367374984448513685566514152147362e-0009L,
 314 };
 315 
 316 /*
 317  * long double precision coefficients for computing log(x)-1 in tgamma.
 318  *  See "algorithm" for details
 319  *
 320  *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 321  *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 322  *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 323  *       T2(j) = T2[2j,2j+1] = log(z[j]),
 324  *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15
 325  *  Note
 326  *  (1) the leading entries are truncated to 24 binary point.
 327  *  (2) Remez error for T3(s) is bounded by 2**(-136.54)
 328  */
 329 static const long double T1[] = {
 330         -1.000000000000000000000000000000000000000000e+00L,
 331         +0.000000000000000000000000000000000000000000e+00L,
 332         -3.068528175354003906250000000000000000000000e-01L,
 333         -1.904654299957767878541823431924500011926579e-09L,
 334         +3.862943053245544433593750000000000000000000e-01L,
 335         +5.579533617547508924291635313615100141107647e-08L,
 336         +1.079441487789154052734375000000000000000000e+00L,
 337         +5.389068187551732136437452970422650211661470e-08L,
 338         +1.772588670253753662109375000000000000000000e+00L,
 339         +5.198602757555955348583270627230200282215294e-08L,
 340         +2.465735852718353271484375000000000000000000e+00L,
 341         +5.008137327560178560729088284037750352769117e-08L,
 342         +3.158883035182952880859375000000000000000000e+00L,
 343         +4.817671897564401772874905940845299849351090e-08L,
 344         +3.852030217647552490234375000000000000000000e+00L,
 345         +4.627206467568624985020723597652849919904913e-08L,
 346         +4.545177400112152099609375000000000000000000e+00L,
 347         +4.436741037572848197166541254460399990458737e-08L,
 348         +5.238324582576751708984375000000000000000000e+00L,
 349         +4.246275607577071409312358911267950061012560e-08L,
 350         +5.931471765041351318359375000000000000000000e+00L,
 351         +4.055810177581294621458176568075500131566384e-08L,
 352 };
 353 


 540         +1.35425554693689272829801474014070273e+00L,
 541         +1.38390988196383195487265952726519287e+00L,
 542         +1.41421356237309504880168872420969798e+00L,
 543         +1.44518080697704662003700624147167095e+00L,
 544         +1.47682614593949931138690748037404985e+00L,
 545         +1.50916442759342273976601955103319352e+00L,
 546         +1.54221082540794082361229186209073479e+00L,
 547         +1.57598084510788648645527016018190504e+00L,
 548         +1.61049033194925430817952066735740067e+00L,
 549         +1.64575547815396484451875672472582254e+00L,
 550         +1.68179283050742908606225095246642969e+00L,
 551         +1.71861929812247791562934437645631244e+00L,
 552         +1.75625216037329948311216061937531314e+00L,
 553         +1.79470907500310718642770324212778174e+00L,
 554         +1.83400808640934246348708318958828892e+00L,
 555         +1.87416763411029990132999894995444645e+00L,
 556         +1.91520656139714729387261127029583086e+00L,
 557         +1.95714412417540026901832225162687149e+00L,
 558 #endif
 559 };
 560 
 561 static const long double S_trail[] = {
 562 #if defined(__x86)
 563         +0.0000000000000000000000000e+00L,
 564         +2.6327965667180882569382524e-20L,
 565         +8.3765863521895191129661899e-20L,
 566         +3.9798705777454504249209575e-20L,
 567         +1.0668046596651558640993042e-19L,
 568         +1.9376009847285360448117114e-20L,
 569         +6.7081819456112953751277576e-21L,
 570         +1.9711680502629186462729727e-20L,
 571         +2.9932584438449523689104569e-20L,
 572         +6.8887754153039109411061914e-20L,
 573         +6.8002718741225378942847820e-20L,
 574         +6.5846917376975403439742349e-20L,
 575         +1.2171958727511372194876001e-20L,
 576         +3.5625253228704087115438260e-20L,
 577         +3.1129551559077560956309179e-20L,
 578         +5.7519192396164779846216492e-20L,
 579         +3.7900651177865141593101239e-20L,
 580         +1.1659262405698741798080115e-20L,
 581         +7.1364385105284695967172478e-20L,
 582         +5.2631003710812203588788949e-20L,
 583         +2.6328853788732632868460580e-20L,
 584         +5.4583950085438242788190141e-20L,
 585         +9.5803254376938269960718656e-20L,
 586         +7.6837733983874245823512279e-21L,
 587         +2.4415965910835093824202087e-20L,
 588         +2.6052966871016580981769728e-20L,
 589         +2.6876456344632553875309579e-21L,
 590         +1.2861930155613700201703279e-20L,
 591         +8.8166633394037485606572294e-20L,
 592         +2.9788615389580190940837037e-20L,
 593         +5.2352341619805098677422139e-20L,
 594         +5.2578463064010463732242363e-20L,
 595 #else
 596         +0.00000000000000000000000000000000000e+00L,
 597         +1.80506787420330954745573333054573786e-35L,
 598         -9.37452029228042742195756741973083214e-35L,
 599         -1.59696844729275877071290963023149997e-35L,
 600         +9.11249341012502297851168610167248666e-35L,
 601         -6.50422820697854828723037477525938871e-35L,
 602         -8.14846884452585113732569176748815532e-35L,
 603         -5.06621457672180031337233074514290335e-35L,
 604         -1.35983097468881697374987563824591912e-35L,
 605         +9.49742763556319647030771056643324660e-35L,
 606         -3.28317052317699860161506596533391526e-36L,
 607         -5.01723570938719041029018653045842895e-35L,
 608         -2.39147479768910917162283430160264014e-35L,
 609         -8.35057135763390881529889073794408385e-36L,
 610         +7.03675688907326504242173719067187644e-35L,
 611         -5.18248485306464645753689301856695619e-35L,
 612         +9.42224254862183206569211673639406488e-35L,
 613         -3.96750082539886230916730613021641828e-35L,
 614         +7.14352899156330061452327361509276724e-35L,
 615         +1.15987125286798512424651783410044433e-35L,
 616         +4.69693347835811549530973921320187447e-35L,
 617         -3.38651317599500471079924198499981917e-35L,
 618         -8.58731877429824706886865593510387445e-35L,
 619         -9.60595154874935050318549936224606909e-35L,
 620         +9.60973393212801278450755869714178581e-35L,
 621         +6.37839792144002843924476144978084855e-35L,
 622         +7.79243078569586424945646112516927770e-35L,
 623         +7.36133776758845652413193083663393220e-35L,
 624         -6.47299514791334723003521457561217053e-35L,
 625         +8.58747441795369869427879806229522962e-35L,
 626         +2.37181542282517483569165122830269098e-35L,
 627         -3.02689168209611877300459737342190031e-37L,
 628 #endif
 629 };

 630 
 631 
 632 /*
 633  * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
 634  *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
 635  *                = L1 + L2 + L3,
 636  */

 637 static struct LDouble
 638 large_gam(long double x, int *m)
 639 {
 640         long double z, t1, t2, t3, z2, t5, w, y, u, r, v;
 641         long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L;
 642         int n2, j2, k, ix, j, i;
 643         struct LDouble zz;
 644         long double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
 645 
 646 /* BEGIN CSTYLED */
 647 /*
 648  * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
 649  *
 650  *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 651  *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 652  *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 653  *       T2(j) = T2[2j,2j+1] = log(z[j]),
 654  *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15
 655  *  Note
 656  *  (1) the leading entries are truncated to 24 binary point.
 657  *  (2) Remez error for T3(s) is bounded by 2**(-72.4)
 658  *                                   2**(-24)
 659  *                           _________V___________________
 660  *               T1(n):     |_________|___________________|
 661  *                             _______ ______________________
 662  *               T2(j):       |_______|______________________|
 663  *                                ____ _______________________
 664  *               2s:             |____|_______________________|
 665  *                                    __________________________
 666  *          +    T3(s)-2s:           |__________________________|
 667  *                       -------------------------------------------
 668  *                          [leading] + [Trailing]
 669  */
 670         /* END CSTYLED */
 671         ix = H0_WORD(x);
 672         n2 = (ix >> 16) - 0x3fff; /* exponent of x, range:3-10 */
 673         y = scalbnl(x, -n2);      /* y = scale x to [1,2] */
 674         n2 += n2;                 /* 2n */
 675         j = (ix >> 10) & 0x3f;          /* j */
 676         z = 1.0078125L + (long double)j * 0.015625L;    /* z[j]=1+j/64+1/128 */
 677         j2 = j + j;
 678         t1 = y + z;
 679         t2 = y - z;
 680         r = one / t1;
 681         u = r * t2;                     /* u = (y-z)/(y+z) */
 682         t1 = CHOPPED(t1);
 683         t4 = T2[j2 + 1] + T1[n2 + 1];
 684         z2 = u * u;
 685         k = H0_WORD(u) & 0x7fffffff;
 686         t3 = T2[j2] + T1[n2];
 687 
 688         for (t5 = T3[6], i = 5; i >= 0; i--)
 689                 t5 = z2 * t5 + T3[i];
 690 
 691         if ((k >> 16) < 0x3fec) {      /* |u|<2**-19 */
 692                 t2 = t4 + u * (two + z2 * t5);
 693         } else {
 694                 t5 = t4 + (u * z2) * t5;
 695                 u2 = u + u;
 696                 v = (long double)((int)(u2 * t24)) * p24;
 697                 t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
 698                 t3 += v;
 699         }
 700 
 701         ss_h = CHOPPED((t2 + t3));
 702         ss_l = t2 - (ss_h - t3);
 703 
 704 /*
 705  * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
 706  * where ss = log(x) - 1 in already in extra precision
 707  */

 708         z = one / x;
 709         r = x - half;
 710         r_h = CHOPPED((r));
 711         w_h = r_h * ss_h + hln2pim1_h;
 712         z2 = z * z;
 713         w = (r - r_h) * ss_h + r * ss_l;
 714         t1 = GP[19];
 715 
 716         for (i = 18; i > 0; i--)
 717                 t1 = z2 * t1 + GP[i];
 718 
 719         w += hln2pim1_l;
 720         w_l = z * (GP[0] + z2 * t1) + w;
 721         k = (int)((w_h + w_l) * invln2_32 + half);
 722 
 723         /* compute the exponential of w_h+w_l */
 724 
 725         j = k & 0x1f;
 726         *m = k >> 5;
 727         t3 = (long double)k;
 728 
 729         /* perform w - k*ln2_32 (represent as w_h - w_l) */
 730         t1 = w_h - t3 * ln2_32hi;
 731         t2 = t3 * ln2_32lo;
 732         w = t2 - w_l;
 733         w_h = t1 - w;
 734         w_l = w - (t1 - w_h);
 735 
 736         /* compute exp(w_h-w_l) */
 737         z = w_h - w_l;
 738 
 739         for (t1 = Et[10], i = 9; i >= 0; i--)
 740                 t1 = z * t1 + Et[i];
 741 
 742         t3 = w_h - (w_l - (z * z) * t1);        /* t3 = expm1(z) */
 743         zz.l = S_trail[j] * (one + t3) + S[j] * t3;
 744         zz.h = S[j];
 745         return (zz);
 746 }
 747 
 748 
 749 /*
 750  * kpsin(x)= sin(pi*x)/pi
 751  *                 3        5        7        9        11                27
 752  *      = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x  + ... + ks[12]*x
 753  */
 754 static const long double ks[] = {
 755         -1.64493406684822643647241516664602518705158902870e+0000L,
 756         +8.11742425283353643637002772405874238094995726160e-0001L,
 757         -1.90751824122084213696472111835337366232282723933e-0001L,
 758         +2.61478478176548005046532613563241288115395517084e-0002L,
 759         -2.34608103545582363750893072647117829448016479971e-0003L,
 760         +1.48428793031071003684606647212534027556262040158e-0004L,
 761         -6.97587366165638046518462722252768122615952898698e-0006L,
 762         +2.53121740413702536928659271747187500934840057929e-0007L,
 763         -7.30471182221385990397683641695766121301933621956e-0009L,
 764         +1.71653847451163495739958249695549313987973589884e-0010L,
 765         -3.34813314714560776122245796929054813458341420565e-0012L,
 766         +5.50724992262622033449487808306969135431411753047e-0014L,
 767         -7.67678132753577998601234393215802221104236979928e-0016L,
 768 };

 769 
 770 /*
 771  * assume x is not tiny and positive
 772  */
 773 static struct LDouble
 774 kpsin(long double x)
 775 {
 776         long double z, t1, t2;
 777         struct LDouble xx;
 778         int i;
 779 
 780         z = x * x;
 781         xx.h = x;
 782 
 783         for (t2 = ks[12], i = 11; i > 0; i--)
 784                 t2 = z * t2 + ks[i];
 785 
 786         t1 = z * x;
 787         t2 *= z * t1;
 788         xx.l = t1 * ks[0] + t2;
 789         return (xx);
 790 }
 791 
 792 
 793 /*
 794  * kpcos(x)= cos(pi*x)/pi
 795  *                     2        4        6        8        10        12
 796  *      = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x  +kc[5]*x
 797  *
 798  *                     2        4        6        8        10            22
 799  *      = 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x  +...+kc[9]*x
 800  *
 801  * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l)
 802  *         =  npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x
 803  *         =  npi_2_h*x_f*x_f + npi_2_h*(x*x-x_f*x_f) + npi_2_l*x*x
 804  *         =  npi_2_h*x_f*x_f + npi_2_h*(x+x_f)*(x-x_f) + npi_2_l*x*x
 805  * Here x_f = (long double) (float)x
 806  * Note that pi/2(in hex) =
 807  *  1.921FB54442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29
 808  * npi_2_h = -pi/2 chopped to 25 bits = -1.921FB50000000000000000000000000 =
 809  *  -1.570796310901641845703125000000000 and
 810  * npi_2_l =
 811  *  -0.0000004442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 =
 812  *  -.0000000158932547735281966916397514420985846996875529104874722961539 =
 813  *  -1.5893254773528196691639751442098584699687552910487472296153e-8
 814  * 1/pi(in hex) =
 815  *  .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 816  * will be splitted into:
 817  *  one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000...  and
 818  *  one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 819  */
 820 
 821 static const long double
 822 #if defined(__x86)
 823         one_pi_h = 0.3183098861481994390487670898437500L,       /* 31 bits */
 824         one_pi_l = 3.559123248900043690127872406891929148e-11L,
 825 #else
 826         one_pi_h = 0.31830988618379052468299050815403461456298828125L,
 827         one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
 828 #endif
 829         npi_2_h = -1.570796310901641845703125000000000L,
 830         npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;

 831 static const long double kc[] = {
 832         +1.29192819501249250731151312779548918765320728489e+0000L,
 833         -4.25027339979557573976029596929319207009444090366e-0001L,
 834         +7.49080661650990096109672954618317623888421628613e-0002L,
 835         -8.21458866111282287985539464173976555436050215120e-0003L,
 836         +6.14202578809529228503205255165761204750211603402e-0004L,
 837         -3.33073432691149607007217330302595267179545908740e-0005L,
 838         +1.36970959047832085796809745461530865597993680204e-0006L,
 839         -4.41780774262583514450246512727201806217271097336e-0008L,
 840         +1.14741409212381858820016567664488123478660705759e-0009L,
 841         -2.44261236114707374558437500654381006300502749632e-0011L,
 842 };

 843 
 844 /*
 845  * assume x is not tiny and positive
 846  */
 847 static struct LDouble
 848 kpcos(long double x)
 849 {
 850         long double z, t1, t2, t3, t4, x4, x8;
 851         int i;
 852         struct LDouble xx;
 853 
 854         z = x * x;
 855         xx.h = one_pi_h;
 856         t1 = (long double)((float)x);
 857         x4 = z * z;
 858         t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
 859 
 860         for (i = 8, t3 = kc[9]; i >= 0; i--)
 861                 t3 = z * t3 + kc[i];
 862 
 863         t3 = one_pi_l + x4 * t3;
 864         t4 = t1 * t1 * npi_2_h;
 865         x8 = t2 + t3;
 866         xx.l = x8 + t4;
 867         return (xx);
 868 }
 869 

 870 static const long double
 871 /* 0.13486180573279076968979393577465291700642511139552429398233 */
 872 #if defined(__x86)
 873         t0z1 = 0.1348618057327907696779385054997035808810L,
 874         t0z1_l = 1.1855430274949336125392717150257379614654e-20L,
 875 #else
 876         t0z1 = 0.1348618057327907696897939357746529168654L,
 877         t0z1_l = 1.4102088588676879418739164486159514674310e-37L,
 878 #endif
 879 /* 0.46163214496836234126265954232572132846819620400644635129599 */
 880 #if defined(__x86)
 881         t0z2 = 0.4616321449683623412538115843295472018326L,
 882         t0z2_l = 8.84795799617412663558532305039261747030640e-21L,
 883 #else
 884         t0z2 = 0.46163214496836234126265954232572132343318L,
 885         t0z2_l = 5.03501162329616380465302666480916271611101e-36L,
 886 #endif
 887 /* 0.81977310110050060178786870492160699631174407846245179119586 */
 888 #if defined(__x86)
 889         t0z3 = 0.81977310110050060178773362329351925836817L,
 890         t0z3_l = 1.350816280877379435658077052534574556256230e-22L
 891 #else
 892         t0z3 = 0.8197731011005006017878687049216069516957449L,
 893         t0z3_l = 4.461599916947014419045492615933551648857380e-35L
 894 #endif
 895 ;

 896 
 897 /*
 898  * gamma(x+i) for 0 <= x < 1
 899  */
 900 static struct LDouble
 901 gam_n(int i, long double x)
 902 {
 903         struct LDouble rr = { 0.0L, 0.0L }, yy;
 904         long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
 905 
 906         /* compute yy = gamma(x+1) */
 907         if (x > 0.2845L) {
 908                 if (x > 0.6374L) {
 909                         r1 = x - t0z3;
 910                         r2 = CHOPPED((r1 - t0z3_l));
 911                         t2 = r1 - r2;
 912                         yy = GT3(r2, t2 - t0z3_l);
 913                 } else {
 914                         r1 = x - t0z2;
 915                         r2 = CHOPPED((r1 - t0z2_l));
 916                         t2 = r1 - r2;
 917                         yy = GT2(r2, t2 - t0z2_l);
 918                 }
 919         } else {
 920                 r1 = x - t0z1;
 921                 r2 = CHOPPED((r1 - t0z1_l));
 922                 t2 = r1 - r2;
 923                 yy = GT1(r2, t2 - t0z1_l);
 924         }
 925 
 926         /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
 927         switch (i) {
 928         case 0:                         /* yy/x */
 929                 r1 = one / x;
 930                 xh = CHOPPED((x));      /* x is not tiny */
 931                 rr.h = CHOPPED(((yy.h + yy.l) * r1));
 932                 rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) - r1 *
 933                     yy.l);
 934                 break;
 935         case 1:                         /* yy */
 936                 rr.h = yy.h;
 937                 rr.l = yy.l;
 938                 break;
 939         case 2:                         /* (x+1)*yy */
 940                 z = x + one;            /* may not be exact */
 941                 zh = CHOPPED((z));
 942                 rr.h = zh * yy.h;
 943                 rr.l = z * yy.l + (x - (zh - one)) * yy.h;
 944                 break;
 945         case 3:                         /* (x+2)*(x+1)*yy */
 946                 z1 = x + one;
 947                 z2 = x + 2.0L;
 948                 z = z1 * z2;
 949                 xh = CHOPPED((z));
 950                 zh = CHOPPED((z1));
 951                 xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
 952 
 953                 rr.h = xh * yy.h;


 981                 z *= z2;
 982                 xh = CHOPPED((z));
 983                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 984                 rr.h = xh * yy.h;
 985                 rr.l = z * yy.l + xl * yy.h;
 986                 break;
 987         case 6:                         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
 988                 z1 = x + 2.0L;
 989                 z2 = x + 3.0L;
 990                 z = z1 * z2;
 991                 zh = CHOPPED((z1));
 992                 yh = CHOPPED((z));
 993                 z1 = x - (zh - 2.0L);
 994                 yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
 995                 z2 = z - 2.0L;
 996                 x5 = x + 5.0L;
 997                 z *= z2;
 998                 xh = CHOPPED(z);
 999                 zh += 3.0;
1000                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
1001 
1002                 /*
1003                  * xh+xl=(x+1)*...*(x+4)
1004                  * wh+wl=(x+5)*yy
1005                  */
1006                 wh = CHOPPED((x5 * (yy.h + yy.l)));
1007                 wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
1008                 rr.h = wh * xh;
1009                 rr.l = z * wl + xl * wh;
1010                 break;
1011         case 7:                 /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
1012                 z1 = x + 3.0L;
1013                 z2 = x + 4.0L;
1014                 z = z2 * z1;
1015                 zh = CHOPPED((z1));
1016                 yh = CHOPPED((z));      /* yh+yl = (x+3)(x+4) */
1017                 yl = (x - (zh - 3.0L)) * (z2 + zh) - (yh - (zh * (zh + one)));
1018                 z1 = x + 6.0L;
1019                 z2 = z - 2.0L;          /* z2 = (x+2)*(x+5) */
1020                 z *= z2;
1021                 xh = CHOPPED((z));
1022                 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
1023 
1024                 /*
1025                  * xh+xl=(x+2)*...*(x+5)
1026                  * wh+wl=(x+1)(x+6)*yy
1027                  */
1028                 z2 -= 4.0L;             /* z2 = (x+1)(x+6) */
1029                 wh = CHOPPED((z2 * (yy.h + yy.l)));
1030                 wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0L) * yy.h);
1031                 rr.h = wh * xh;
1032                 rr.l = z * wl + xl * wh;
1033         }
1034 
1035         return (rr);
1036 }
1037 
1038 long double
1039 tgammal(long double x)
1040 {
1041         struct LDouble ss, ww;
1042         long double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
1043         int i, j, m, ix, hx, xk;
1044         unsigned lx;
1045 
1046         hx = H0_WORD(x);
1047         lx = H3_WORD(x);
1048         ix = hx & 0x7fffffff;
1049         y = x;
1050 
1051         if (ix < 0x3f8e0000)         /* x < 2**-113 */
1052                 return (one / x);
1053 
1054         if (ix >= 0x7fff0000)
1055                 return (x * ((hx < 0) ? zero : x));  /* Inf or NaN */
1056 
1057         if (x > overflow)                            /* overflow threshold */
1058                 return (x * 1.0e4932L);
1059 
1060         if (hx >= 0x40020000) {                              /* x >= 8 */
1061                 ww = large_gam(x, &m);
1062                 w = ww.h + ww.l;
1063                 return (scalbnl(w, m));
1064         }
1065 
1066         if (hx > 0) {                        /* 0 < x < 8 */
1067                 i = (int)x;
1068                 ww = gam_n(i, x - (long double)i);
1069                 return (ww.h + ww.l);
1070         }
1071 
1072         /*
1073          * negative x
1074          */
1075 
1076         /*
1077          * compute xk =
1078          *      -2 ... x is an even int (-inf is considered an even #)
1079          *      -1 ... x is an odd int
1080          *      +0 ... x is not an int but chopped to an even int
1081          *      +1 ... x is not an int but chopped to an odd int
1082          */

1083         xk = 0;
1084 #if defined(__x86)
1085         if (ix >= 0x403e0000) {              /* x >= 2**63 } */
1086                 if (ix >= 0x403f0000)
1087                         xk = -2;
1088                 else
1089                         xk = -2 + (lx & 1);
1090 #else
1091         if (ix >= 0x406f0000) {              /* x >= 2**112 */
1092                 if (ix >= 0x40700000)
1093                         xk = -2;
1094                 else
1095                         xk = -2 + (lx & 1);
1096 #endif
1097         } else if (ix >= 0x3fff0000) {
1098                 w = -x;
1099                 t1 = floorl(w);
1100                 t2 = t1 * half;
1101                 t3 = floorl(t2);
1102 
1103                 if (t1 == w) {
1104                         if (t2 == t3)
1105                                 xk = -2;
1106                         else
1107                                 xk = -1;
1108                 } else {
1109                         if (t2 == t3)
1110                                 xk = 0;
1111                         else
1112                                 xk = 1;
1113                 }
1114         }
1115 
1116         if (xk < 0) {
1117                 /* return NaN. Ideally gamma(-n)= (-1)**(n+1) * inf */
1118                 return ((x - x) / (x - x));
1119         }
1120 
1121         /*
1122          * negative underflow thresold -(1774+9ulp)
1123          */
1124         if (x < -1774.0000000000000000000000000000017749370L) {
1125                 z = tiny / x;
1126 
1127                 if (xk == 1)
1128                         z = -z;
1129 
1130                 return (z * tiny);
1131         }
1132 
1133 
1134         /*
1135          * now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
1136          */
1137 
1138         /*
1139          * First compute ss = -sin(pi*y)/pi so that
1140          * gamma(x) = 1/(ss*gamma(1+y))
1141          */

1142         y = -x;
1143         j = (int)y;
1144         z = y - (long double)j;
1145 
1146         if (z > 0.3183098861837906715377675L) {
1147                 if (z > 0.6816901138162093284622325L)
1148                         ss = kpsin(one - z);
1149                 else
1150                         ss = kpcos(0.5L - z);
1151         } else {
1152                 ss = kpsin(z);
1153         }
1154 
1155         if (xk == 0) {
1156                 ss.h = -ss.h;
1157                 ss.l = -ss.l;
1158         }
1159 
1160         /* Then compute ww = gamma(1+y), note that result scale to 2**m */
1161         m = 0;
1162 
1163         if (j < 7) {
1164                 ww = gam_n(j + 1, z);
1165         } else {
1166                 w = y + one;
1167 
1168                 if ((lx & 1) == 0) {        /* y+1 exact (note that y<184) */
1169                         ww = large_gam(w, &m);
1170                 } else {
1171                         t = w - one;
1172 
1173                         if (t == y) {   /* y+one exact */
1174                                 ww = large_gam(w, &m);
1175                         } else {        /* use y*gamma(y) */
1176                                 if (j == 7)
1177                                         ww = gam_n(j, z);
1178                                 else
1179                                         ww = large_gam(y, &m);
1180 
1181                                 t4 = ww.h + ww.l;
1182                                 t1 = CHOPPED((y));
1183                                 t2 = CHOPPED((t4));
1184                                 /* t4 will not be too large */
1185                                 ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
1186                                 ww.h = t1 * t2;
1187                         }
1188                 }
1189         }
1190 
1191         /* compute 1/(ss*ww) */
1192         t3 = ss.h + ss.l;
1193         t4 = ww.h + ww.l;
1194         t1 = CHOPPED((t3));
1195         t2 = CHOPPED((t4));
1196         z1 = ss.l - (t1 - ss.h);        /* (t1,z1) = ss */
1197         z2 = ww.l - (t2 - ww.h);        /* (t2,z2) = ww */
1198         t3 = t3 * t4;                   /* t3 = ss*ww */
1199         z3 = one / t3;                  /* z3 = 1/(ss*ww) */
1200         t5 = t1 * t2;