Mercurial > octave

annotate liboctave/numeric/randpoisson.cc @ 23475:d691ed308237

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Clean up #includes in liboctave/numeric directory. * build-aux/mk-opts.pl: Change Perl to generate "" around local include libraries rather than <>. Include "lo-math.h" rather than <cmath>. * CollocWt.cc, DAERTFunc.h, DASPK.cc, DASPK.h, DASRT.cc, DASRT.h, DASSL.cc, DASSL.h, DET.h, EIG.cc, EIG.h, LSODE.cc, LSODE.h, ODE.h, ODES.cc, ODESFunc.h, Quad.cc, aepbalance.cc, base-de.h, base-min.h, bsxfun-decl.h, bsxfun-defs.cc, bsxfun.h, chol.cc, eigs-base.cc, fEIG.cc, fEIG.h, gepbalance.cc, gsvd.cc, hess.cc, lo-blas-proto.h, lo-lapack-proto.h, lo-mappers.cc, lo-mappers.h, lo-qrupdate-proto.h, lo-slatec-proto.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, randgamma.cc, randpoisson.cc, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-dmsolve.cc, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc: Rationalize #includes. Use forward declarations of just classes where possible. Reformat some long lines < 80 characters. Reformat some comments for readabliity. * mx-inlines.cc: Rationalize #includes for this file in liboctave/operators used by many in liboctave/numeric.
author Rik <rik@octave.org>
date Tue, 09 May 2017 08:46:07 -0700
parents 092078913d54
children 864ac8dabeff
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
1 /*
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
2
23219
3ac9f9ecfae5 maint: Update copyright dates.
John W. Eaton <jwe@octave.org>
parents: 23083
diff changeset
3 Copyright (C) 2006-2017 John W. Eaton
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
4
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
5 This file is part of Octave.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
6
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
8 under the terms of the GNU General Public License as published by
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
9 the Free Software Foundation; either version 3 of the License, or
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
10 (at your option) any later version.
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
11
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
12 Octave is distributed in the hope that it will be useful, but
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
13 WITHOUT ANY WARRANTY; without even the implied warranty of
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
15 GNU General Public License for more details.
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
16
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
17 You should have received a copy of the GNU General Public License
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
18 along with Octave; see the file COPYING. If not, see
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
19 <http://www.gnu.org/licenses/>.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
20
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
21 */
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
22
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
23 /* Original version written by Paul Kienzle distributed as free
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
24 software in the in the public domain. */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
25
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
26 /* Needs the following defines:
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
27 * NAN: value to return for Not-A-Number
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
28 * RUNI: uniform generator on (0,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
29 * RNOR: normal generator
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
30 * LGAMMA: log gamma function
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
31 * INFINITE: function to test whether a value is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
32 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
33
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
34 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21202
diff changeset
35 # include "config.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
36 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
37
23475
d691ed308237 maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents: 23220
diff changeset
38 #include <cstddef>
d691ed308237 maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents: 23220
diff changeset
39
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
40 #include "f77-fcn.h"
7231
2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe]
jwe
parents: 7019
diff changeset
41 #include "lo-error.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
42 #include "lo-ieee.h"
7231
2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe]
jwe
parents: 7019
diff changeset
43 #include "lo-math.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
44 #include "randmtzig.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
45 #include "randpoisson.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
46
6096
40212bfaed1d [project @ 2006-10-26 01:57:49 by jwe]
jwe
parents: 5775
diff changeset
47 #undef INFINITE
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
48 #define INFINITE lo_ieee_isinf
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
49 #define RUNI oct_randu()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
50 #define RNOR oct_randn()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
51 #define LGAMMA xlgamma
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
52
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
53 static double
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
54 xlgamma (double x)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
55 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
56 double result;
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21723
diff changeset
57 #if defined (HAVE_LGAMMA)
7601
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
58 result = lgamma (x);
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
59 #else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
60 double sgngam;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
61
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
62 if (lo_ieee_isnan (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
63 result = x;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
64 else if (x <= 0 || lo_ieee_isinf (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
65 result = octave_Inf;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
66 else
22950
f2fae9bd5cec * randpoisson.cc: Fix dlgams call to match prototype.
John W. Eaton <jwe@octave.org>
parents: 22755
diff changeset
67 F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
7601
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
68 #endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
69 return result;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
70 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
71
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
72 /* ---- pprsc.c from Stadloeber's winrand --- */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
73
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
74 /* flogfak(k) = ln(k!) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
75 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
76 flogfak (double k)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
77 {
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
78 #define C0 9.18938533204672742e-01
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
79 #define C1 8.33333333333333333e-02
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
80 #define C3 -2.77777777777777778e-03
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
81 #define C5 7.93650793650793651e-04
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
82 #define C7 -5.95238095238095238e-04
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
83
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
84 static double logfak[30L] =
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
85 {
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
86 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
87 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
89 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
90 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
91 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
92 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
93 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
94 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
95 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
96 };
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
97
21562
6c2fd62db1f7 maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents: 21301
diff changeset
98 double r, rr;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
99
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
100 if (k >= 30.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
101 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
102 r = 1.0 / k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
103 rr = r * r;
23475
d691ed308237 maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents: 23220
diff changeset
104 return ((k + 0.5)*std::log (k) - k + C0
d691ed308237 maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents: 23220
diff changeset
105 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
106 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
107 else
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
108 return (logfak[static_cast<int> (k)]);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
109 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
110
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
111 /******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112 * *
21562
6c2fd62db1f7 maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents: 21301
diff changeset
113 * Poisson Distribution - Patchwork Rejection/Inversion *
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
114 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
115 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
116 * *
21562
6c2fd62db1f7 maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents: 21301
diff changeset
117 * For parameter my < 10, Tabulated Inversion is applied. *
6c2fd62db1f7 maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents: 21301
diff changeset
118 * For my >= 10, Patchwork Rejection is employed: *
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
119 * The area below the histogram function f(x) is rearranged in *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
120 * its body by certain point reflections. Within a large center *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
121 * interval variates are sampled efficiently by rejection from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 * uniform hats. Rectangular immediate acceptance regions speed *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
123 * up the generation. The remaining tails are covered by *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
124 * exponential functions. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
125 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
126 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
127 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
128 * FUNCTION : - pprsc samples a random number from the Poisson *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
129 * distribution with parameter my > 0. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
130 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
131 * continuous and discrete unimodal distributions, *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
132 * Doctoral Dissertation, 156 pp., Technical *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
133 * University Graz, Austria. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
134 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
135 * unsigned long integer *seed. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
136 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
137 * Implemented by H. Zechner, January 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
138 * Revised by F. Niederl, July 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
140 ******************************************************************/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
141
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
142 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
143 f (double k, double l_nu, double c_pm)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
144 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
145 return exp (k * l_nu - flogfak (k) - c_pm);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
146 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
147
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
148 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
149 pprsc (double my)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
150 {
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
151 static double my_last = -1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
152 static double m, k2, k4, k1, k5;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
153 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
154 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
155 double Dk, X, Y;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
156 double Ds, U, V, W;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
157
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
158 if (my != my_last)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
159 { /* set-up */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
160 my_last = my;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
161 /* approximate deviation of reflection points k2, k4 from my - 1/2 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
162 Ds = sqrt (my + 0.25);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
163
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
164 /* mode m, reflection points k2 and k4, and points k1 and k5, */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
165 /* which delimit the centre region of h(x) */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
166 m = std::floor (my);
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
167 k2 = ceil (my - 0.5 - Ds);
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
168 k4 = std::floor (my - 0.5 + Ds);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
169 k1 = k2 + k2 - m + 1L;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
170 k5 = k4 + k4 - m;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
171
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
172 /* range width of the critical left and right centre region */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
173 dl = (k2 - k1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
174 dr = (k5 - k4);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
175
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
176 /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
177 r1 = my / k1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
178 r2 = my / k2;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
179 r4 = my / (k4 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
180 r5 = my / (k5 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
181
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
182 /* reciprocal values of the scale parameters of exp. tail envelope */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
183 ll = std::log (r1); /* expon. tail left */
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
184 lr = -std::log (r5); /* expon. tail right*/
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
185
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 /* Poisson constants, necessary for computing function values f(k) */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
187 l_my = std::log (my);
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
188 c_pm = m * l_my - flogfak (m);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
189
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
190 /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
191 f2 = f (k2, l_my, c_pm);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
192 f4 = f (k4, l_my, c_pm);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
193 f1 = f (k1, l_my, c_pm);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
194 f5 = f (k5, l_my, c_pm);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
195
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
196 /* area of the two centre and the two exponential tail regions */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
197 /* area of the two immediate acceptance regions between k2, k4 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
198 p1 = f2 * (dl + 1.0); /* immed. left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
199 p2 = f2 * dl + p1; /* centre left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
200 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
201 p4 = f4 * dr + p3; /* centre right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
202 p5 = f1 / ll + p4; /* exp. tail left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
203 p6 = f5 / lr + p5; /* exp. tail right*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
204 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
205
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
206 for (;;)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
207 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
208 /* generate uniform number U -- U(0, p6) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
209 /* case distinction corresponding to U */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
210 if ((U = RUNI * p6) < p2)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
211 { /* centre left */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
212
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
213 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
214 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
215 if ((V = U - p1) < 0.0) return (k2 + std::floor (U/f2));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
216 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
217 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
218 if ((W = V / dl) < f1 ) return (k1 + std::floor (V/f1));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
219
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
220 /* computation of candidate X < k2, and its counterpart Y > k2 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
221 /* either squeeze-acceptance of X or acceptance-rejection of Y */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
222 Dk = std::floor (dl * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
223 if (W <= f2 - Dk * (f2 - f2/r2))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
224 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
225 return (k2 - Dk); /* X = k2 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
226 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
227 if ((V = f2 + f2 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
228 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
229 Y = k2 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
230 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
231 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
232 return (Y); /* Y = k2 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
233 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
234 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
235 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
236 X = k2 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
237 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
238 else if (U < p4)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
239 { /* centre right */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
240 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
241 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
242 if ((V = U - p3) < 0.0) return (k4 - std::floor ((U - p2)/f4));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
243 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
244 R4 = [k4+1, k5+1)*[0, f5) */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
245 if ((W = V / dr) < f5 ) return (k5 - std::floor (V/f5));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
246
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
247 /* computation of candidate X > k4, and its counterpart Y < k4 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
248 /* either squeeze-acceptance of X or acceptance-rejection of Y */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
249 Dk = std::floor (dr * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
250 if (W <= f4 - Dk * (f4 - f4*r4))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
251 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
252 return (k4 + Dk); /* X = k4 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
253 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
254 if ((V = f4 + f4 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
255 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
256 Y = k4 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
257 if (V <= f4 + Dk * (1.0 - f4)/ dr)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
258 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
259 return (Y); /* Y = k4 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
260 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
261 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
262 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
263 X = k4 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
264 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
265 else
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
266 {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
267 W = RUNI;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
268 if (U < p5)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
269 { /* expon. tail left */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
270 Dk = std::floor (1.0 - std::log (W)/ll);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
271 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
272 W *= (U - p4) * ll; /* W -- U(0, h(x)) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
273 if (W <= f1 - Dk * (f1 - f1/r1))
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
274 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
275 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
276 else
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
277 { /* expon. tail right*/
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
278 Dk = std::floor (1.0 - std::log (W)/lr);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
279 X = k5 + Dk; /* X >= k5 + 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
280 W *= (U - p5) * lr; /* W -- U(0, h(x)) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
281 if (W <= f5 - Dk * (f5 - f5*r5))
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
282 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
283 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
284 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
285
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
286 /* acceptance-rejection test of candidate X from the original area */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
287 /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
288 /* log f(X) = (X - m)*log(my) - log X! + log m! */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
289 if (std::log (W) <= X * l_my - flogfak (X) - c_pm) return (X);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
290 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
291 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
292 /* ---- pprsc.c end ------ */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
293
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
294 /* The remainder of the file is by Paul Kienzle */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
295
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
296 /* Given uniform u, find x such that CDF(L,x)==u. Return x. */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
297 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
298 poisson_cdf_lookup (double lambda, double *p, size_t n)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
299 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
300 /* Table size is predicated on the maximum value of lambda
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
301 * we want to store in the table, and the maximum value of
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
302 * returned by the uniform random number generator on [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
303 * With lambda==10 and u_max = 1 - 1/(2^32+1), we
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
304 * have poisson_pdf(lambda,36) < 1-u_max. If instead our
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
305 * generator uses more bits of mantissa or returns a value
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
306 * in the range [0,1], then for lambda==10 we need a table
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
307 * size of 46 instead. For long doubles, the table size
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
308 * will need to be longer still. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
309 #define TABLESIZE 46
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
310 double t[TABLESIZE];
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
311
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
312 /* Precompute the table for the u up to and including 0.458.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
313 * We will almost certainly need it. */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
314 int intlambda = static_cast<int> (std::floor (lambda));
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
315 double P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
316 int tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
317 size_t i = n;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
318
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
319 t[0] = P = exp (-lambda);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
320 for (tableidx = 1; tableidx <= intlambda; tableidx++)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
321 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
322 P = P*lambda/static_cast<double> (tableidx);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
323 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
324 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
325
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
326 while (i-- > 0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
327 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
328 double u = RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
329
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
330 /* If u > 0.458 we know we can jump to floor(lambda) before
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
331 * comparing (this observation is based on Stadlober's winrand
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
332 * code). For lambda >= 1, this will be a win. Lambda < 1
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
333 * is already fast, so adding an extra comparison is not a
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
334 * problem. */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
335 int k = (u > 0.458 ? intlambda : 0);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
336
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
337 /* We aren't using a for loop here because when we find the
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
338 * right k we want to jump to the next iteration of the
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
339 * outer loop, and the continue statement will only work for
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
340 * the inner loop. */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
341 nextk:
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
342 if (u <= t[k])
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
343 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
344 p[i] = static_cast<double> (k);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
345 continue;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
346 }
21168
26f85aa072de maint: Replace instances of goto in liboctave where convenient.
Rik <rik@octave.org>
parents: 21066
diff changeset
347 if (++k < tableidx)
26f85aa072de maint: Replace instances of goto in liboctave where convenient.
Rik <rik@octave.org>
parents: 21066
diff changeset
348 goto nextk;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
349
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
350 /* We only need high values of the table very rarely so we
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
351 * don't automatically compute the entire table. */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
352 while (tableidx < TABLESIZE)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
353 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
354 P = P*lambda/static_cast<double> (tableidx);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
355 t[tableidx] = t[tableidx-1] + P;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
356 /* Make sure we converge to 1.0 just in case u is uniform
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
357 * on [0,1] rather than [0,1). */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
358 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
359 tableidx++;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
360 if (u <= t[tableidx-1]) break;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
361 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
362
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
363 /* We are assuming that the table size is big enough here.
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
364 * This should be true even if RUNI is returning values in
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
365 * the range [0,1] rather than [0,1). */
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
366 p[i] = static_cast<double> (tableidx-1);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
367 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
368 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
369
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
370 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
371 poisson_cdf_lookup_float (double lambda, float *p, size_t n)
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
372 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
373 double t[TABLESIZE];
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
374
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
375 /* Precompute the table for the u up to and including 0.458.
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
376 * We will almost certainly need it. */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
377 int intlambda = static_cast<int> (std::floor (lambda));
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
378 double P;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
379 int tableidx;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
380 size_t i = n;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
381
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
382 t[0] = P = exp (-lambda);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
383 for (tableidx = 1; tableidx <= intlambda; tableidx++)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
384 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
385 P = P*lambda/static_cast<double> (tableidx);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
386 t[tableidx] = t[tableidx-1] + P;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
387 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
388
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
389 while (i-- > 0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
390 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
391 double u = RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
392 int k = (u > 0.458 ? intlambda : 0);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
393 nextk:
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
394 if (u <= t[k])
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
395 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
396 p[i] = static_cast<float> (k);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
397 continue;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
398 }
21168
26f85aa072de maint: Replace instances of goto in liboctave where convenient.
Rik <rik@octave.org>
parents: 21066
diff changeset
399 if (++k < tableidx)
26f85aa072de maint: Replace instances of goto in liboctave where convenient.
Rik <rik@octave.org>
parents: 21066
diff changeset
400 goto nextk;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
401
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
402 while (tableidx < TABLESIZE)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
403 {
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
404 P = P*lambda/static_cast<double> (tableidx);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
405 t[tableidx] = t[tableidx-1] + P;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
406 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
407 tableidx++;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
408 if (u <= t[tableidx-1]) break;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
409 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
410
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21568
diff changeset
411 p[i] = static_cast<float> (tableidx-1);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
412 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
413 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
414
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
415 /* From Press, et al., Numerical Recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
416 static void
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
417 poisson_rejection (double lambda, double *p, size_t n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
418 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
419 double sq = sqrt (2.0*lambda);
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
420 double alxm = std::log (lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
421 double g = lambda*alxm - LGAMMA(lambda+1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
422 size_t i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
423
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
424 for (i = 0; i < n; i++)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
425 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
426 double y, em, t;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
427 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
428 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
429 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
430 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
431 y = tan (M_PI*RUNI);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
432 em = sq * y + lambda;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
433 } while (em < 0.0);
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
434 em = std::floor (em);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
435 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
436 } while (RUNI > t);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
437 p[i] = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
438 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
439 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
440
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
441 /* From Press, et al., Numerical Recipes */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
442 static void
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
443 poisson_rejection_float (double lambda, float *p, size_t n)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
444 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
445 double sq = sqrt (2.0*lambda);
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
446 double alxm = std::log (lambda);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
447 double g = lambda*alxm - LGAMMA(lambda+1.0);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
448 size_t i;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
449
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
450 for (i = 0; i < n; i++)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
451 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
452 double y, em, t;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
453 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
454 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
455 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
456 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
457 y = tan (M_PI*RUNI);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
458 em = sq * y + lambda;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
459 } while (em < 0.0);
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
460 em = std::floor (em);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
461 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
462 } while (RUNI > t);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
463 p[i] = em;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
464 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
465 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
466
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
467 /* The cutoff of L <= 1e8 in the following two functions before using
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
468 * the normal approximation is based on:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
469 * > L=1e8; x=floor(linspace(0,2*L,1000));
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
470 * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
21568
3d60ed163b70 maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents: 21562
diff changeset
471 * ans = 1.1376e-28
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
472 * For L=1e7, the max is around 1e-9, which is within the step size of RUNI.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
473 * For L>1e10 the pprsc function breaks down, as I saw from the histogram
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
474 * of a large sample, so 1e8 is both small enough and large enough. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
475
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
476 /* Generate a set of poisson numbers with the same distribution */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
477 void
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
478 oct_fill_randp (double L, octave_idx_type n, double *p)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
479 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
480 octave_idx_type i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
481 if (L < 0.0 || INFINITE(L))
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
482 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
483 for (i=0; i<n; i++)
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
484 p[i] = octave::numeric_limits<double>::NaN ();
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
485 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
486 else if (L <= 10.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
487 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
488 poisson_cdf_lookup (L, p, n);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
489 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
490 else if (L <= 1e8)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
491 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
492 for (i=0; i<n; i++)
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
493 p[i] = pprsc (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
494 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
495 else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
496 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
497 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
498 const double sqrtL = sqrt (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
499 for (i = 0; i < n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
500 {
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
501 p[i] = std::floor (RNOR*sqrtL + L + 0.5);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
502 if (p[i] < 0.0)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
503 p[i] = 0.0; /* will probably never happen */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
504 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
505 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
506 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
507
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
508 /* Generate one poisson variate */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
509 double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
510 oct_randp (double L)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
511 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
512 double ret;
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
513 if (L < 0.0) ret = octave::numeric_limits<double>::NaN ();
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
514 else if (L <= 12.0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
515 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
516 /* From Press, et al. Numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
517 double g = exp (-L);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
518 int em = -1;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
519 double t = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
520 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
521 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
522 ++em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
523 t *= RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
524 } while (t > g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
525 ret = em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
526 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
527 else if (L <= 1e8)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
528 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
529 /* numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
530 poisson_rejection (L, &ret, 1);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
531 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
532 else if (INFINITE(L))
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
533 {
21066
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 19697
diff changeset
534 /* FIXME: R uses NaN, but the normal approximation suggests that
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 19697
diff changeset
535 * limit should be Inf. Which is correct? */
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
536 ret = octave::numeric_limits<double>::NaN ();
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
537 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
538 else
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
539 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
540 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
541 ret = std::floor (RNOR*sqrt (L) + L + 0.5);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
542 if (ret < 0.0) ret = 0.0; /* will probably never happen */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
543 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
544 return ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
545 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
546
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
547 /* Generate a set of poisson numbers with the same distribution */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
548 void
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
549 oct_fill_float_randp (float FL, octave_idx_type n, float *p)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
550 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
551 double L = FL;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
552 octave_idx_type i;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
553 if (L < 0.0 || INFINITE(L))
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
554 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
555 for (i=0; i<n; i++)
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
556 p[i] = octave::numeric_limits<double>::NaN ();
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
557 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
558 else if (L <= 10.0)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
559 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
560 poisson_cdf_lookup_float (L, p, n);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
561 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
562 else if (L <= 1e8)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
563 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
564 for (i=0; i<n; i++)
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
565 p[i] = pprsc (L);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
566 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
567 else
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
568 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
569 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
570 const double sqrtL = sqrt (L);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
571 for (i = 0; i < n; i++)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
572 {
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
573 p[i] = std::floor (RNOR*sqrtL + L + 0.5);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
574 if (p[i] < 0.0)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
575 p[i] = 0.0; /* will probably never happen */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
576 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
577 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
578 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
579
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
580 /* Generate one poisson variate */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
581 float
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
582 oct_float_randp (float FL)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
583 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
584 double L = FL;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
585 float ret;
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
586 if (L < 0.0) ret = octave::numeric_limits<float>::NaN ();
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
587 else if (L <= 12.0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
588 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
589 /* From Press, et al. Numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
590 double g = exp (-L);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
591 int em = -1;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
592 double t = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
593 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
594 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
595 ++em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
596 t *= RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
597 } while (t > g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
598 ret = em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
599 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
600 else if (L <= 1e8)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
601 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
602 /* numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
603 poisson_rejection_float (L, &ret, 1);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
604 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
605 else if (INFINITE(L))
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
606 {
21066
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 19697
diff changeset
607 /* FIXME: R uses NaN, but the normal approximation suggests that
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 19697
diff changeset
608 * limit should be Inf. Which is correct? */
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
609 ret = octave::numeric_limits<float>::NaN ();
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
610 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
611 else
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
612 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
613 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
21942
aab79a1885cc limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
614 ret = std::floor (RNOR*sqrt (L) + L + 0.5);
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
615 if (ret < 0.0) ret = 0.0; /* will probably never happen */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
616 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
617 return ret;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
618 }