diff options
author | Gwenolé Beauchesne <gbeauchesne@mandriva.org> | 2003-06-04 18:44:09 +0000 |
---|---|---|
committer | Gwenolé Beauchesne <gbeauchesne@mandriva.org> | 2003-06-04 18:44:09 +0000 |
commit | 4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b (patch) | |
tree | acd4001a266a8713495af7f1b2102b61e67113b0 /mdk-stage1/dietlibc/libm/gamma.c | |
parent | 71b111ec6c4671667a19c6fbe0023d33422535d7 (diff) | |
download | drakx-4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b.tar drakx-4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b.tar.gz drakx-4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b.tar.bz2 drakx-4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b.tar.xz drakx-4cd6a4a5d7e49d54d53dcf4a6f3393d50bd88e8b.zip |
Import dietlibc 0.22 + other fixes for AMD64
Diffstat (limited to 'mdk-stage1/dietlibc/libm/gamma.c')
-rw-r--r-- | mdk-stage1/dietlibc/libm/gamma.c | 98 |
1 files changed, 98 insertions, 0 deletions
diff --git a/mdk-stage1/dietlibc/libm/gamma.c b/mdk-stage1/dietlibc/libm/gamma.c new file mode 100644 index 000000000..d5f3e4275 --- /dev/null +++ b/mdk-stage1/dietlibc/libm/gamma.c @@ -0,0 +1,98 @@ +#include "dietlibm.h" + +/*--------------------------------------------------------------------------* + +Name gamma, lgamma - gamma function + +Usage double gamma (double x); + double lgamma(double x); + extern int signgam; + +Prototype in math.h + +Description gamma returns the logarithm of the absolute value of the + gamma function. So it is possible â(x) for very large x. + The sign is stored in signgam, a extern variable + overwritten during every call to gamma(). lgamma() is + a synonym for gamma(). + You can calculate â(x) by the following sequence: + + double gammafunction(double x) + { double y=exp(gamma(x)); + + return signgam ? -y : +y; + } + +Return value gamma returns a value in range (-0.1208, +oo). For a input + value of zero, it returns +oo and errno is set to: + + ERANGE Result out of range + +*---------------------------------------------------------------------------*/ + +#include <stdlib.h> +#include <math.h> + +#define B0 + 1.0l/ 6/ 1/ 2 +#define B1 - 1.0l/ 30/ 3/ 4 +#define B2 + 1.0l/ 42/ 5/ 6 +#define B3 - 1.0l/ 30/ 7/ 8 +#define B4 + 5.0l/ 66/ 9/10 +#define B5 - 691.0l/2730/11/12 +#define B6 + 7.0l/ 6/13/14 +#define B7 - 3617.0l/ 510/15/16 +#define B8 + 43867.0l/ 798/17/18 +#define B9 - 174611.0l/ 330/19/20 +#define B10 + 854513.0l/ 138/21/22 +#define B11 - 236364091.0l/2730/23/24 +#define B12 + 8553103.0l/ 6/25/26 + +static const double coeff[] = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, B10 }; +int signgam; + +#define EXPL(x) (((short *)&x)[4] & 0x7FFF) + +static double logfact ( long double x ) +{ + long double z = 2. * M_PI * x; + register int e = EXPL (x); + + static unsigned char list [] = { 6, 4, 3, 3, 2, 2 }; + + return (log(x) - 1) * x + 0.5*log(z) + __poly (1./(x*x), e<0x4003 ? 10 : (e>0x4008 ? 1 : list [e-0x4003] ), coeff) / x; +} + + +double lgamma ( double x ) +{ + register int k = floor (x); + long double w; + long double y; + long double z; + + signgam = 0; + + if ( k >= 7 ) + return logfact (x-1); + + if ( k == x ) + switch (k) { + case 1 : + case 2 : return 0.000000000000000000000000000l; + case 3 : return 0.693147180559945309432805516l; + case 4 : return 1.791759469228055000858148560l; + case 5 : return 3.178053830347945619723759592l; + case 6 : return 4.787491742782045994244981560l; + default: return 1./0.; /* ignore the gcc warning, this is intentional */ + } + + z = logfact (y = x - k + 7.0 - 1); + w = 1; + for ( k = 7 - k; k--; ) + w *= y, y -= 1.; + + signgam = k >= 0 ? 0 : k & 1; + return z - log (w); +} + +double gamma ( double val ) __attribute__ ((weak,alias("lgamma"))); |