htestu-0.1.0.0: cbits/testu/src/num2.c
/*************************************************************************\
*
* Package: MyLib
* File: num2.c
* Environment: ANSI C
*
* Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
* e-mail: lecuyer@iro.umontreal.ca
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted without a fee for private, research,
* academic, or other non-commercial purposes.
* Any use of this software in a commercial environment requires a
* written licence from the copyright owner.
*
* Any changes made to this package must be clearly identified as such.
*
* In scientific publications which used this software, a reference to it
* would be appreciated.
*
* Redistributions of source code must retain this copyright notice
* and the following disclaimer.
*
* THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
\*************************************************************************/
#include "num2.h"
#include "util.h"
#include "num.h"
#include <math.h>
#include <float.h>
#include <stdlib.h>
#define EPSILON 1.0E-15
#define MAXI 50
/* The factorials n! from n = 0 to n = 170 */
static double Factorials[] = {
1,
1,
2,
6,
24,
120,
720,
5040,
40320,
362880,
3628800,
39916800,
479001600,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
1.21645100408832e+17,
2.43290200817664e+18,
5.109094217170944e+19,
1.124000727777608e+21,
2.585201673888498e+22,
6.204484017332394e+23,
1.551121004333099e+25,
4.032914611266057e+26,
1.088886945041835e+28,
3.048883446117138e+29,
8.841761993739701e+30,
2.65252859812191e+32,
8.222838654177922e+33,
2.631308369336935e+35,
8.683317618811886e+36,
2.952327990396041e+38,
1.033314796638614e+40,
3.719933267899012e+41,
1.376375309122634e+43,
5.23022617466601e+44,
2.039788208119744e+46,
8.159152832478977e+47,
3.34525266131638e+49,
1.40500611775288e+51,
6.041526306337383e+52,
2.658271574788449e+54,
1.196222208654802e+56,
5.502622159812088e+57,
2.586232415111682e+59,
1.241391559253607e+61,
6.082818640342675e+62,
3.041409320171338e+64,
1.551118753287382e+66,
8.065817517094388e+67,
4.274883284060025e+69,
2.308436973392414e+71,
1.269640335365828e+73,
7.109985878048635e+74,
4.052691950487722e+76,
2.350561331282879e+78,
1.386831185456899e+80,
8.320987112741392e+81,
5.075802138772248e+83,
3.146997326038794e+85,
1.98260831540444e+87,
1.268869321858842e+89,
8.247650592082472e+90,
5.443449390774431e+92,
3.647111091818868e+94,
2.480035542436831e+96,
1.711224524281413e+98,
1.197857166996989e+100,
8.504785885678622e+101,
6.123445837688608e+103,
4.470115461512683e+105,
3.307885441519386e+107,
2.480914081139539e+109,
1.88549470166605e+111,
1.451830920282858e+113,
1.132428117820629e+115,
8.946182130782973e+116,
7.156945704626378e+118,
5.797126020747366e+120,
4.75364333701284e+122,
3.945523969720657e+124,
3.314240134565352e+126,
2.817104114380549e+128,
2.422709538367272e+130,
2.107757298379527e+132,
1.854826422573984e+134,
1.650795516090845e+136,
1.485715964481761e+138,
1.352001527678402e+140,
1.24384140546413e+142,
1.156772507081641e+144,
1.087366156656742e+146,
1.032997848823905e+148,
9.916779348709491e+149,
9.619275968248206e+151,
9.426890448883242e+153,
9.33262154439441e+155,
9.33262154439441e+157,
9.425947759838354e+159,
9.614466715035121e+161,
9.902900716486175e+163,
1.029901674514562e+166,
1.08139675824029e+168,
1.146280563734708e+170,
1.226520203196137e+172,
1.324641819451828e+174,
1.443859583202493e+176,
1.588245541522742e+178,
1.762952551090244e+180,
1.974506857221073e+182,
2.231192748659812e+184,
2.543559733472186e+186,
2.925093693493014e+188,
3.393108684451897e+190,
3.969937160808719e+192,
4.684525849754288e+194,
5.574585761207603e+196,
6.689502913449124e+198,
8.09429852527344e+200,
9.875044200833598e+202,
1.214630436702532e+205,
1.50614174151114e+207,
1.882677176888925e+209,
2.372173242880046e+211,
3.012660018457658e+213,
3.856204823625803e+215,
4.974504222477285e+217,
6.466855489220472e+219,
8.471580690878817e+221,
1.118248651196004e+224,
1.487270706090685e+226,
1.992942746161518e+228,
2.69047270731805e+230,
3.659042881952547e+232,
5.01288874827499e+234,
6.917786472619486e+236,
9.615723196941086e+238,
1.346201247571752e+241,
1.89814375907617e+243,
2.695364137888161e+245,
3.854370717180071e+247,
5.550293832739301e+249,
8.047926057471987e+251,
1.17499720439091e+254,
1.727245890454638e+256,
2.556323917872864e+258,
3.808922637630567e+260,
5.71338395644585e+262,
8.627209774233235e+264,
1.311335885683452e+267,
2.006343905095681e+269,
3.089769613847349e+271,
4.789142901463391e+273,
7.471062926282891e+275,
1.172956879426414e+278,
1.853271869493734e+280,
2.946702272495037e+282,
4.714723635992059e+284,
7.590705053947215e+286,
1.229694218739449e+289,
2.004401576545302e+291,
3.287218585534294e+293,
5.423910666131586e+295,
9.003691705778433e+297,
1.503616514864998e+300,
2.526075744973197e+302,
4.269068009004703e+304,
7.257415615307994e+306
};
#define MLIM 50
/* The natural logarithm of factorials n! from n = 0 to 50 */
static double LnFactorials[MLIM + 1] = {
0.,
0.,
0.6931471805599453,
1.791759469228055,
3.178053830347946,
4.787491742782046,
6.579251212010101,
8.525161361065415,
10.60460290274525,
12.80182748008147,
15.10441257307552,
17.50230784587389,
19.98721449566188,
22.55216385312342,
25.19122118273868,
27.89927138384088,
30.67186010608066,
33.50507345013688,
36.39544520803305,
39.33988418719949,
42.33561646075348,
45.3801388984769,
48.47118135183522,
51.60667556776437,
54.7847293981123,
58.00360522298051,
61.26170176100199,
64.55753862700632,
67.88974313718154,
71.257038967168,
74.65823634883016,
78.09222355331529,
81.55795945611503,
85.05446701758153,
88.58082754219767,
92.13617560368708,
95.7196945421432,
99.3306124547874,
102.9681986145138,
106.6317602606434,
110.3206397147574,
114.0342117814617,
117.7718813997451,
121.5330815154386,
125.3172711493569,
129.1239336391272,
132.9525750356163,
136.8027226373264,
140.6739236482343,
144.5657439463449,
148.477766951773
};
/*=========================================================================*/
double num2_Factorial (int n)
{
util_Assert (n >= 0, "num2_Factorial: n < 0");
if (n <= 170)
return Factorials[n];
util_Warning (1, "num2_Factorial: n > 170: return inf");
return 1.0 / 0.0;
}
/*=========================================================================*/
double num2_LnFactorial (int n)
{
util_Assert (n >= 0, "num2_LnFactorial: n < 0");
if (n <= MLIM) {
return LnFactorials[n];
} else {
double x = (double) (n + 1);
double y = 1.0 / (x * x);
double z = ((-(5.95238095238E-4 * y) + 7.936500793651E-4) * y -
2.7777777777778E-3) * y + 8.3333333333333E-2;
z = ((x - 0.5) * log (x) - x) + 9.1893853320467E-1 + z / x;
return z;
}
}
/*=========================================================================*/
#ifndef HAVE_LGAMMA
/* The new standard ISO_C99 includes the lgamma function in math.h;
otherwise, we shall have to use our own. */
double num2_LnGamma (double x)
{
const double xlimbig = 1.0 / DBL_EPSILON;
const double xlim1 = 18.0;
const double dk2 = 0.91893853320467274178; /* Ln (sqrt (2 Pi)) */
const double dk1 = 0.9574186990510627;
const int N = 15; /* Degree of Chebyshev polynomial */
double y = 0, z = 0;
int i, k;
/* Chebyshev coefficients for lnGamma (x + 3), 0 <= x <= 1 In Yudell Luke:
The special functions and their approximations, Vol. II, Academic Press,
p. 301, 1969. There is an error in the additive constant in the formula:
(Ln (2)). */
static const double A[] = {
0.52854303698223459887,
0.54987644612141411418,
0.02073980061613665136,
-0.00056916770421543842,
0.00002324587210400169,
-0.00000113060758570393,
0.00000006065653098948,
-0.00000000346284357770,
0.00000000020624998806,
-0.00000000001266351116,
0.00000000000079531007,
-0.00000000000005082077,
0.00000000000000329187,
-0.00000000000000021556,
0.00000000000000001424,
-0.00000000000000000095
};
util_Assert (x > 0.0, "num2_LnGamma: accepts only x > 0");
if (x > xlim1) {
if (x > xlimbig)
y = 0.0;
else
y = 1.0 / (x * x);
z = ((-(5.95238095238E-4 * y) + 7.936500793651E-4) * y -
2.7777777777778E-3) * y + 8.3333333333333E-2;
z = ((x - 0.5) * log (x) - x) + dk2 + z / x;
return z;
} else if (x > 4.0) {
k = (int) x;
z = x - k;
y = 1.0;
for (i = 3; i < k; i++)
y *= z + i;
y = log (y);
} else if (x <= 0.0) {
return DBL_MAX;
} else if (x < 3.0) {
k = (int) x;
z = x - k;
y = 1.0;
for (i = 2; i >= k; i--)
y *= z + i;
y = -log (y);
} else { /* 3 <= x <= 4 */
z = x - 3.0;
y = 0.0;
}
z = num2_EvalCheby (A, N, 2.0 * z - 1.0);
return z + dk1 + y;
}
#endif
/*=========================================================================*/
#define NLIM 100 /* pour eviter les debordements */
double num2_Combination (int n, int s)
{
double Res;
int i;
int Diff;
if (s == 0 || s == n)
return 1.0;
if (s < 0) {
util_Warning (1, "num2_Combination: s < 0");
return 0.0;
}
if (s > n) {
util_Warning (1, "num2_Combination: s > n");
return 0.0;
}
if (s > (n / 2))
s = n - s;
if (n <= NLIM) {
Res = 1.0;
Diff = n - s;
for (i = 1; i <= s; i++) {
Res = (Res * (double) (Diff + i)) / (double) (i);
}
return Res;
} else {
Res = (num2_LnFactorial (n) - num2_LnFactorial (s))
- num2_LnFactorial (n - s);
return exp (Res);
}
}
/*=========================================================================*/
#ifndef HAVE_LOG1P
double num2_log1p (double x)
{
/* returns a value equivalent to log (1 + x) accurate also for small x. */
if (fabs (x) > 0.1) {
return log (1.0 + x);
} else {
double term = x;
double sum = x;
int s = 2;
while (fabs (term) > EPSILON * fabs (sum) && s < MAXI) {
term *= -x;
sum += term / s;
s++;
}
return sum;
}
}
#endif
/*=========================================================================*/
void num2_CalcMatStirling (double ***M, int m, int n)
/* Calcul des elements de la matrice MatStirling [0..m, 0..n]. */
{
int i, j, k;
*M = (double **) util_Calloc ((size_t) (m + 1), sizeof (double *));
for (i = 0; i <= m; i++)
(*M)[i] = (double *) util_Calloc ((size_t) (n + 1), sizeof (double));
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++) {
(*M)[i][j] = 0.0;
}
}
(*M)[0][0] = 1.0;
for (j = 1; j <= n; j++) {
(*M)[0][j] = 0.0;
if (j <= m) {
k = j - 1;
(*M)[j][j] = 1.0;
} else
k = m;
for (i = 1; i <= k; i++) {
(*M)[i][j] = (double) (i) * (*M)[i][j - 1] + (*M)[i - 1][j - 1];
}
}
}
/*=========================================================================*/
void num2_FreeMatStirling (double ***M, int m)
{
int i;
for (i = 0; i <= m; i++)
free ((*M)[i]);
free (*M);
*M = NULL;
}
/*=========================================================================*/
double num2_VolumeSphere (double pLR, int k)
/* Returns volume of unit sphere in dimension k, norm p */
{
const double eps = 2.0 * DBL_EPSILON;
int p = pLR;
double kLR = (double) k;
double Vol;
int s;
util_Assert (pLR >= 0.0, "num2_VolumeSphere: p < 0");
if (fabs (pLR - p) <= eps) {
switch (p) {
case 0:
return num_TwoExp[k];
break;
case 1:
return num_TwoExp[k] / num2_Factorial (k);
break;
case 2:
if ((k % 2) == 0) {
return pow (num_Pi, kLR / 2.0) / num2_Factorial (k / 2);
} else {
s = (k + 1) / 2;
return pow (num_Pi, (double) (s) - 1.0) * num2_Factorial (s) *
num_TwoExp[2 * s] / num2_Factorial (2 * s);
}
break;
default:
break;
}
}
Vol = kLR * (num_Ln2 + num2_LnGamma (1.0 + 1.0 / pLR)) -
num2_LnGamma (1.0 + kLR / pLR);
return exp (Vol);
}
/*=========================================================================*/
double num2_EvalCheby (const double A[], int N, double x)
{
int j;
double xx;
double b0, b1, b2;
util_Warning (fabs (x) > 1.0,
"Chebychev polynomial evaluated at x outside [-1, 1]");
xx = 2.0 * x;
b0 = 0.0;
b1 = 0.0;
for (j = N; j >= 0; j--) {
b2 = b1;
b1 = b0;
b0 = (xx * b1 - b2) + A[j];
}
return (b0 - b2) / 2.0;
}
/*=========================================================================*/
#define DEGREE 6
double num2_BesselK025 (double x)
{
double rac;
double xx;
double temp;
double Res;
double C;
double B;
int j;
static const double c[8] = {
32177591145.0,
2099336339520.0,
16281990144000.0,
34611957596160.0,
26640289628160.0,
7901666082816.0,
755914244096.0
};
static const double b[8] = {
75293843625.0,
2891283595200.0,
18691126272000.0,
36807140966400.0,
27348959232000.0,
7972533043200.0,
755914244096.0
};
if (x < 1.E-300)
return DBL_MAX;
/*------------------------------------------------------------------*/
/* x > 0.6 => approximation asymptotique rationnelle dans Luke: */
/* Yudell L.Luke "Mathematical functions and their approximations", */
/* Academic Press Inc. New York, 1975, p.371 */
/*------------------------------------------------------------------*/
if (x >= 0.6) {
B = b[DEGREE];
C = c[DEGREE];
for (j = DEGREE; j >= 1; j--) {
B = B * x + b[j - 1];
C = C * x + c[j - 1];
}
Res = sqrt (num_Pi / (2.0 * x)) * exp (-x) * (C / B);
return Res;
}
/*------------------------------------------------------------------*/
/* x < 0.6 => la serie de K_{1/4} = Pi/Sqrt(2) [I_{-1/4} - I_{1/4}] */
/*------------------------------------------------------------------*/
xx = x * x;
rac = pow (x / 2.0, 0.25);
Res = (((xx / 1386.0 + 1.0 / 42.0) * xx + 1.0 / 3.0) * xx + 1.0) /
(1.225416702465177 * rac);
temp = (((xx / 3510.0 + 1.0 / 90.0) * xx + 0.2) * xx + 1.0) * rac /
0.906402477055477;
Res = num_Pi * (Res - temp) / num_Rac2;
return Res;
}
#undef DEGREE
/*=========================================================================*/
double num2_Digamma (double x)
{
static const double C7[] = {
1.3524999667726346383e4, 4.5285601699547289655e4,
4.5135168469736662555e4, 1.8529011818582610168e4,
3.3291525149406935532e3, 2.4068032474357201831e2,
5.1577892000139084710, 6.2283506918984745826e-3
};
static const double D7[] = {
6.9389111753763444376e-7, 1.9768574263046736421e4,
4.1255160835353832333e4, 2.9390287119932681918e4,
9.0819666074855170271e3, 1.2447477785670856039e3,
6.7429129516378593773e1, 1.0
};
static const double C4[] = {
-2.728175751315296783e-15, -6.481571237661965099e-1,
-4.486165439180193579, -7.016772277667586642, -2.129404451310105168
};
static const double D4[] = {
7.777885485229616042, 5.461177381032150702e1,
8.929207004818613702e1, 3.227034937911433614e1, 1.0
};
double prodPj = 0.0;
double prodQj = 0.0;
double digX = 0.0;
if (x >= 3.0) {
double x2 = 1.0 / (x * x);
int j;
for (j = 4; j >= 0; j--) {
prodPj = prodPj * x2 + C4[j];
prodQj = prodQj * x2 + D4[j];
}
digX = log (x) - (0.5 / x) + (prodPj / prodQj);
} else if (x >= 0.5) {
const double X0 = 1.46163214496836234126;
int j;
for (j = 7; j >= 0; j--) {
prodPj = x * prodPj + C7[j];
prodQj = x * prodQj + D7[j];
}
digX = (x - X0) * (prodPj / prodQj);
} else {
double f = (1.0 - x) - floor (1.0 - x);
digX = num2_Digamma (1.0 - x) + num_Pi / tan (num_Pi * f);
}
return digX;
}