summaryrefslogtreecommitdiff
path: root/modules/randlib/src/c/igngeom.c
blob: 355cbce3475f9213ac12858eb1ab05f4b317535c (plain)
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
/*
 *  PURPOSE
 *     generate a random deviate from G(p) : the geometric
 *     law. If a r.v. X ~ G(p), X is the number of Bernouilli trials
 *     (B(p)) until succes is met. So X take its values in
 *
 *         {1, 2, 3, ...., }
 *
 *     and  P(X=i) = p * (1-p)^(i-1)
 *
 *  METHOD
 *     inversion of the cdf leads to :
 *
 *     (1)   X = 1 + floor( log(1-u) / log(1-p) )
 *
 *     u being a random deviate from U[0,1).
 *
 *     by taking into account that 1-u follows also U(0,1)) this may be
 *     replaced with X = ceil( log(u) / log(1-p) ) or 1 + floor(log(u)/log(1-p))
 *     which needs less work. But as ranf() provides number in [0,1[ , 0 may be
 *     gotten and these formulae may give then +oo.
 *
 *     With ranf() the max number is 1 - 2^(-32). This let us choose a safe min
 *     value for p (to avoid a +oo due to log(1-p)) in the following manner :
 *
 *        the max is gotten for  M = log(2^(-32)) / (-p)
 *
 *     (for very small |x|, the accurate func logp1(x):=log(1+x) return simply x)
 *
 *     and we want  M <= Rmax (near 1.798+308 in ieee 754 double)
 *
 *     so p >= 32 log(2)/Rmax which is near 1.234e-307 ; Says pmin = 1.3e-307.
 *     (anyway the results gotten for such small values of p are certainly
 *      not meaningful...)
 *
 *  NOTE
 *     this function returns a double instead of an int type : this is
 *     to avoid an extra conversion because in scilab it will be a double.
 *
 *  ASSUMPTION
 *     p must be in [pmin,1]  (to do at the calling level).
 *
 *  AUTHOR
 *     Bruno Pincon (<Bruno.Pincon@iecn.u-nancy.fr>)
 *
 */
#include <math.h>
#include "grand.h"

/* the external functions used  here : */
double F2C(logp1)(double *x);  /* celle-ci est ds SCI/modules/elementary_functions/src/fortran/watan.f */
double igngeom(double p)
{
    static double p_save = 1.0, ln_1_m_p = 0.0;
    double u;

    if ( p == 1 )
    {
        return ( 1.0 );
    }
    else if ( p != p_save )   /* => recompute log(1-p) */
    {
        p_save = p;
        u = -p;
        ln_1_m_p = F2C(logp1)(&u);
    };

    u = -C2F(ranf)();
    return ( floor( 1.0 + F2C(logp1)(&u) / ln_1_m_p) );
}