mirrors/cosmopolitan
1
0
Fork 0
mirror of https://github.com/jart/cosmopolitan.git synced 2025-03-03 07:29:23 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Introduce lgammal, tgammal, erfl, and erfcl

Browse source 
git://git.musl-libc.org/musl
79bdacff83a6bd5b70ff5ae5eb8b6de82c2f7c30
This commit is contained in:
Justine Tunney 2023-10-31 21:58:29 -07:00
parent b4084dd6c6
commit 3b791d2f44
No known key found for this signature in database
GPG key ID: BE714B4575D6E328
8 changed files with 1100 additions and 319 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
libc/tinymath/erf.c
View file

@ -334,3 +334,8 @@ double erfc(double x)
}
return sign ? 2 - 0x1p-1022 : 0x1p-1022*0x1p-1022;
}
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
__weak_reference(erf, erfl);
__weak_reference(erfc, erfcl);
#endif

381
libc/tinymath/erfl.c Normal file
View file

@ -0,0 +1,381 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│
╚──────────────────────────────────────────────────────────────────────────────╝
│ │
│ Musl Libc │
│ Copyright © 2005-2014 Rich Felker, et al. │
│ │
│ Permission is hereby granted, free of charge, to any person obtaining │
│ a copy of this software and associated documentation files (the │
│ "Software"), to deal in the Software without restriction, including │
│ without limitation the rights to use, copy, modify, merge, publish, │
│ distribute, sublicense, and/or sell copies of the Software, and to │
│ permit persons to whom the Software is furnished to do so, subject to │
│ the following conditions: │
│ │
│ The above copyright notice and this permission notice shall be │
│ included in all copies or substantial portions of the Software. │
│ │
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
│ │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
#include "libc/tinymath/ldshape.internal.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
asm(".ident\t\"\\n\\n\
OpenBSD libm (ISC License)\\n\
Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_erfl.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. For |x| in [0, 0.84375]
* erf(x) = x + x*R(x^2)
* erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
* = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
* guarantee the error is less than one ulp for erf.
*
* 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(x) = sign(x) * (c + P1(s)/Q1(s))
* erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
*
* 3. For x in [1.25,1/0.35(~2.857143)],
* erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
* z=1/x^2
* erf(x) = 1 - erfc(x)
*
* 4. For x in [1/0.35,107]
* erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
* if -6.666<x<0
* = 2.0 - tiny (if x <= -6.666)
* z=1/x^2
* erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
* erf(x) = sign(x)*(1.0 - tiny)
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
*
* 5. For inf > x >= 107
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
static const long double
erx = 0.845062911510467529296875L,
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
/* 8 * (2/sqrt(pi) - 1) */
efx8 = 1.0270333367641005911692712249723613735048E0L,
pp[6] = {
1.122751350964552113068262337278335028553E6L,
-2.808533301997696164408397079650699163276E6L,
-3.314325479115357458197119660818768924100E5L,
-6.848684465326256109712135497895525446398E4L,
-2.657817695110739185591505062971929859314E3L,
-1.655310302737837556654146291646499062882E2L,
},
qq[6] = {
8.745588372054466262548908189000448124232E6L,
3.746038264792471129367533128637019611485E6L,
7.066358783162407559861156173539693900031E5L,
7.448928604824620999413120955705448117056E4L,
4.511583986730994111992253980546131408924E3L,
1.368902937933296323345610240009071254014E2L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
/* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
-0.15625 <= x <= +.25
Peak relative error 8.5e-22 */
pa[8] = {
-1.076952146179812072156734957705102256059E0L,
1.884814957770385593365179835059971587220E2L,
-5.339153975012804282890066622962070115606E1L,
4.435910679869176625928504532109635632618E1L,
1.683219516032328828278557309642929135179E1L,
-2.360236618396952560064259585299045804293E0L,
1.852230047861891953244413872297940938041E0L,
9.394994446747752308256773044667843200719E-2L,
},
qa[7] = {
4.559263722294508998149925774781887811255E2L,
3.289248982200800575749795055149780689738E2L,
2.846070965875643009598627918383314457912E2L,
1.398715859064535039433275722017479994465E2L,
6.060190733759793706299079050985358190726E1L,
2.078695677795422351040502569964299664233E1L,
4.641271134150895940966798357442234498546E0L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
1/2.85711669921875 < 1/x < 1/1.25
Peak relative error 3.1e-21 */
ra[] = {
1.363566591833846324191000679620738857234E-1L,
1.018203167219873573808450274314658434507E1L,
1.862359362334248675526472871224778045594E2L,
1.411622588180721285284945138667933330348E3L,
5.088538459741511988784440103218342840478E3L,
8.928251553922176506858267311750789273656E3L,
7.264436000148052545243018622742770549982E3L,
2.387492459664548651671894725748959751119E3L,
2.220916652813908085449221282808458466556E2L,
},
sa[] = {
-1.382234625202480685182526402169222331847E1L,
-3.315638835627950255832519203687435946482E2L,
-2.949124863912936259747237164260785326692E3L,
-1.246622099070875940506391433635999693661E4L,
-2.673079795851665428695842853070996219632E4L,
-2.880269786660559337358397106518918220991E4L,
-1.450600228493968044773354186390390823713E4L,
-2.874539731125893533960680525192064277816E3L,
-1.402241261419067750237395034116942296027E2L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erfc in [1/.35,107]
*/
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
1/6.6666259765625 < 1/x < 1/2.85711669921875
Peak relative error 4.2e-22 */
rb[] = {
-4.869587348270494309550558460786501252369E-5L,
-4.030199390527997378549161722412466959403E-3L,
-9.434425866377037610206443566288917589122E-2L,
-9.319032754357658601200655161585539404155E-1L,
-4.273788174307459947350256581445442062291E0L,
-8.842289940696150508373541814064198259278E0L,
-7.069215249419887403187988144752613025255E0L,
-1.401228723639514787920274427443330704764E0L,
},
sb[] = {
4.936254964107175160157544545879293019085E-3L,
1.583457624037795744377163924895349412015E-1L,
1.850647991850328356622940552450636420484E0L,
9.927611557279019463768050710008450625415E0L,
2.531667257649436709617165336779212114570E1L,
2.869752886406743386458304052862814690045E1L,
1.182059497870819562441683560749192539345E1L,
/* 1.000000000000000000000000000000000000000E0 */
},
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
1/107 <= 1/x <= 1/6.6666259765625
Peak relative error 1.1e-21 */
rc[] = {
-8.299617545269701963973537248996670806850E-5L,
-6.243845685115818513578933902532056244108E-3L,
-1.141667210620380223113693474478394397230E-1L,
-7.521343797212024245375240432734425789409E-1L,
-1.765321928311155824664963633786967602934E0L,
-1.029403473103215800456761180695263439188E0L,
},
sc[] = {
8.413244363014929493035952542677768808601E-3L,
2.065114333816877479753334599639158060979E-1L,
1.639064941530797583766364412782135680148E0L,
4.936788463787115555582319302981666347450E0L,
5.005177727208955487404729933261347679090E0L,
/* 1.000000000000000000000000000000000000000E0 */
};
static long double erfc1(long double x)
{
long double s,P,Q;
s = fabsl(x) - 1;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
return 1 - erx - P / Q;
}
static long double erfc2(uint32_t ix, long double x)
{
union ldshape u;
long double s,z,R,S;
if (ix < 0x3fffa000) /* 0.84375 <= |x| < 1.25 */
return erfc1(x);
x = fabsl(x);
s = 1 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.857 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
} else if (ix < 0x4001d555) { /* 2.857 <= |x| < 6.6666259765625 */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 6.666 <= |x| < 107 (erfc only) */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
u.f = x;
u.i.m &= -1ULL << 40;
z = u.f;
return expl(-z*z - 0.5625) * expl((z - x) * (z + x) + R / S) / x;
}
long double erfl(long double x)
{
long double r, s, z, y;
union ldshape u = {x};
uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
int sign = u.i.se >> 15;
if (ix >= 0x7fff0000)
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1 - 2*sign + 1/x;
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fde8000) { /* |x| < 2**-33 */
return 0.125 * (8 * x + efx8 * x); /* avoid underflow */
}
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1] +
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
return x + x * y;
}
if (ix < 0x4001d555) /* |x| < 6.6666259765625 */
y = 1 - erfc2(ix,x);
else
y = 1 - 0x1p-16382L;
return sign ? -y : y;
}
long double erfcl(long double x)
{
long double r, s, z, y;
union ldshape u = {x};
uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
int sign = u.i.se >> 15;
if (ix >= 0x7fff0000)
/* erfc(nan) = nan, erfc(+-inf) = 0,2 */
return 2*sign + 1/x;
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fbe0000) /* |x| < 2**-65 */
return 1.0 - x;
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1] +
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */
return 1.0 - (x + x * y);
return 0.5 - (x - 0.5 + x * y);
}
if (ix < 0x4005d600) /* |x| < 107 */
return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
y = 0x1p-16382L;
return sign ? 2 - y : y*y;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
// TODO: broken implementation to make things compile
long double erfl(long double x)
{
return erf(x);
}
long double erfcl(long double x)
{
return erfc(x);
}
#endif /* long double is long */

