1 /*
2 * Copyright (c) 1996-2001 Lucent Technologies.
3 * See README file for details.
4 *
5 *
6 miscellaneous functions that may not be defined in the math
7 libraries. The implementations are crude.
8 lflgamma(x) -- log(gamma(x))
9 lferf(x) -- erf(x)
10 lferfc(x) -- erfc(x)
11 lfdaws(x) -- dawson's function
12
13 where required, these must be #define'd in local.h.
14
15 also includes
16 ptail(x) -- exp(x*x/2)*int_{-\infty}^x exp(-u^2/2)du for x < -6.
17 logit(x) -- logistic function.
18 expit(x) -- inverse of logit.
19 */
20
21 #include "local.h"
22
23 double lferfc();
24
lferf(x)25 double lferf(x)
26 double x;
27 { static double val[] = { 0.0, 0.52049987781304674,
28 0.84270079294971501, 0.96610514647531076, 0.99532226501895282,
29 0.99959304798255499, 0.99997790950300125 };
30 double h, xx, y, z, f0, f1, f2;
31 int m, j;
32 if (x<0) return(-lferf(-x));
33 if (x>3.2) return(1-lferfc(x));
34 m = (int) (2*x+0.5);
35 xx= ((double)m)/2;
36 h = x-xx; y = h;
37 f0 = val[m];
38 f1 = 2*exp(-xx*xx)/SQRPI;
39 z = f0+h*f1;
40 j = 0;
41 while (fabs(y)>1.0e-12)
42 { f2 = -2*j*f0-2*xx*f1;
43 f0 = f1; f1 = f2;
44 y *= h/(j+2);
45 z += y*f2;
46 j++;
47 }
48 return(z);
49 }
50
lferfc(x)51 double lferfc(x)
52 double x;
53 { if (x<0) return(1+lferf(-x));
54 if (x<2.5) return(1-lferf(x));
55 return(exp(-x*x)/(x*SQRPI));
56 }
57
lflgamma(x)58 double lflgamma(x)
59 double x;
60 { double x1;
61 static double ilg[] = { 0.0, 0.0, 0.69314718055994529,
62 1.791759469228055, 3.1780538303479458, 4.7874917427820458, 6.5792512120101012,
63 8.5251613610654147, 10.604602902745251, 12.801827480081469 };
64 static double hlg[] = { 0.57236494292470008, -0.12078223763524520,
65 0.28468287047291918, 1.20097360234707430, 2.45373657084244230,
66 3.95781396761871650, 5.66256205985714270, 7.53436423675873360,
67 9.54926725730099870, 11.68933342079726900 };
68
69 if (x<=0.0) return(0.0);
70 if (x<10)
71 { if (x==(int)x) return(ilg[(int)x-1]);
72 if ((x-0.5)==(int)(x-0.5)) return(hlg[(int)(x-0.5)]);
73 }
74 if (x<3) return(lflgamma(x+1)-log(x));
75
76 x1 = x-1;
77 return(HL2PI+(x1+0.5)*log(x1)-x1+1/(12*x1));
78 }
79
lfdaws(x)80 double lfdaws(x)
81 double x;
82 { static double val[] = {
83 0, 0.24485619356002, 0.46034428261948, 0.62399959848185, 0.72477845900708,
84 0.76388186132749, 0.75213621001998, 0.70541701910853, 0.63998807456541,
85 0.56917098836654, 0.50187821196415, 0.44274283060424, 0.39316687916687,
86 0.35260646480842, 0.31964847250685, 0.29271122077502, 0.27039629581340,
87 0.25160207761769, 0.23551176224443, 0.22153505358518, 0.20924575719548,
88 0.19833146819662, 0.18855782729305, 0.17974461154688, 0.17175005072385 };
89 double h, f0, f1, f2, y, z, xx;
90 int j, m;
91 if (x<0) return(-daws(-x));
92 if (x>6)
93 { /* Tail series: 1/x + 1/x^3 + 1.3/x^5 + 1.3.5/x^7 + ... */
94 y = z = 1/x;
95 j = 0;
96 while (((f0=(2*j+1)/(x*x))<1) && (y>1.0e-10*z))
97 { y *= f0;
98 z += y;
99 j++;
100 }
101 return(z);
102 }
103 m = (int) (4*x);
104 h = x-0.25*m;
105 if (h>0.125)
106 { m++;
107 h = h-0.25;
108 }
109 xx = 0.25*m;
110 f0 = val[m];
111 f1 = 1-xx*f0;
112 z = f0+h*f1;
113 y = h;
114 j = 2;
115 while (fabs(y)>z*1.0e-10)
116 { f2 = -(j-1)*f0-xx*f1;
117 y *= h/j;
118 z += y*f2;
119 f0 = f1; f1 = f2;
120 j++;
121 }
122 return(z);
123 }
124
ptail(x)125 double ptail(x) /* exp(x*x/2)*int_{-\infty}^x exp(-u^2/2)du for x < -6 */
126 double x;
127 { double y, z, f0;
128 int j;
129 y = z = -1.0/x;
130 j = 0;
131 while ((fabs(f0= -(2*j+1)/(x*x))<1) && (fabs(y)>1.0e-10*z))
132 { y *= f0;
133 z += y;
134 j++;
135 }
136 return(z);
137 }
138
logit(x)139 double logit(x)
140 double x;
141 { return(log(x/(1-x)));
142 }
143
expit(x)144 double expit(x)
145 double x;
146 { double u;
147 if (x<0)
148 { u = exp(x);
149 return(u/(1+u));
150 }
151 return(1/(1+exp(-x)));
152 }
153
factorial(n)154 int factorial(n)
155 int n;
156 { if (n<=1) return(1.0);
157 return(n*factorial(n-1));
158 }
159