1755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* @(#)e_lgamma_r.c 1.3 95/01/18 */ 2755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 3755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * ==================================================== 4755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * 6755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Developed at SunSoft, a Sun Microsystems, Inc. business. 7755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Permission to use, copy, modify, and distribute this 8755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * software is freely granted, provided that this notice 9755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * is preserved. 10755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * ==================================================== 11755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 12755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 13755318548d787c7ca36901be064c69efba41c6cfElliott Hughes#include <sys/cdefs.h> 14755318548d787c7ca36901be064c69efba41c6cfElliott Hughes__FBSDID("$FreeBSD$"); 15755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 16755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 17755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * See e_lgamma_r.c for complete comments. 18755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * 19755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Converted to long double by Steven G. Kargl. 20755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 21755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 22755318548d787c7ca36901be064c69efba41c6cfElliott Hughes#include "fpmath.h" 23755318548d787c7ca36901be064c69efba41c6cfElliott Hughes#include "math.h" 24755318548d787c7ca36901be064c69efba41c6cfElliott Hughes#include "math_private.h" 25755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 26755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const volatile double vzero = 0; 27755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 28755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const double 29755318548d787c7ca36901be064c69efba41c6cfElliott Hugheszero= 0, 30755318548d787c7ca36901be064c69efba41c6cfElliott Hugheshalf= 0.5, 31755318548d787c7ca36901be064c69efba41c6cfElliott Hughesone = 1; 32755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 33755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 34755318548d787c7ca36901be064c69efba41c6cfElliott Hughespi = 3.14159265358979323846264338327950288e+00L; 35755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 36755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Domain y in [0x1p-119, 0.28], range ~[-1.4065e-36, 1.4065e-36]: 37755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-119.1 38755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 39755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 40755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa0 = 7.72156649015328606065120900824024296e-02L, 41755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa1 = 3.22467033424113218236207583323018498e-01L, 42755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa2 = 6.73523010531980951332460538330282217e-02L, 43755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa3 = 2.05808084277845478790009252803463129e-02L, 44755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa4 = 7.38555102867398526627292839296001626e-03L, 45755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa5 = 2.89051033074152328576829509522483468e-03L, 46755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa6 = 1.19275391170326097618357349881842913e-03L, 47755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa7 = 5.09669524743042462515256340206203019e-04L, 48755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa8 = 2.23154758453578096143609255559576017e-04L, 49755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa9 = 9.94575127818397632126978731542755129e-05L, 50755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa10 = 4.49262367375420471287545895027098145e-05L, 51755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa11 = 2.05072127845117995426519671481628849e-05L, 52755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa12 = 9.43948816959096748454087141447939513e-06L, 53755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa13 = 4.37486780697359330303852050718287419e-06L, 54755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa14 = 2.03920783892362558276037363847651809e-06L, 55755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa15 = 9.55191070057967287877923073200324649e-07L, 56755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa16 = 4.48993286185740853170657139487620560e-07L, 57755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa17 = 2.13107543597620911675316728179563522e-07L, 58755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa18 = 9.70745379855304499867546549551023473e-08L, 59755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa19 = 5.61889970390290257926487734695402075e-08L, 60755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa20 = 6.42739653024130071866684358960960951e-09L, 61755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa21 = 3.34491062143649291746195612991870119e-08L, 62755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa22 = -1.57068547394315223934653011440641472e-08L, 63755318548d787c7ca36901be064c69efba41c6cfElliott Hughesa23 = 1.30812825422415841213733487745200632e-08L; 64755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 65755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Domain x in [tc-0.24, tc+0.28], range ~[-6.3201e-37, 6.3201e-37]: 66755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * |(lgamma(x) - tf) - t(x - tc)| < 2**-120.3. 67755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 68755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 69755318548d787c7ca36901be064c69efba41c6cfElliott Hughestc = 1.46163214496836234126265954232572133e+00L, 70755318548d787c7ca36901be064c69efba41c6cfElliott Hughestf = -1.21486290535849608095514557177691584e-01L, 71755318548d787c7ca36901be064c69efba41c6cfElliott Hughestt = 1.57061739945077675484237837992951704e-36L, 72755318548d787c7ca36901be064c69efba41c6cfElliott Hughest0 = -1.99238329499314692728655623767019240e-36L, 73755318548d787c7ca36901be064c69efba41c6cfElliott Hughest1 = -6.08453430711711404116887457663281416e-35L, 74755318548d787c7ca36901be064c69efba41c6cfElliott Hughest2 = 4.83836122723810585213722380854828904e-01L, 75755318548d787c7ca36901be064c69efba41c6cfElliott Hughest3 = -1.47587722994530702030955093950668275e-01L, 76755318548d787c7ca36901be064c69efba41c6cfElliott Hughest4 = 6.46249402389127526561003464202671923e-02L, 77755318548d787c7ca36901be064c69efba41c6cfElliott Hughest5 = -3.27885410884813055008502586863748063e-02L, 78755318548d787c7ca36901be064c69efba41c6cfElliott Hughest6 = 1.79706751152103942928638276067164935e-02L, 79755318548d787c7ca36901be064c69efba41c6cfElliott Hughest7 = -1.03142230366363872751602029672767978e-02L, 80755318548d787c7ca36901be064c69efba41c6cfElliott Hughest8 = 6.10053602051788840313573150785080958e-03L, 81755318548d787c7ca36901be064c69efba41c6cfElliott Hughest9 = -3.68456960831637325470641021892968954e-03L, 82755318548d787c7ca36901be064c69efba41c6cfElliott Hughest10 = 2.25976482322181046611440855340968560e-03L, 83755318548d787c7ca36901be064c69efba41c6cfElliott Hughest11 = -1.40225144590445082933490395950664961e-03L, 84755318548d787c7ca36901be064c69efba41c6cfElliott Hughest12 = 8.78232634717681264035014878172485575e-04L, 85755318548d787c7ca36901be064c69efba41c6cfElliott Hughest13 = -5.54194952796682301220684760591403899e-04L, 86755318548d787c7ca36901be064c69efba41c6cfElliott Hughest14 = 3.51912956837848209220421213975000298e-04L, 87755318548d787c7ca36901be064c69efba41c6cfElliott Hughest15 = -2.24653443695947456542669289367055542e-04L, 88755318548d787c7ca36901be064c69efba41c6cfElliott Hughest16 = 1.44070395420840737695611929680511823e-04L, 89755318548d787c7ca36901be064c69efba41c6cfElliott Hughest17 = -9.27609865550394140067059487518862512e-05L, 90755318548d787c7ca36901be064c69efba41c6cfElliott Hughest18 = 5.99347334438437081412945428365433073e-05L, 91755318548d787c7ca36901be064c69efba41c6cfElliott Hughest19 = -3.88458388854572825603964274134801009e-05L, 92755318548d787c7ca36901be064c69efba41c6cfElliott Hughest20 = 2.52476631610328129217896436186551043e-05L, 93755318548d787c7ca36901be064c69efba41c6cfElliott Hughest21 = -1.64508584981658692556994212457518536e-05L, 94755318548d787c7ca36901be064c69efba41c6cfElliott Hughest22 = 1.07434583475987007495523340296173839e-05L, 95755318548d787c7ca36901be064c69efba41c6cfElliott Hughest23 = -7.03070407519397260929482550448878399e-06L, 96755318548d787c7ca36901be064c69efba41c6cfElliott Hughest24 = 4.60968590693753579648385629003100469e-06L, 97755318548d787c7ca36901be064c69efba41c6cfElliott Hughest25 = -3.02765473778832036018438676945512661e-06L, 98755318548d787c7ca36901be064c69efba41c6cfElliott Hughest26 = 1.99238771545503819972741288511303401e-06L, 99755318548d787c7ca36901be064c69efba41c6cfElliott Hughest27 = -1.31281299822614084861868817951788579e-06L, 100755318548d787c7ca36901be064c69efba41c6cfElliott Hughest28 = 8.60844432267399655055574642052370223e-07L, 101755318548d787c7ca36901be064c69efba41c6cfElliott Hughest29 = -5.64535486432397413273248363550536374e-07L, 102755318548d787c7ca36901be064c69efba41c6cfElliott Hughest30 = 3.99357783676275660934903139592727737e-07L, 103755318548d787c7ca36901be064c69efba41c6cfElliott Hughest31 = -2.95849029193433121795495215869311610e-07L, 104755318548d787c7ca36901be064c69efba41c6cfElliott Hughest32 = 1.37790144435073124976696250804940384e-07L; 105755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 106755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Domain y in [-0.1, 0.232], range ~[-1.4046e-37, 1.4181e-37]: 107755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-122.8 108755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 109755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 110755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu0 = -7.72156649015328606065120900824024311e-02L, 111755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu1 = 4.24082772271938167430983113242482656e-01L, 112755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu2 = 2.96194003481457101058321977413332171e+00L, 113755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu3 = 6.49503267711258043997790983071543710e+00L, 114755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu4 = 7.40090051288150177152835698948644483e+00L, 115755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu5 = 4.94698036296756044610805900340723464e+00L, 116755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu6 = 2.00194224610796294762469550684947768e+00L, 117755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu7 = 4.82073087750608895996915051568834949e-01L, 118755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu8 = 6.46694052280506568192333848437585427e-02L, 119755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu9 = 4.17685526755100259316625348933108810e-03L, 120755318548d787c7ca36901be064c69efba41c6cfElliott Hughesu10 = 9.06361003550314327144119307810053410e-05L, 121755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv1 = 5.15937098592887275994320496999951947e+00L, 122755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv2 = 1.14068418766251486777604403304717558e+01L, 123755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv3 = 1.41164839437524744055723871839748489e+01L, 124755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv4 = 1.07170702656179582805791063277960532e+01L, 125755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv5 = 5.14448694179047879915042998453632434e+00L, 126755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv6 = 1.55210088094585540637493826431170289e+00L, 127755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv7 = 2.82975732849424562719893657416365673e-01L, 128755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv8 = 2.86424622754753198010525786005443539e-02L, 129755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv9 = 1.35364253570403771005922441442688978e-03L, 130755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv10 = 1.91514173702398375346658943749580666e-05L, 131755318548d787c7ca36901be064c69efba41c6cfElliott Hughesv11 = -3.25364686890242327944584691466034268e-08L; 132755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 133755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Domain x in (2, 3], range ~[-1.3341e-36, 1.3536e-36]: 134755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-120.1 135755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * with y = x - 2. 136755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 137755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 138755318548d787c7ca36901be064c69efba41c6cfElliott Hughess0 = -7.72156649015328606065120900824024297e-02L, 139755318548d787c7ca36901be064c69efba41c6cfElliott Hughess1 = 1.23221687850916448903914170805852253e-01L, 140755318548d787c7ca36901be064c69efba41c6cfElliott Hughess2 = 5.43673188699937239808255378293820020e-01L, 141755318548d787c7ca36901be064c69efba41c6cfElliott Hughess3 = 6.31998137119005233383666791176301800e-01L, 142755318548d787c7ca36901be064c69efba41c6cfElliott Hughess4 = 3.75885340179479850993811501596213763e-01L, 143755318548d787c7ca36901be064c69efba41c6cfElliott Hughess5 = 1.31572908743275052623410195011261575e-01L, 144755318548d787c7ca36901be064c69efba41c6cfElliott Hughess6 = 2.82528453299138685507186287149699749e-02L, 145755318548d787c7ca36901be064c69efba41c6cfElliott Hughess7 = 3.70262021550340817867688714880797019e-03L, 146755318548d787c7ca36901be064c69efba41c6cfElliott Hughess8 = 2.83374000312371199625774129290973648e-04L, 147755318548d787c7ca36901be064c69efba41c6cfElliott Hughess9 = 1.15091830239148290758883505582343691e-05L, 148755318548d787c7ca36901be064c69efba41c6cfElliott Hughess10 = 2.04203474281493971326506384646692446e-07L, 149755318548d787c7ca36901be064c69efba41c6cfElliott Hughess11 = 9.79544198078992058548607407635645763e-10L, 150755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr1 = 2.58037466655605285937112832039537492e+00L, 151755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr2 = 2.86289413392776399262513849911531180e+00L, 152755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr3 = 1.78691044735267497452847829579514367e+00L, 153755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr4 = 6.89400381446725342846854215600008055e-01L, 154755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr5 = 1.70135865462567955867134197595365343e-01L, 155755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr6 = 2.68794816183964420375498986152766763e-02L, 156755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr7 = 2.64617234244861832870088893332006679e-03L, 157755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr8 = 1.52881761239180800640068128681725702e-04L, 158755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr9 = 4.63264813762296029824851351257638558e-06L, 159755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr10 = 5.89461519146957343083848967333671142e-08L, 160755318548d787c7ca36901be064c69efba41c6cfElliott Hughesr11 = 1.79027678176582527798327441636552968e-10L; 161755318548d787c7ca36901be064c69efba41c6cfElliott Hughes/* 162755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * Domain z in [8, 0x1p70], range ~[-9.8214e-35, 9.8214e-35]: 163755318548d787c7ca36901be064c69efba41c6cfElliott Hughes * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-113.0 164755318548d787c7ca36901be064c69efba41c6cfElliott Hughes */ 165755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic const long double 166755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw0 = 4.18938533204672741780329736405617738e-01L, 167755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw1 = 8.33333333333333333333333333332852026e-02L, 168755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw2 = -2.77777777777777777777777727810123528e-03L, 169755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw3 = 7.93650793650793650791708939493907380e-04L, 170755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw4 = -5.95238095238095234390450004444370959e-04L, 171755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw5 = 8.41750841750837633887817658848845695e-04L, 172755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw6 = -1.91752691752396849943172337347259743e-03L, 173755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw7 = 6.41025640880333069429106541459015557e-03L, 174755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw8 = -2.95506530801732133437990433080327074e-02L, 175755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw9 = 1.79644237328444101596766586979576927e-01L, 176755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw10 = -1.39240539108367641920172649259736394e+00L, 177755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw11 = 1.33987701479007233325288857758641761e+01L, 178755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw12 = -1.56363596431084279780966590116006255e+02L, 179755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw13 = 2.14830978044410267201172332952040777e+03L, 180755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw14 = -3.28636067474227378352761516589092334e+04L, 181755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw15 = 5.06201257747865138432663574251462485e+05L, 182755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw16 = -6.79720123352023636706247599728048344e+06L, 183755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw17 = 6.57556601705472106989497289465949255e+07L, 184755318548d787c7ca36901be064c69efba41c6cfElliott Hughesw18 = -3.26229058141181783534257632389415580e+08L; 185755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 186755318548d787c7ca36901be064c69efba41c6cfElliott Hughesstatic long double 187755318548d787c7ca36901be064c69efba41c6cfElliott Hughessin_pil(long double x) 188755318548d787c7ca36901be064c69efba41c6cfElliott Hughes{ 189755318548d787c7ca36901be064c69efba41c6cfElliott Hughes volatile long double vz; 190755318548d787c7ca36901be064c69efba41c6cfElliott Hughes long double y,z; 191755318548d787c7ca36901be064c69efba41c6cfElliott Hughes uint64_t lx, n; 192755318548d787c7ca36901be064c69efba41c6cfElliott Hughes uint16_t hx; 193755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 194755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y = -x; 195755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 196755318548d787c7ca36901be064c69efba41c6cfElliott Hughes vz = y+0x1.p112; 197755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z = vz-0x1.p112; 198755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if (z == y) 199755318548d787c7ca36901be064c69efba41c6cfElliott Hughes return zero; 200755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 201755318548d787c7ca36901be064c69efba41c6cfElliott Hughes vz = y+0x1.p110; 202755318548d787c7ca36901be064c69efba41c6cfElliott Hughes EXTRACT_LDBL128_WORDS(hx,lx,n,vz); 203755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z = vz-0x1.p110; 204755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if (z > y) { 205755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z -= 0.25; 206755318548d787c7ca36901be064c69efba41c6cfElliott Hughes n--; 207755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 208755318548d787c7ca36901be064c69efba41c6cfElliott Hughes n &= 7; 209488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes y = y - z + n * 0.25; 210755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 211755318548d787c7ca36901be064c69efba41c6cfElliott Hughes switch (n) { 212755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 0: y = __kernel_sinl(pi*y,zero,0); break; 213755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 1: 214755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; 215488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes case 3: 216755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; 217755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 5: 218755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; 219755318548d787c7ca36901be064c69efba41c6cfElliott Hughes default: y = __kernel_sinl(pi*(y-2.0),zero,0); break; 220755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 221755318548d787c7ca36901be064c69efba41c6cfElliott Hughes return -y; 222755318548d787c7ca36901be064c69efba41c6cfElliott Hughes} 223755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 224755318548d787c7ca36901be064c69efba41c6cfElliott Hugheslong double 225755318548d787c7ca36901be064c69efba41c6cfElliott Hugheslgammal_r(long double x, int *signgamp) 226755318548d787c7ca36901be064c69efba41c6cfElliott Hughes{ 227755318548d787c7ca36901be064c69efba41c6cfElliott Hughes long double nadj,p,p1,p2,p3,q,r,t,w,y,z; 228755318548d787c7ca36901be064c69efba41c6cfElliott Hughes uint64_t llx,lx; 229755318548d787c7ca36901be064c69efba41c6cfElliott Hughes int i; 230488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes uint16_t hx,ix; 231755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 232488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes EXTRACT_LDBL128_WORDS(hx,lx,llx,x); 233755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 234488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* purge +-Inf and NaNs */ 235755318548d787c7ca36901be064c69efba41c6cfElliott Hughes *signgamp = 1; 236488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes ix = hx&0x7fff; 237488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(ix==0x7fff) return x*x; 238755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 239488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* purge +-0 and tiny arguments */ 240488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes *signgamp = 1-2*(hx>>15); 241488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(ix<0x3fff-116) { /* |x|<2**-(p+3), return -log(|x|) */ 242488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if((ix|lx|llx)==0) 243488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes return one/vzero; 244488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes return -logl(fabsl(x)); 245755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 246488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes 247488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* purge negative integers and start evaluation for other x < 0 */ 248755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if(hx&0x8000) { 249488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes *signgamp = 1; 250488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */ 251755318548d787c7ca36901be064c69efba41c6cfElliott Hughes return one/vzero; 252755318548d787c7ca36901be064c69efba41c6cfElliott Hughes t = sin_pil(x); 253755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if(t==zero) return one/vzero; 254755318548d787c7ca36901be064c69efba41c6cfElliott Hughes nadj = logl(pi/fabsl(t*x)); 255755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if(t<zero) *signgamp = -1; 256755318548d787c7ca36901be064c69efba41c6cfElliott Hughes x = -x; 257755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 258755318548d787c7ca36901be064c69efba41c6cfElliott Hughes 259488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* purge 1 and 2 */ 260488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0; 261488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* for x < 2.0 */ 262488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes else if(ix<0x4000) { 263488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(x<=8.9999961853027344e-01) { 264755318548d787c7ca36901be064c69efba41c6cfElliott Hughes r = -logl(x); 265488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} 266488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} 267755318548d787c7ca36901be064c69efba41c6cfElliott Hughes else {y = x; i=2;} 268755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } else { 269488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes r = 0; 270488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes if(x>=1.7316312789916992e+00) {y=2-x;i=0;} 271488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} 272755318548d787c7ca36901be064c69efba41c6cfElliott Hughes else {y=x-1;i=2;} 273755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 274755318548d787c7ca36901be064c69efba41c6cfElliott Hughes switch(i) { 275755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 0: 276755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z = y*y; 277755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+ 278755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z*(a18+z*(a20+z*a22)))))))))); 279755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ 280755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z*(a17+z*(a19+z*(a21+z*a23))))))))))); 281755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p = y*p1+p2; 282488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes r += p-y/2; break; 283755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 1: 284755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ 285755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ 286755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+ 287755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+ 288755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(t31+y*t32)))))))))))))))))))))))))))))); 289488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes r += tf + p; break; 290755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 2: 291755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+ 292755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(u8+y*(u9+y*u10)))))))))); 293755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+ 294755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(v8+y*(v9+y*(v10+y*v11)))))))))); 295488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes r += p1/p2-y/2; 296755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 297755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 298488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* x < 8.0 */ 299488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes else if(ix<0x4002) { 300755318548d787c7ca36901be064c69efba41c6cfElliott Hughes i = x; 301755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y = x-i; 302755318548d787c7ca36901be064c69efba41c6cfElliott Hughes p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ 303755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(s9+y*(s10+y*s11))))))))))); 304755318548d787c7ca36901be064c69efba41c6cfElliott Hughes q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+ 305755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(r9+y*(r10+y*r11)))))))))); 306755318548d787c7ca36901be064c69efba41c6cfElliott Hughes r = y/2+p/q; 307755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ 308755318548d787c7ca36901be064c69efba41c6cfElliott Hughes switch(i) { 309755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 7: z *= (y+6); /* FALLTHRU */ 310755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 6: z *= (y+5); /* FALLTHRU */ 311755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 5: z *= (y+4); /* FALLTHRU */ 312755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 4: z *= (y+3); /* FALLTHRU */ 313755318548d787c7ca36901be064c69efba41c6cfElliott Hughes case 3: z *= (y+2); /* FALLTHRU */ 314755318548d787c7ca36901be064c69efba41c6cfElliott Hughes r += logl(z); break; 315755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } 316488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* 8.0 <= x < 2**(p+3) */ 317488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes } else if (ix<0x3fff+116) { 318755318548d787c7ca36901be064c69efba41c6cfElliott Hughes t = logl(x); 319755318548d787c7ca36901be064c69efba41c6cfElliott Hughes z = one/x; 320755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y = z*z; 321755318548d787c7ca36901be064c69efba41c6cfElliott Hughes w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*(w8+ 322755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+ 323755318548d787c7ca36901be064c69efba41c6cfElliott Hughes y*(w17+y*w18))))))))))))))))); 324755318548d787c7ca36901be064c69efba41c6cfElliott Hughes r = (x-half)*(t-one)+w; 325488268b134723c7a6598338bb253be5f64d53be4Elliott Hughes /* 2**(p+3) <= x <= inf */ 326755318548d787c7ca36901be064c69efba41c6cfElliott Hughes } else 327755318548d787c7ca36901be064c69efba41c6cfElliott Hughes r = x*(logl(x)-1); 328755318548d787c7ca36901be064c69efba41c6cfElliott Hughes if(hx&0x8000) r = nadj - r; 329755318548d787c7ca36901be064c69efba41c6cfElliott Hughes return r; 330755318548d787c7ca36901be064c69efba41c6cfElliott Hughes} 331