318
libc/tinymath/gamma.c
View file

@ -1,318 +0,0 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│
╚──────────────────────────────────────────────────────────────────────────────╝
│ │
│ Musl Libc │
│ Copyright © 2005-2014 Rich Felker, et al. │
│ │
│ Permission is hereby granted, free of charge, to any person obtaining │
│ a copy of this software and associated documentation files (the │
│ "Software"), to deal in the Software without restriction, including │
│ without limitation the rights to use, copy, modify, merge, publish, │
│ distribute, sublicense, and/or sell copies of the Software, and to │
│ permit persons to whom the Software is furnished to do so, subject to │
│ the following conditions: │
│ │
│ The above copyright notice and this permission notice shall be │
│ included in all copies or substantial portions of the Software. │
│ │
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
│ │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
#include "libc/tinymath/kernel.internal.h"
asm(".ident\t\"\\n\\n\
fdlibm (fdlibm license)\\n\
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
/* clang-format off */
/* origin: FreeBSD /usr/src/lib/msun/src/e_lgamma_r.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
/* lgamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3)
* = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
* 2. Polynomial approximation of lgamma around its
* minimum ymin=1.461632144968362245 to maintain monotonicity.
* On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
* Let z = x-ymin;
* lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
* where
* poly(z) is a 14 degree polynomial.
* 2. Rational approximation in the primary interval [2,3]
* We use the following approximation:
* s = x-2.0;
* lgamma(x) = 0.5*s + s*P(s)/Q(s)
* with accuracy
* |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
* Our algorithms are based on the following observation
*
* zeta(2)-1 2 zeta(3)-1 3
* lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
* 2 3
*
* where Euler = 0.5771... is the Euler constant, which is very
* close to 0.5.
*
* 3. For x>=8, we have
* lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
* (better formula:
* lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
* Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
* where
* |w - f(z)| < 2**-58.74
*
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
* Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
*
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1) = lgamma(2) = 0
* lgamma(x) ~ -log(|x|) for tiny x
* lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
* lgamma(inf) = inf
* lgamma(-inf) = inf (bug for bug compatible with C99!?)
*
*/
static const double
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
/* tt = -(tail of tf) */
tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
/* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */
static double sin_pi(double x)
{
int n;
/* spurious inexact if odd int */
x = 2.0*(x*0.5 - floor(x*0.5)); /* x mod 2.0 */
n = (int)(x*4.0);
n = (n+1)/2;
x -= n*0.5f;
x *= pi;
switch (n) {
default: /* case 4: */
case 0: return __sin(x, 0.0, 0);
case 1: return __cos(x, 0.0);
case 2: return __sin(-x, 0.0, 0);
case 3: return -__cos(x, 0.0);
}
}
double lgamma_r(double x, int *signgamp)
{
union {double f; uint64_t i;} u = {x};
double_t t,y,z,nadj,p,p1,p2,p3,q,r,w;
uint32_t ix;
int sign,i;
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
*signgamp = 1;
sign = u.i>>63;
ix = u.i>>32 & 0x7fffffff;
if (ix >= 0x7ff00000)
return x*x;
if (ix < (0x3ff-70)<<20) { /* |x|<2**-70, return -log(|x|) */
if(sign) {
x = -x;
*signgamp = -1;
}
return -log(x);
}
if (sign) {
x = -x;
t = sin_pi(x);
if (t == 0.0) /* -integer */
return 1.0/(x-x);
if (t > 0.0)
*signgamp = -1;
else
t = -t;
nadj = log(pi/(t*x));
}
/* purge off 1 and 2 */
if ((ix == 0x3ff00000 || ix == 0x40000000) && (uint32_t)u.i == 0)
r = 0;
/* for x < 2.0 */
else if (ix < 0x40000000) {
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x);
if (ix >= 0x3FE76944) {
y = 1.0 - x;
i = 0;
} else if (ix >= 0x3FCDA661) {
y = x - (tc-1.0);
i = 1;
} else {
y = x;
i = 2;
}
} else {
r = 0.0;
if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */
y = 2.0 - x;
i = 0;
} else if(ix >= 0x3FF3B4C4) { /* [1.23,1.73] */
y = x - tc;
i = 1;
} else {
y = x - 1.0;
i = 2;
}
}
switch (i) {
case 0:
z = y*y;
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
p = y*p1+p2;
r += (p-0.5*y);
break;
case 1:
z = y*y;
w = z*y;
p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
p = z*p1-(tt-w*(p2+y*p3));
r += tf + p;
break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5*y + p1/p2;
}
} else if (ix < 0x40200000) { /* x < 8.0 */
i = (int)x;
y = x - (double)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = 0.5*y+p/q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7: z *= y + 6.0; /* FALLTHRU */
case 6: z *= y + 5.0; /* FALLTHRU */
case 5: z *= y + 4.0; /* FALLTHRU */
case 4: z *= y + 3.0; /* FALLTHRU */
case 3: z *= y + 2.0; /* FALLTHRU */
r += log(z);
break;
}
} else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */
t = log(x);
z = 1.0/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-0.5)*(t-1.0)+w;
} else /* 2**58 <= x <= inf */
r = x*(log(x)-1.0);
if (sign)
r = nadj - r;
return r;
}

