1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
|
#include <stdio.h>
#include <math.h>
#include "../cephes.h"
#undef fabs
#include "misc.h"
/*
Inverse of the (regularised) incomplete Gamma integral.
Given a, find x such that igam(a, x) = y.
For y not small, we just use igami(a, 1-y) (inverse of the complemented
incomplete Gamma integral). For y small, however, 1-y is about 1, and we
lose digits.
*/
extern double MACHEP;
static double
gammainc(double x, double params[2])
{
return cephes_igam(params[0], x) - params[1];
}
double
gammaincinv(double a, double y)
{
double lo = 0.0, hi;
double flo = -y, fhi = 0.25 - y;
double params[2];
double best_x, best_f;
fsolve_result_t r;
if (a <= 0.0 || y <= 0.0 || y > 0.25) {
return cephes_igami(a, 1-y);
}
params[0] = a;
params[1] = y;
hi = cephes_igami(a, 0.75);
/* I found Newton to be unreliable. Also, after we generate a small
interval by bisection above, false position will do a large step
from an interval of width ~1e-4 to ~1e-14 in one step (a=10, x=0.05,
but similiar for other values).
*/
r = false_position(&lo, &flo, &hi, &fhi,
(objective_function)gammainc, params,
MACHEP, MACHEP, 1e-2*a,
&best_x, &best_f);
if (r == FSOLVE_NOT_BRACKET) {
best_x = 0.0;
}
return best_x;
}
|
