4
Copyright (C) 2005 Gabor Csardi <csardi@rmki.kfki.hu>
5
MTA RMKI, Konkoly-Thege Miklos st. 29-33, Budapest 1121, Hungary
7
This program is free software; you can redistribute it and/or modify
8
it under the terms of the GNU General Public License as published by
9
the Free Software Foundation; either version 2 of the License, or
10
(at your option) any later version.
12
This program is distributed in the hope that it will be useful,
13
but WITHOUT ANY WARRANTY; without even the implied warranty of
14
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15
GNU General Public License for more details.
17
You should have received a copy of the GNU General Public License
18
along with this program; if not, write to the Free Software
19
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
29
#include "igraph_math.h"
31
int igraph_rng_inited = 0;
34
#ifndef USING_R /* R provides a replacement */
35
/* expm1 replacement */
36
static double expm1 (double x)
40
/* Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ... */
44
double term = x / 1.0;
51
while (fabs(term) > fabs(sum) * 2.22e-16);
56
return expl(x) - 1.0L;
62
#ifndef USING_R /* R provides a replacement */
63
/* rint replacement */
64
static double rint (double x)
66
return ( (x<0.) ? -floor(-x+.5) : floor(x+.5) );
72
static float rintf (float x)
74
return ( (x<(float)0.) ? -(float)floor(-x+.5) : (float)floor(x+.5) );
81
* This function appends the rest of the needed random number to the
85
int igraph_random_sample_alga(igraph_vector_t *res, igraph_integer_t l, igraph_integer_t h,
86
igraph_integer_t length) {
87
igraph_real_t N=h-l+1;
88
igraph_real_t n=length;
90
igraph_real_t top=N-n;
91
igraph_real_t Nreal=N;
93
igraph_real_t V, quot;
105
quot=(quot*top)/Nreal;
108
igraph_vector_push_back(res, l); /* allocated */
109
Nreal=-1.0+Nreal; n=-1+n;
112
S=floor(round(Nreal)*RNG_UNIF01());
114
igraph_vector_push_back(res, l); /* allocated */
121
* \function igraph_random_sample
122
* \brief Generates an increasing random sequence of integers.
125
* This function generates an incresing sequence of random integer
126
* numbers from a given interval. The algorithm is taken literally
127
* from Jeffrey Scott Vitter: 'An Efficient Algorithm for Sequential
128
* Random Sampling', ACM Transactions on Mathematical Software, 13/1,
129
* 58--67. This method can be used for generating numbers from a
130
* \em very large interval, it is primilarly created for randomly
131
* selecting some edges from the sometimes huge set of possible edges
133
* \param res Pointer to an initialized vector, this will hold the
134
* result. It will be resized to the proper size.
135
* \param l The lower limit of the generation interval (inclusive).
136
* \param h The upper limit of the generation interval (inclusive).
137
* \param length The number of random integers to generate.
138
* \return Error code.
140
* Time complexity: according to the referenced paper, the expected
141
* running time is O(length).
144
int igraph_random_sample(igraph_vector_t *res, igraph_integer_t l, igraph_integer_t h,
145
igraph_integer_t length) {
146
igraph_real_t N=h-l+1;
147
igraph_real_t n=length;
150
igraph_real_t nreal=length;
151
igraph_real_t ninv=1.0/nreal;
152
igraph_real_t Nreal=N;
153
igraph_real_t Vprime;
154
igraph_real_t qu1=-n+1+N;
155
igraph_real_t qu1real=-nreal+1.0+Nreal;
156
igraph_real_t negalphainv=-13;
157
igraph_real_t threshold=-negalphainv*n;
160
igraph_vector_clear(res);
161
IGRAPH_CHECK(igraph_vector_reserve(res, length));
165
Vprime=exp(log(RNG_UNIF01())*ninv);
167
while (n>1 && threshold < N) {
169
igraph_real_t limit, t;
170
igraph_real_t negSreal, y1, y2, top, bottom;
171
igraph_real_t nmin1inv=1.0/(-1.0+nreal);
174
X=Nreal*(-Vprime+1.0);
177
if (S <qu1) { break; }
178
Vprime = exp(log(RNG_UNIF01())*ninv);
183
y1=exp(log(U*Nreal/qu1real)*nmin1inv);
184
Vprime=y1*(-X/Nreal+1.0)*(qu1real/(negSreal+qu1real));
185
if (Vprime <= 1.0) { break; }
193
bottom=-1.0+negSreal+Nreal;
196
for (t=-1+N; t>=limit; t--) {
201
if (Nreal/(-X+Nreal) >= y1*exp(log(y2)*nmin1inv)) {
202
Vprime=exp(log(RNG_UNIF01())*nmin1inv);
205
Vprime=exp(log(RNG_UNIF01())*ninv);
209
igraph_vector_push_back(res, l); /* allocated */
210
N=-S+(-1+N); Nreal=negSreal+(-1.0+Nreal);
211
n=-1+n; nreal=-1.0+nreal; ninv=nmin1inv;
212
qu1=-S+qu1; qu1real=negSreal+qu1real;
213
threshold=threshold+negalphainv;
217
retval=igraph_random_sample_alga(res, l, h, n);
222
igraph_vector_push_back(res, l); /* allocated */
235
* Mathlib : A C Library of Special Functions
236
* Copyright (C) 1998 Ross Ihaka
237
* Copyright (C) 2000 The R Development Core Team
238
* based on AS 111 (C) 1977 Royal Statistical Society
239
* and on AS 241 (C) 1988 Royal Statistical Society
241
* This program is free software; you can redistribute it and/or modify
242
* it under the terms of the GNU General Public License as published by
243
* the Free Software Foundation; either version 2 of the License, or
244
* (at your option) any later version.
246
* This program is distributed in the hope that it will be useful,
247
* but WITHOUT ANY WARRANTY; without even the implied warranty of
248
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
249
* GNU General Public License for more details.
251
* You should have received a copy of the GNU General Public License
252
* along with this program; if not, write to the Free Software
253
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
257
* double qnorm5(double p, double mu, double sigma,
258
* int lower_tail, int log_p)
259
* {qnorm (..) is synonymous and preferred inside R}
263
* Compute the quantile function for the normal distribution.
265
* For small to moderate probabilities, algorithm referenced
266
* below is used to obtain an initial approximation which is
267
* polished with a final Newton step.
269
* For very large arguments, an algorithm of Wichura is used.
273
* Beasley, J. D. and S. G. Springer (1977).
274
* Algorithm AS 111: The percentage points of the normal distribution,
275
* Applied Statistics, 26, 118-121.
277
* Wichura, M.J. (1988).
278
* Algorithm AS 241: The Percentage Points of the Normal Distribution.
279
* Applied Statistics, 37, 477-484.
283
* Mathlib : A C Library of Special Functions
284
* Copyright (C) 1998-2004 The R Development Core Team
286
* This program is free software; you can redistribute it and/or modify
287
* it under the terms of the GNU General Public License as published by
288
* the Free Software Foundation; either version 2 of the License, or
289
* (at your option) any later version.
291
* This program is distributed in the hope that it will be useful,
292
* but WITHOUT ANY WARRANTY; without even the implied warranty of
293
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
294
* GNU General Public License for more details.
296
* You should have received a copy of the GNU General Public License
297
* along with this program; if not, write to the Free Software
298
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
302
/* Private header file for use during compilation of Mathlib */
303
#ifndef MATHLIB_PRIVATE_H
304
#define MATHLIB_PRIVATE_H
307
# define ML_POSINF IGRAPH_INFINITY
308
# define ML_NEGINF -IGRAPH_INFINITY
309
# define ML_NAN IGRAPH_NAN
311
# define ML_POSINF (1.0 / 0.0)
312
# define ML_NEGINF ((-1.0) / 0.0)
313
# define ML_NAN (0.0 / 0.0)
316
#define ML_ERROR(x) /* nothing */
317
#define ML_UNDERFLOW (DBL_MIN * DBL_MIN)
318
#define ML_VALID(x) (!ISNAN(x))
323
/* argument out of domain */
325
/* value out of range */
327
/* process did not converge */
328
#define ME_PRECISION 8
329
/* does not have "full" precision */
330
#define ME_UNDERFLOW 16
331
/* and underflow occured (important for IEEE)*/
333
#define ML_ERR_return_NAN { ML_ERROR(ME_DOMAIN); return ML_NAN; }
335
/* Wilcoxon Rank Sum Distribution */
337
#define WILCOX_MAX 50
339
/* Wilcoxon Signed Rank Distribution */
341
#define SIGNRANK_MAX 50
343
/* Formerly private part of Mathlib.h */
345
/* always remap internal functions */
347
#define chebyshev_eval Rf_chebyshev_eval
348
#define chebyshev_init Rf_chebyshev_init
349
#define i1mach Rf_i1mach
350
#define gammalims Rf_gammalims
351
#define lfastchoose Rf_lfastchoose
352
#define lgammacor Rf_lgammacor
353
#define stirlerr Rf_stirlerr
355
/* Chebyshev Series */
357
int chebyshev_init(double*, int, double);
358
double chebyshev_eval(double, const double *, const int);
360
/* Gamma and Related Functions */
362
void gammalims(double*, double*);
363
double lgammacor(double); /* log(gamma) correction */
364
double stirlerr(double); /* Stirling expansion "error" */
366
double lfastchoose(double, double);
368
double bd0(double, double);
370
/* Consider adding these two to the API (Rmath.h): */
371
double dbinom_raw(double, double, double, double, int);
372
double dpois_raw (double, double, int);
373
double pnchisq_raw(double, double, double, double, double, int);
378
void bratio(double a, double b, double x, double y,
379
double *w, double *w1, int *ierr);
382
#endif /* MATHLIB_PRIVATE_H */
385
/* Utilities for `dpq' handling (density/probability/quantile) */
387
/* give_log in "d"; log_p in "p" & "q" : */
388
#define give_log log_p
391
#define R_D__0 (log_p ? ML_NEGINF : 0.) /* 0 */
392
#define R_D__1 (log_p ? 0. : 1.) /* 1 */
393
#define R_DT_0 (lower_tail ? R_D__0 : R_D__1) /* 0 */
394
#define R_DT_1 (lower_tail ? R_D__1 : R_D__0) /* 1 */
396
#define R_D_Lval(p) (lower_tail ? (p) : (1 - (p))) /* p */
397
#define R_D_Cval(p) (lower_tail ? (1 - (p)) : (p)) /* 1 - p */
399
#define R_D_val(x) (log_p ? log(x) : (x)) /* x in pF(x,..) */
400
#define R_D_qIv(p) (log_p ? exp(p) : (p)) /* p in qF(p,..) */
401
#define R_D_exp(x) (log_p ? (x) : exp(x)) /* exp(x) */
402
#define R_D_log(p) (log_p ? (p) : log(p)) /* log(p) */
403
#define R_D_Clog(p) (log_p ? log1p(-(p)) : (1 - (p)))/* [log](1-p) */
405
/* log(1-exp(x)): R_D_LExp(x) == (log1p(- R_D_qIv(x))) but even more stable:*/
406
#define R_D_LExp(x) (log_p ? R_Log1_Exp(x) : log1p(-x))
409
* #define R_DT_val(x) R_D_val(R_D_Lval(x))
410
* #define R_DT_Cval(x) R_D_val(R_D_Cval(x)) */
411
#define R_DT_val(x) (lower_tail ? R_D_val(x) : R_D_Clog(x))
412
#define R_DT_Cval(x) (lower_tail ? R_D_Clog(x) : R_D_val(x))
414
/*#define R_DT_qIv(p) R_D_Lval(R_D_qIv(p)) * p in qF ! */
415
#define R_DT_qIv(p) (log_p ? (lower_tail ? exp(p) : - expm1(p)) \
418
/*#define R_DT_CIv(p) R_D_Cval(R_D_qIv(p)) * 1 - p in qF */
419
#define R_DT_CIv(p) (log_p ? (lower_tail ? -expm1(p) : exp(p)) \
422
#define R_DT_exp(x) R_D_exp(R_D_Lval(x)) /* exp(x) */
423
#define R_DT_Cexp(x) R_D_exp(R_D_Cval(x)) /* exp(1 - x) */
425
#define R_DT_log(p) (lower_tail? R_D_log(p) : R_D_LExp(p))/* log(p) in qF */
426
#define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* log(1-p) in qF*/
427
#define R_DT_Log(p) (lower_tail? (p) : R_Log1_Exp(p))
428
/* == R_DT_log when we already "know" log_p == TRUE :*/
430
#define R_Q_P01_check(p) \
431
if ((log_p && p > 0) || \
432
(!log_p && (p < 0 || p > 1)) ) \
435
/* additions for density functions (C.Loader) */
436
#define R_D_fexp(f,x) (give_log ? -0.5*log(f)+(x) : exp(x)/sqrt(f))
437
#define R_D_forceint(x) floor((x) + 0.5)
438
#define R_D_nonint(x) (fabs((x) - floor((x)+0.5)) > 1e-7)
439
/* [neg]ative or [non int]eger : */
440
#define R_D_negInonint(x) (x < 0. || R_D_nonint(x))
442
#define R_D_nonint_check(x) \
443
if(R_D_nonint(x)) { \
444
MATHLIB_WARNING("non-integer x = %f", x); \
448
double igraph_qnorm5(double p, double mu, double sigma, int lower_tail, int log_p)
450
double p_, q, r, val;
453
if (ISNAN(p) || ISNAN(mu) || ISNAN(sigma))
454
return p + mu + sigma;
456
if (p == R_DT_0) return ML_NEGINF;
457
if (p == R_DT_1) return ML_POSINF;
460
if(sigma < 0) ML_ERR_return_NAN;
461
if(sigma == 0) return mu;
463
p_ = R_DT_qIv(p);/* real lower_tail prob. p */
466
/*-- use AS 241 --- */
467
/* double ppnd16_(double *p, long *ifault)*/
468
/* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
470
Produces the normal deviate Z corresponding to a given lower
471
tail area of P; Z is accurate to about 1 part in 10**16.
473
(original fortran code used PARAMETER(..) for the coefficients
474
and provided hash codes for checking them...)
476
if (fabs(q) <= .425) {/* 0.075 <= p <= 0.925 */
479
q * (((((((r * 2509.0809287301226727 +
480
33430.575583588128105) * r + 67265.770927008700853) * r +
481
45921.953931549871457) * r + 13731.693765509461125) * r +
482
1971.5909503065514427) * r + 133.14166789178437745) * r +
483
3.387132872796366608)
484
/ (((((((r * 5226.495278852854561 +
485
28729.085735721942674) * r + 39307.89580009271061) * r +
486
21213.794301586595867) * r + 5394.1960214247511077) * r +
487
687.1870074920579083) * r + 42.313330701600911252) * r + 1.);
489
else { /* closer than 0.075 from {0,1} boundary */
491
/* r = min(p, 1-p) < 0.075 */
493
r = R_DT_CIv(p);/* 1-p */
495
r = p_;/* = R_DT_Iv(p) ^= p */
497
r = sqrt(- ((log_p &&
498
((lower_tail && q <= 0) || (!lower_tail && q > 0))) ?
499
p : /* else */ log(r)));
500
/* r = sqrt(-log(r)) <==> min(p, 1-p) = exp( - r^2 ) */
502
if (r <= 5.) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
504
val = (((((((r * 7.7454501427834140764e-4 +
505
.0227238449892691845833) * r + .24178072517745061177) *
506
r + 1.27045825245236838258) * r +
507
3.64784832476320460504) * r + 5.7694972214606914055) *
508
r + 4.6303378461565452959) * r +
509
1.42343711074968357734)
511
1.05075007164441684324e-9 + 5.475938084995344946e-4) *
512
r + .0151986665636164571966) * r +
513
.14810397642748007459) * r + .68976733498510000455) *
514
r + 1.6763848301838038494) * r +
515
2.05319162663775882187) * r + 1.);
517
else { /* very close to 0 or 1 */
519
val = (((((((r * 2.01033439929228813265e-7 +
520
2.71155556874348757815e-5) * r +
521
.0012426609473880784386) * r + .026532189526576123093) *
522
r + .29656057182850489123) * r +
523
1.7848265399172913358) * r + 5.4637849111641143699) *
524
r + 6.6579046435011037772)
526
2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
527
r + 1.8463183175100546818e-5) * r +
528
7.868691311456132591e-4) * r + .0148753612908506148525)
529
* r + .13692988092273580531) * r +
530
.59983220655588793769) * r + 1.);
535
/* return (q >= 0.)? r : -r ;*/
537
return mu + sigma * val;
540
double fsign(double x, double y)
543
if (ISNAN(x) || ISNAN(y))
546
return ((y >= 0) ? fabs(x) : -fabs(x));
549
int imax2(int x, int y)
551
return (x < y) ? y : x;
554
int imin2(int x, int y)
556
return (x < y) ? x : y;
559
#ifdef HAVE_WORKING_ISFINITE
560
/* isfinite is defined in <math.h> according to C99 */
561
# define R_FINITE(x) isfinite(x)
562
#elif HAVE_WORKING_FINITE
563
/* include header needed to define finite() */
564
# ifdef HAVE_IEEE754_H
565
# include <ieee754.h> /* newer Linuxen */
567
# ifdef HAVE_IEEEFP_H
568
# include <ieeefp.h> /* others [Solaris], .. */
571
# define R_FINITE(x) finite(x)
573
# define R_FINITE(x) R_finite(x)
576
int R_finite(double x)
578
#ifdef HAVE_WORKING_ISFINITE
580
#elif HAVE_WORKING_FINITE
583
/* neither finite nor isfinite work. Do we really need the AIX exception? */
587
# elif defined(_MSC_VER)
590
return (!isnan(x) & (x != 1/0.0) & (x != -1.0/0.0));
595
int R_isnancpp(double x)
597
return (isnan(x)!=0);
601
int R_isnancpp(double); /* in arithmetic.c */
602
# define ISNAN(x) R_isnancpp(x)
604
# define ISNAN(x) (isnan(x)!=0)
607
double igraph_norm_rand(void) {
611
#define BIG 134217728 /* 2^27 */
612
/* unif_rand() alone is not of high enough precision */
614
u1 = (int)(BIG*u1) + RNG_UNIF(0,1);
615
return igraph_qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
619
* Mathlib : A C Library of Special Functions
620
* Copyright (C) 1998 Ross Ihaka
621
* Copyright (C) 2000-2002 the R Development Core Team
623
* This program is free software; you can redistribute it and/or modify
624
* it under the terms of the GNU General Public License as published by
625
* the Free Software Foundation; either version 2 of the License, or
626
* (at your option) any later version.
628
* This program is distributed in the hope that it will be useful,
629
* but WITHOUT ANY WARRANTY; without even the implied warranty of
630
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
631
* GNU General Public License for more details.
633
* You should have received a copy of the GNU General Public License
634
* along with this program; if not, write to the Free Software
635
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
640
* double exp_rand(void);
644
* Random variates from the standard exponential distribution.
648
* Ahrens, J.H. and Dieter, U. (1972).
649
* Computer methods for sampling from the exponential and
650
* normal distributions.
651
* Comm. ACM, 15, 873-882.
654
double igraph_exp_rand(void)
656
/* q[k-1] = sum(log(2)^k / k!) k=1,..,n, */
657
/* The highest n (here 8) is determined by q[n-1] = 1.0 */
658
/* within standard precision */
678
double a, u, ustar, umin;
682
/* precaution if u = 0 is ever returned */
684
while(u <= 0.0 || u >= 1.0) u = RNG_UNIF01();
697
ustar = RNG_UNIF01();
700
ustar = RNG_UNIF01();
705
return a + umin * q[0];
709
* Mathlib : A C Library of Special Functions
710
* Copyright (C) 1998 Ross Ihaka
711
* Copyright (C) 2000-2001 The R Development Core Team
713
* This program is free software; you can redistribute it and/or modify
714
* it under the terms of the GNU General Public License as published by
715
* the Free Software Foundation; either version 2 of the License, or
716
* (at your option) any later version.
718
* This program is distributed in the hope that it will be useful,
719
* but WITHOUT ANY WARRANTY; without even the implied warranty of
720
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
721
* GNU General Public License for more details.
723
* You should have received a copy of the GNU General Public License
724
* along with this program; if not, write to the Free Software
725
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
730
* double rpois(double lambda)
734
* Random variates from the Poisson distribution.
738
* Ahrens, J.H. and Dieter, U. (1982).
739
* Computer generation of Poisson deviates
740
* from modified normal distributions.
741
* ACM Trans. Math. Software 8, 163-179.
746
#define a2 -0.2500068
748
#define a4 -0.1661269
750
#define a6 -0.1384794
753
#define one_7 0.1428571428571428571
754
#define one_12 0.0833333333333333333
755
#define one_24 0.0416666666666666667
757
#define repeat for(;;)
761
#define M_1_SQRT_2PI 0.398942280401432677939946059934 /* 1/sqrt(2pi) */
763
double igraph_rpois(double mu)
765
/* Factorial Table (0:9)! */
766
const double fact[10] =
768
1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880.
771
/* These are static --- persistent between calls for same mu : */
774
static double b1, b2, c, c0, c1, c2, c3;
775
static double pp[36], p0, p, q, s, d, omega;
776
static double big_l;/* integer "w/o overflow" */
777
static double muprev = 0., muprev2 = 0.;/*, muold = 0.*/
779
/* Local Vars [initialize some for -Wall]: */
780
double del, difmuk= 0., E= 0., fk= 0., fx, fy, g, px, py, t, u= 0., v, x;
782
int k, kflag, big_mu, new_big_mu = FALSE;
794
if (!(big_mu && mu == muprev)) {/* maybe compute new persistent par.s */
798
/* Case A. (recalculation of s,d,l because mu has changed):
799
* The poisson probabilities pk exceed the discrete normal
800
* probabilities fk whenever k >= m(mu).
805
big_l = floor(mu - 1.1484);
806
/* = an upper bound to m(mu) for all mu >= 10.*/
808
else { /* Small mu ( < 10) -- not using normal approx. */
810
/* Case B. (start new table and calculate p0 if necessary) */
812
/*muprev = 0.;-* such that next time, mu != muprev ..*/
815
m = imax2(1, (int) mu);
816
l = 0; /* pp[] is already ok up to pp[l] */
817
q = p0 = p = exp(-mu);
821
/* Step U. uniform sample for inversion method */
826
/* Step T. table comparison until the end pp[l] of the
827
pp-table of cumulative poisson probabilities
828
(0.458 > ~= pp[9](= 0.45792971447) for mu=10 ) */
830
for (k = (u <= 0.458) ? 1 : imin2(l, m); k <= l; k++)
833
if (l == 35) /* u > pp[35] */
836
/* Step C. creation of new poisson
837
probabilities p[l..] and their cumulatives q =: pp[k] */
839
for (k = l; k <= 35; k++) {
852
} /* end {initialize persistent vars} */
854
/* Only if mu >= 10 : ----------------------- */
856
/* Step N. normal sample */
857
g = mu + s * igraph_norm_rand();/* norm_rand() ~ N(0,1), standard normal */
861
/* Step I. immediate acceptance if pois is large enough */
864
/* Step S. squeeze acceptance */
867
u = RNG_UNIF01(); /* ~ U(0,1) - sample */
868
if (d * u >= difmuk * difmuk * difmuk)
872
/* Step P. preparations for steps Q and H.
873
(recalculations of parameters if necessary) */
875
if (new_big_mu || mu != muprev2) {
876
/* Careful! muprev2 is not always == muprev
877
because one might have exited in step I or S
880
omega = M_1_SQRT_2PI / s;
881
/* The quantities b1, b2, c3, c2, c1, c0 are for the Hermite
882
* approximations to the discrete normal probabilities fk. */
886
c3 = one_7 * b1 * b2;
888
c1 = b1 - 6. * b2 + 45. * c3;
889
c0 = 1. - b1 + 3. * b2 - 15. * c3;
890
c = 0.1069 / mu; /* guarantees majorization by the 'hat'-function. */
894
/* 'Subroutine' F is called (kflag=0 for correct return) */
901
/* Step E. Exponential Sample */
903
E = igraph_exp_rand(); /* ~ Exp(1) (standard exponential) */
905
/* sample t from the laplace 'hat'
906
(if t <= -0.6744 then pk < fk for all mu >= 10.) */
907
u = 2 * RNG_UNIF01() - 1.;
908
t = 1.8 + fsign(E, u);
910
pois = floor(mu + s * t);
914
/* 'subroutine' F is called (kflag=1 for correct return) */
917
Step_F: /* 'subroutine' F : calculation of px,py,fx,fy. */
919
if (pois < 10) { /* use factorials from table fact[] */
921
py = pow(mu, pois) / fact[(int)pois];
924
/* Case pois >= 10 uses polynomial approximation
925
a0-a7 for accuracy when advisable */
927
del = del * (1. - 4.8 * del * del);
930
px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) *
931
v + a3) * v + a2) * v + a1) * v + a0)
934
px = fk * log(1. + v) - difmuk - del;
935
py = M_1_SQRT_2PI / sqrt(fk);
937
x = (0.5 - difmuk) / s;
940
fy = omega * (((c3 * x + c2) * x + c1) * x + c0);
942
/* Step H. Hat acceptance (E is repeated on rejection) */
943
if (c * fabs(u) <= py * exp(px + E) - fy * exp(fx + E))
946
/* Step Q. Quotient acceptance (rare case) */
947
if (fy - u * fy <= py * exp(px - fx))
954
double igraph_rgeom(double p) {
955
if (ISNAN(p) || p <= 0 || p > 1) ML_ERR_return_NAN;
957
return igraph_rpois(igraph_exp_rand() * ((1 - p) / p));
960
/* This is from nmath/rbinom.c */
962
#define repeat for(;;)
964
double igraph_rbinom(double nin, double pp)
966
/* FIXME: These should become THREAD_specific globals : */
968
static double c, fm, npq, p1, p2, p3, p4, qn;
969
static double xl, xll, xlr, xm, xr;
971
static double psave = -1.0;
972
static int nsave = -1;
975
double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
976
double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
979
if (!R_FINITE(nin)) ML_ERR_return_NAN;
980
n = floor(nin + 0.5);
981
if (n != nin) ML_ERR_return_NAN;
984
/* n=0, p=0, p=1 are not errors <TSL>*/
985
n < 0 || pp < 0. || pp > 1.) ML_ERR_return_NAN;
987
if (n == 0 || pp == 0.) return 0;
988
if (pp == 1.) return n;
990
p = fmin(pp, 1. - pp);
996
/* Setup, perform only when parameters change [using static (globals): */
998
/* FIXING: Want this thread safe
999
-- use as little (thread globals) as possible
1001
if (pp != psave || n != nsave) {
1005
/* inverse cdf logic for mean less than 30 */
1006
qn = pow(q, (double) n);
1013
p1 = (int)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
1017
c = 0.134 + 20.5 / (15.3 + fm);
1018
al = (ffm - xl) / (ffm - xl * p);
1019
xll = al * (1.0 + 0.5 * al);
1020
al = (xr - ffm) / (xr * q);
1021
xlr = al * (1.0 + 0.5 * al);
1022
p2 = p1 * (1.0 + c + c);
1026
} else if (n == nsave) {
1031
/*-------------------------- np = n*p >= 30 : ------------------- */
1033
u = RNG_UNIF01() * p4;
1035
/* triangular region */
1037
ix = xm - p1 * v + u;
1040
/* parallelogram region */
1042
x = xl + (u - p1) / c;
1043
v = v * c + 1.0 - fabs(xm - x) / p1;
1044
if (v > 1.0 || v <= 0.)
1048
if (u > p3) { /* right tail */
1049
ix = xr - log(v) / xlr;
1052
v = v * (u - p3) * xlr;
1053
} else {/* left tail */
1054
ix = xl + log(v) / xll;
1057
v = v * (u - p2) * xll;
1060
/* determine appropriate way to perform accept/reject test */
1062
if (k <= 20 || k >= npq / 2 - 1) {
1063
/* explicit evaluation */
1066
for (i = m + 1; i <= ix; i++)
1068
} else if (m != ix) {
1069
for (i = ix + 1; i <= m; i++)
1075
/* squeezing using upper and lower bounds on log(f(x)) */
1076
amaxp = (k / npq) * ((k * (k / 3. + 0.625) + 0.1666666666666) / npq + 0.5);
1077
ynorm = -k * k / (2.0 * npq);
1079
if (alv < ynorm - amaxp)
1081
if (alv <= ynorm + amaxp) {
1082
/* stirling's formula to machine accuracy */
1083
/* for the final acceptance/rejection test */
1092
if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.)
1099
/*---------------------- np = n*p < 30 : ------------------------- */
1123
/**********************************************************
1124
* Testing purposes *
1125
*********************************************************/
1131
/* for (i=0; i<100; i++) { */
1132
/* printf("%i ", RNG_INTEGER(0,10)); */
1136
/* for (i=0; i<100; i++) { */
1137
/* printf("%f ", RNG_UNIF(0,1)); */
1141
/* for (i=0; i<100; i++) { */
1142
/* printf("%f ", RNG_NORMAL(0,5)); */
added
removed
src/random.c