5
libc/tinymath/lgamma.c
View file

@ -17,7 +17,6 @@
│ PERFORMANCE OF THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
#include "libc/tinymath/kernel.internal.h"
/**
* Returns natural logarithm of absolute value of gamma function.
@ -25,3 +24,7 @@
double lgamma(double x) {
return lgamma_r(x, &signgam);
}
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
__weak_reference(lgamma, lgammal);
#endif

4
libc/tinymath/lgamma_r.c
View file

@ -319,3 +319,7 @@ double lgamma_r(double x, int *signgamp)
r = nadj - r;
return r;
}
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
__weak_reference(lgamma_r, lgammal_r);
#endif

388
libc/tinymath/lgammal.c Normal file
View file

@ -0,0 +1,388 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│
╚──────────────────────────────────────────────────────────────────────────────╝
│ │
│ Musl Libc │
│ Copyright © 2005-2014 Rich Felker, et al. │
│ │
│ Permission is hereby granted, free of charge, to any person obtaining │
│ a copy of this software and associated documentation files (the │
│ "Software"), to deal in the Software without restriction, including │
│ without limitation the rights to use, copy, modify, merge, publish, │
│ distribute, sublicense, and/or sell copies of the Software, and to │
│ permit persons to whom the Software is furnished to do so, subject to │
│ the following conditions: │
│ │
│ The above copyright notice and this permission notice shall be │
│ included in all copies or substantial portions of the Software. │
│ │
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
│ │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
#include "libc/tinymath/kernel.internal.h"
#include "libc/tinymath/ldshape.internal.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
asm(".ident\t\"\\n\\n\
OpenBSD libm (ISC License)\\n\
Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_lgammal.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* lgammal(x)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3)
* = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
* 2. Polynomial approximation of lgamma around its
* minimun ymin=1.461632144968362245 to maintain monotonicity.
* On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
* Let z = x-ymin;
* lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
* 2. Rational approximation in the primary interval [2,3]
* We use the following approximation:
* s = x-2.0;
* lgamma(x) = 0.5*s + s*P(s)/Q(s)
* Our algorithms are based on the following observation
*
* zeta(2)-1 2 zeta(3)-1 3
* lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
* 2 3
*
* where Euler = 0.5771... is the Euler constant, which is very
* close to 0.5.
*
* 3. For x>=8, we have
* lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
* (better formula:
* lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
* Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
*
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
* Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
*
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1)=lgamma(2)=0
* lgamma(x) ~ -log(x) for tiny x
* lgamma(0) = lgamma(inf) = inf
* lgamma(-integer) = +-inf
*
*/
static const long double
pi = 3.14159265358979323846264L,
/* lgam(1+x) = 0.5 x + x a(x)/b(x)
-0.268402099609375 <= x <= 0
peak relative error 6.6e-22 */
a0 = -6.343246574721079391729402781192128239938E2L,
a1 = 1.856560238672465796768677717168371401378E3L,
a2 = 2.404733102163746263689288466865843408429E3L,
a3 = 8.804188795790383497379532868917517596322E2L,
a4 = 1.135361354097447729740103745999661157426E2L,
a5 = 3.766956539107615557608581581190400021285E0L,
b0 = 8.214973713960928795704317259806842490498E3L,
b1 = 1.026343508841367384879065363925870888012E4L,
b2 = 4.553337477045763320522762343132210919277E3L,
b3 = 8.506975785032585797446253359230031874803E2L,
b4 = 6.042447899703295436820744186992189445813E1L,
/* b5 = 1.000000000000000000000000000000000000000E0 */
tc = 1.4616321449683623412626595423257213284682E0L,
tf = -1.2148629053584961146050602565082954242826E-1, /* double precision */
/* tt = (tail of tf), i.e. tf + tt has extended precision. */
tt = 3.3649914684731379602768989080467587736363E-18L,
/* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
-1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
/* lgam (x + tc) = tf + tt + x g(x)/h(x)
-0.230003726999612341262659542325721328468 <= x
<= 0.2699962730003876587373404576742786715318
peak relative error 2.1e-21 */
g0 = 3.645529916721223331888305293534095553827E-18L,
g1 = 5.126654642791082497002594216163574795690E3L,
g2 = 8.828603575854624811911631336122070070327E3L,
g3 = 5.464186426932117031234820886525701595203E3L,
g4 = 1.455427403530884193180776558102868592293E3L,
g5 = 1.541735456969245924860307497029155838446E2L,
g6 = 4.335498275274822298341872707453445815118E0L,
h0 = 1.059584930106085509696730443974495979641E4L,
h1 = 2.147921653490043010629481226937850618860E4L,
h2 = 1.643014770044524804175197151958100656728E4L,
h3 = 5.869021995186925517228323497501767586078E3L,
h4 = 9.764244777714344488787381271643502742293E2L,
h5 = 6.442485441570592541741092969581997002349E1L,
/* h6 = 1.000000000000000000000000000000000000000E0 */
/* lgam (x+1) = -0.5 x + x u(x)/v(x)
-0.100006103515625 <= x <= 0.231639862060546875
peak relative error 1.3e-21 */
u0 = -8.886217500092090678492242071879342025627E1L,
u1 = 6.840109978129177639438792958320783599310E2L,
u2 = 2.042626104514127267855588786511809932433E3L,
u3 = 1.911723903442667422201651063009856064275E3L,
u4 = 7.447065275665887457628865263491667767695E2L,
u5 = 1.132256494121790736268471016493103952637E2L,
u6 = 4.484398885516614191003094714505960972894E0L,
v0 = 1.150830924194461522996462401210374632929E3L,
v1 = 3.399692260848747447377972081399737098610E3L,
v2 = 3.786631705644460255229513563657226008015E3L,
v3 = 1.966450123004478374557778781564114347876E3L,
v4 = 4.741359068914069299837355438370682773122E2L,
v5 = 4.508989649747184050907206782117647852364E1L,
/* v6 = 1.000000000000000000000000000000000000000E0 */
/* lgam (x+2) = .5 x + x s(x)/r(x)
0 <= x <= 1
peak relative error 7.2e-22 */
s0 = 1.454726263410661942989109455292824853344E6L,
s1 = -3.901428390086348447890408306153378922752E6L,
s2 = -6.573568698209374121847873064292963089438E6L,
s3 = -3.319055881485044417245964508099095984643E6L,
s4 = -7.094891568758439227560184618114707107977E5L,
s5 = -6.263426646464505837422314539808112478303E4L,
s6 = -1.684926520999477529949915657519454051529E3L,
r0 = -1.883978160734303518163008696712983134698E7L,
r1 = -2.815206082812062064902202753264922306830E7L,
r2 = -1.600245495251915899081846093343626358398E7L,
r3 = -4.310526301881305003489257052083370058799E6L,
r4 = -5.563807682263923279438235987186184968542E5L,
r5 = -3.027734654434169996032905158145259713083E4L,
r6 = -4.501995652861105629217250715790764371267E2L,
/* r6 = 1.000000000000000000000000000000000000000E0 */
/* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
x >= 8
Peak relative error 1.51e-21
w0 = LS2PI - 0.5 */
w0 = 4.189385332046727417803e-1L,
w1 = 8.333333333333331447505E-2L,
w2 = -2.777777777750349603440E-3L,
w3 = 7.936507795855070755671E-4L,
w4 = -5.952345851765688514613E-4L,
w5 = 8.412723297322498080632E-4L,
w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L;
/* sin(pi*x) assuming x > 2^-1000, if sin(pi*x)==0 the sign is arbitrary */
static long double sin_pi(long double x)
{
int n;
/* spurious inexact if odd int */
x *= 0.5;
x = 2.0*(x - floorl(x)); /* x mod 2.0 */
n = (int)(x*4.0);
n = (n+1)/2;
x -= n*0.5f;
x *= pi;
switch (n) {
default: /* case 4: */
case 0: return __sinl(x, 0.0, 0);
case 1: return __cosl(x, 0.0);
case 2: return __sinl(-x, 0.0, 0);
case 3: return -__cosl(x, 0.0);
}
}
long double lgammal_r(long double x, int *sg) {
long double t, y, z, nadj, p, p1, p2, q, r, w;
union ldshape u = {x};
uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
int sign = u.i.se >> 15;
int i;
*sg = 1;
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
if (ix >= 0x7fff0000)
return x * x;
if (ix < 0x3fc08000) { /* |x|<2**-63, return -log(|x|) */
if (sign) {
*sg = -1;
x = -x;
}
return -logl(x);
}
if (sign) {
x = -x;
t = sin_pi(x);
if (t == 0.0)
return 1.0 / (x-x); /* -integer */
if (t > 0.0)
*sg = -1;
else
t = -t;
nadj = logl(pi / (t * x));
}
/* purge off 1 and 2 (so the sign is ok with downward rounding) */
if ((ix == 0x3fff8000 || ix == 0x40008000) && u.i.m == 0) {
r = 0;
} else if (ix < 0x40008000) { /* x < 2.0 */
if (ix <= 0x3ffee666) { /* 8.99993896484375e-1 */
/* lgamma(x) = lgamma(x+1) - log(x) */
r = -logl(x);
if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */
y = x - 1.0;
i = 0;
} else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */
y = x - (tc - 1.0);
i = 1;
} else { /* x < 0.23 */
y = x;
i = 2;
}
} else {
r = 0.0;
if (ix >= 0x3fffdda6) { /* 1.73162841796875 */
/* [1.7316,2] */
y = x - 2.0;
i = 0;
} else if (ix >= 0x3fff9da6) { /* 1.23162841796875 */
/* [1.23,1.73] */
y = x - tc;
i = 1;
} else {
/* [0.9, 1.23] */
y = x - 1.0;
i = 2;
}
}
switch (i) {
case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
r += 0.5 * y + y * p1/p2;
break;
case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
p = tt + y * p1/p2;
r += (tf + p);
break;
case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
r += (-0.5 * y + p1 / p2);
}
} else if (ix < 0x40028000) { /* 8.0 */
/* x < 8.0 */
i = (int)x;
y = x - (double)i;
p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
r = 0.5 * y + p / q;
z = 1.0;
/* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7:
z *= (y + 6.0); /* FALLTHRU */
case 6:
z *= (y + 5.0); /* FALLTHRU */
case 5:
z *= (y + 4.0); /* FALLTHRU */
case 4:
z *= (y + 3.0); /* FALLTHRU */
case 3:
z *= (y + 2.0); /* FALLTHRU */
r += logl(z);
break;
}
} else if (ix < 0x40418000) { /* 2^66 */
/* 8.0 <= x < 2**66 */
t = logl(x);
z = 1.0 / x;
y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
r = (x - 0.5) * (t - 1.0) + w;
} else /* 2**66 <= x <= inf */
r = x * (logl(x) - 1.0);
if (sign)
r = nadj - r;
return r;
}
long double lgammal(long double x)
{
return lgammal_r(x, &signgam);
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
// TODO: broken implementation to make things compile
long double lgammal_r(long double x, int *sg)
{
return lgamma_r(x, sg);
}
long double lgammal(long double x)
{
return lgamma_r(x, &signgam);
}
#endif /* long double is long */

4
libc/tinymath/tgamma.c
View file

@ -208,3 +208,7 @@ double tgamma(double x)
y = r * z * z;
return y;
}
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
__weak_reference(tgamma, tgammal);
#endif

314
libc/tinymath/tgammal.c Normal file
View file

@ -0,0 +1,314 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│
╚──────────────────────────────────────────────────────────────────────────────╝
│ │
│ Musl Libc │
│ Copyright © 2005-2014 Rich Felker, et al. │
│ │
│ Permission is hereby granted, free of charge, to any person obtaining │
│ a copy of this software and associated documentation files (the │
│ "Software"), to deal in the Software without restriction, including │
│ without limitation the rights to use, copy, modify, merge, publish, │
│ distribute, sublicense, and/or sell copies of the Software, and to │
│ permit persons to whom the Software is furnished to do so, subject to │
│ the following conditions: │
│ │
│ The above copyright notice and this permission notice shall be │
│ included in all copies or substantial portions of the Software. │
│ │
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
│ │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
#include "libc/tinymath/internal.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
asm(".ident\t\"\\n\\n\
OpenBSD libm (ISC License)\\n\
Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Gamma function
*
*
* SYNOPSIS:
*
* long double x, y, tgammal();
*
* y = tgammal( x );
*
*
* DESCRIPTION:
*
* Returns gamma function of the argument. The result is
* correctly signed.
*
* Arguments |x| <= 13 are reduced by recurrence and the function
* approximated by a rational function of degree 7/8 in the
* interval (2,3). Large arguments are handled by Stirling's
* formula. Large negative arguments are made positive using
* a reflection formula.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -40,+40 10000 3.6e-19 7.9e-20
* IEEE -1755,+1755 10000 4.8e-18 6.5e-19
*
* Accuracy for large arguments is dominated by error in powl().
*
*/
/*
tgamma(x+2) = tgamma(x+2) P(x)/Q(x)
0 <= x <= 1
Relative error
n=7, d=8
Peak error = 1.83e-20
Relative error spread = 8.4e-23
*/
static const long double P[8] = {
4.212760487471622013093E-5L,
4.542931960608009155600E-4L,
4.092666828394035500949E-3L,
2.385363243461108252554E-2L,
1.113062816019361559013E-1L,
3.629515436640239168939E-1L,
8.378004301573126728826E-1L,
1.000000000000000000009E0L,
};
static const long double Q[9] = {
-1.397148517476170440917E-5L,
2.346584059160635244282E-4L,
-1.237799246653152231188E-3L,
-7.955933682494738320586E-4L,
2.773706565840072979165E-2L,
-4.633887671244534213831E-2L,
-2.243510905670329164562E-1L,
4.150160950588455434583E-1L,
9.999999999999999999908E-1L,
};
/*
static const long double P[] = {
-3.01525602666895735709e0L,
-3.25157411956062339893e1L,
-2.92929976820724030353e2L,
-1.70730828800510297666e3L,
-7.96667499622741999770e3L,
-2.59780216007146401957e4L,
-5.99650230220855581642e4L,
-7.15743521530849602425e4L
};
static const long double Q[] = {
1.00000000000000000000e0L,
-1.67955233807178858919e1L,
8.85946791747759881659e1L,
5.69440799097468430177e1L,
-1.98526250512761318471e3L,
3.31667508019495079814e3L,
1.60577839621734713377e4L,
-2.97045081369399940529e4L,
-7.15743521530849602412e4L
};
*/
#define MAXGAML 1755.455L
/*static const long double LOGPI = 1.14472988584940017414L;*/
/* Stirling's formula for the gamma function
tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
z(x) = x
13 <= x <= 1024
Relative error
n=8, d=0
Peak error = 9.44e-21
Relative error spread = 8.8e-4
*/
static const long double STIR[9] = {
7.147391378143610789273E-4L,
-2.363848809501759061727E-5L,
-5.950237554056330156018E-4L,
6.989332260623193171870E-5L,
7.840334842744753003862E-4L,
-2.294719747873185405699E-4L,
-2.681327161876304418288E-3L,
3.472222222230075327854E-3L,
8.333333333333331800504E-2L,
};
#define MAXSTIR 1024.0L
static const long double SQTPI = 2.50662827463100050242E0L;
/* 1/tgamma(x) = z P(z)
* z(x) = 1/x
* 0 < x < 0.03125
* Peak relative error 4.2e-23
*/
static const long double S[9] = {
-1.193945051381510095614E-3L,
7.220599478036909672331E-3L,
-9.622023360406271645744E-3L,
-4.219773360705915470089E-2L,
1.665386113720805206758E-1L,
-4.200263503403344054473E-2L,
-6.558780715202540684668E-1L,
5.772156649015328608253E-1L,
1.000000000000000000000E0L,
};
/* 1/tgamma(-x) = z P(z)
* z(x) = 1/x
* 0 < x < 0.03125
* Peak relative error 5.16e-23
* Relative error spread = 2.5e-24
*/
static const long double SN[9] = {
1.133374167243894382010E-3L,
7.220837261893170325704E-3L,
9.621911155035976733706E-3L,
-4.219773343731191721664E-2L,
-1.665386113944413519335E-1L,
-4.200263503402112910504E-2L,
6.558780715202536547116E-1L,
5.772156649015328608727E-1L,
-1.000000000000000000000E0L,
};
static const long double PIL = 3.1415926535897932384626L;
/* Gamma function computed by Stirling's formula.
*/
static long double stirf(long double x)
{
long double y, w, v;
w = 1.0/x;
/* For large x, use rational coefficients from the analytical expansion. */
if (x > 1024.0)
w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w
+ 1.0;
else
w = 1.0 + w * __polevll(w, STIR, 8);
y = expl(x);
if (x > MAXSTIR) { /* Avoid overflow in pow() */
v = powl(x, 0.5L * x - 0.25L);
y = v * (v / y);
} else {
y = powl(x, x - 0.5L) / y;
}
y = SQTPI * y * w;
return y;
}
long double tgammal(long double x)
{
long double p, q, z;
if (!isfinite(x))
return x + INFINITY;
q = fabsl(x);
if (q > 13.0) {
if (x < 0.0) {
p = floorl(q);
z = q - p;
if (z == 0)
return 0 / z;
if (q > MAXGAML) {
z = 0;
} else {
if (z > 0.5) {
p += 1.0;
z = q - p;
}
z = q * sinl(PIL * z);
z = fabsl(z) * stirf(q);
z = PIL/z;
}
if (0.5 * p == floorl(q * 0.5))
z = -z;
} else if (x > MAXGAML) {
z = x * 0x1p16383L;
} else {
z = stirf(x);
}
return z;
}
z = 1.0;
while (x >= 3.0) {
x -= 1.0;
z *= x;
}
while (x < -0.03125L) {
z /= x;
x += 1.0;
}
if (x <= 0.03125L)
goto small;
while (x < 2.0) {
z /= x;
x += 1.0;
}
if (x == 2.0)
return z;
x -= 2.0;
p = __polevll(x, P, 7);
q = __polevll(x, Q, 8);
z = z * p / q;
return z;
small:
/* z==1 if x was originally +-0 */
if (x == 0 && z != 1)
return x / x;
if (x < 0.0) {
x = -x;
q = z / (x * __polevll(x, SN, 8));
} else
q = z / (x * __polevll(x, S, 8));
return q;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
// TODO: broken implementation to make things compile
long double tgammal(long double x)
{
return tgamma(x);
}
#endif /* long double is long */