Mercurial > hg > octave-nkf

annotate liboctave/numeric/randpoisson.c @ 18518:694e8f0b3862 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
freqz: Make frequency response plot visually compatible with Matlab (bug #41464) * freqz_plot.m: Eliminate the extraneous "pass band" subplot for visual compatibility with Matlab. Use y-axis labels instead of legends. Always autoscale the y-axis. Label the frequency axis in units of normalized radians or Hz depending on a new optional argument. * freqz.m: Tell freqz_plot whether the frequency vector is in normalized radians or Hz.
author Mike Miller <mtmiller@ieee.org>
date Thu, 13 Feb 2014 00:00:10 -0500
parents 49a5a4be04a1
children 4197fc428c7d
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
17744
d63878346099 maint: Update copyright notices for release.
John W. Eaton <jwe@octave.org>
parents: 15851
diff changeset
3 Copyright (C) 2006-2013 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
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
8 under the terms of the GNU General Public License as published by the
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
9 Free Software Foundation; either version 3 of the License, or (at your
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
10 option) any later version.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
11
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
15 for more details.
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)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
35 #include <config.h>
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
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
38 #include <stdio.h>
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
39
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
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
47 #undef NAN
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
48 #define NAN octave_NaN
6096
40212bfaed1d [project @ 2006-10-26 01:57:49 by jwe]
jwe
parents: 5775
diff changeset
49 #undef INFINITE
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
50 #define INFINITE lo_ieee_isinf
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
51 #define RUNI oct_randu()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
52 #define RNOR oct_randn()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
53 #define LGAMMA xlgamma
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
54
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
55 F77_RET_T
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
56 F77_FUNC (dlgams, DLGAMS) (const double *, double *, double *);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
57
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
58 static double
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
59 xlgamma (double x)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
60 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
61 double result;
7601
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
62 #ifdef HAVE_LGAMMA
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
63 result = lgamma (x);
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
64 #else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
65 double sgngam;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
66
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
67 if (lo_ieee_isnan (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
68 result = x;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
69 else if (x <= 0 || lo_ieee_isinf (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
70 result = octave_Inf;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
71 else
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
72 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
73 #endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
74 return result;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
75 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
76
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
77 /* ---- pprsc.c from Stadloeber's winrand --- */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
78
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
79 /* flogfak(k) = ln(k!) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
80 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
81 flogfak (double k)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
82 {
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
83 #define C0 9.18938533204672742e-01
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
84 #define C1 8.33333333333333333e-02
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
85 #define C3 -2.77777777777777778e-03
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
86 #define C5 7.93650793650793651e-04
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
87 #define C7 -5.95238095238095238e-04
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
89 static double logfak[30L] =
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
90 {
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
91 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
92 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
93 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
94 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
95 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
96 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
97 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
98 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
99 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
100 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
101 };
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
102
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
103 double r, rr;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
104
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
105 if (k >= 30.0)
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 r = 1.0 / k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
108 rr = r * r;
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
109 return ((k + 0.5)*log (k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
5742
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 else
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112 return (logfak[(int)k]);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
113 }
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 /******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
117 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
118 * Poisson Distribution - Patchwork Rejection/Inversion *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
119 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
120 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
121 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 * For parameter my < 10 Tabulated Inversion is applied. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
123 * For my >= 10 Patchwork Rejection is employed: *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
124 * 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
125 * its body by certain point reflections. Within a large center *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
126 * interval variates are sampled efficiently by rejection from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
127 * uniform hats. Rectangular immediate acceptance regions speed *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
128 * up the generation. The remaining tails are covered by *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
129 * exponential functions. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
130 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
131 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
132 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
133 * FUNCTION : - pprsc samples a random number from the Poisson *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
134 * distribution with parameter my > 0. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
135 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
136 * continuous and discrete unimodal distributions, *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
137 * Doctoral Dissertation, 156 pp., Technical *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
138 * University Graz, Austria. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
140 * unsigned long integer *seed. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
141 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
142 * Implemented by H. Zechner, January 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
143 * Revised by F. Niederl, July 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
144 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
145 ******************************************************************/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
146
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
147 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
148 f (double k, double l_nu, double c_pm)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
149 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
150 return exp (k * l_nu - flogfak (k) - c_pm);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
151 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
152
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
153 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
154 pprsc (double my)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
155 {
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
156 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
157 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
158 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
159 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
160 double Dk, X, Y;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
161 double Ds, U, V, W;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
162
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
163 if (my != my_last)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
164 { /* set-up */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
165 my_last = my;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
166 /* 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
167 Ds = sqrt (my + 0.25);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
168
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
169 /* 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
170 /* which delimit the centre region of h(x) */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
171 m = floor (my);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
172 k2 = ceil (my - 0.5 - Ds);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
173 k4 = floor (my - 0.5 + Ds);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
174 k1 = k2 + k2 - m + 1L;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
175 k5 = k4 + k4 - m;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
176
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
177 /* range width of the critical left and right centre region */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
178 dl = (k2 - k1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
179 dr = (k5 - k4);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
180
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
181 /* 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
182 r1 = my / k1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
183 r2 = my / k2;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184 r4 = my / (k4 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185 r5 = my / (k5 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
187 /* reciprocal values of the scale parameters of exp. tail envelope */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
188 ll = log (r1); /* expon. tail left */
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
189 lr = -log (r5); /* expon. tail right*/
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
190
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
191 /* Poisson constants, necessary for computing function values f(k) */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
192 l_my = log (my);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
193 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
194
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
195 /* 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
196 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
197 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
198 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
199 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
200
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
201 /* 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
202 /* area of the two immediate acceptance regions between k2, k4 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
203 p1 = f2 * (dl + 1.0); /* immed. left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
204 p2 = f2 * dl + p1; /* centre left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
205 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
206 p4 = f4 * dr + p3; /* centre right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
207 p5 = f1 / ll + p4; /* exp. tail left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
208 p6 = f5 / lr + p5; /* exp. tail right*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
209 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
210
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
211 for (;;)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
212 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
213 /* generate uniform number U -- U(0, p6) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
214 /* case distinction corresponding to U */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
215 if ((U = RUNI * p6) < p2)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
216 { /* centre left */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
217
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
218 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
219 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
220 if ((V = U - p1) < 0.0) return (k2 + floor (U/f2));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
221 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
222 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
223 if ((W = V / dl) < f1 ) return (k1 + floor (V/f1));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
224
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
225 /* 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
226 /* either squeeze-acceptance of X or acceptance-rejection of Y */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
227 Dk = floor (dl * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
228 if (W <= f2 - Dk * (f2 - f2/r2))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
229 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
230 return (k2 - Dk); /* X = k2 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
231 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
232 if ((V = f2 + f2 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
233 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
234 Y = k2 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
235 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
236 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
237 return (Y); /* Y = k2 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
238 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
239 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
240 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
241 X = k2 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
242 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
243 else if (U < p4)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
244 { /* centre right */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
245 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
246 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
247 if ((V = U - p3) < 0.0) return (k4 - floor ((U - p2)/f4));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
248 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
249 R4 = [k4+1, k5+1)*[0, f5) */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
250 if ((W = V / dr) < f5 ) return (k5 - floor (V/f5));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
251
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
252 /* 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
253 /* either squeeze-acceptance of X or acceptance-rejection of Y */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
254 Dk = floor (dr * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
255 if (W <= f4 - Dk * (f4 - f4*r4))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
256 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
257 return (k4 + Dk); /* X = k4 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
258 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
259 if ((V = f4 + f4 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
260 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
261 Y = k4 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
262 if (V <= f4 + Dk * (1.0 - f4)/ dr)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
263 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
264 return (Y); /* Y = k4 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
265 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
266 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
267 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
268 X = k4 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
269 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
270 else
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
271 {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
272 W = RUNI;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
273 if (U < p5)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
274 { /* expon. tail left */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
275 Dk = floor (1.0 - log (W)/ll);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
276 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
277 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
278 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
279 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
280 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
281 else
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
282 { /* expon. tail right*/
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
283 Dk = floor (1.0 - log (W)/lr);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
284 X = k5 + Dk; /* X >= k5 + 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
285 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
286 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
287 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
288 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
289 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
290
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
291 /* acceptance-rejection test of candidate X from the original area */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
292 /* 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
293 /* log f(X) = (X - m)*log(my) - log X! + log m! */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
294 if (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
295 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
296 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
297 /* ---- pprsc.c end ------ */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
298
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 /* The remainder of the file is by Paul Kienzle */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
301
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
302 /* 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
303 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
304 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
305 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
306 /* Table size is predicated on the maximum value of lambda
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
307 * 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
308 * returned by the uniform random number generator on [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
309 * 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
310 * 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
311 * 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
312 * 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
313 * 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
314 * will need to be longer still. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
315 #define TABLESIZE 46
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
316 double t[TABLESIZE];
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
317
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
318 /* 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
319 * We will almost certainly need it. */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
320 int intlambda = (int)floor (lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
321 double P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
322 int tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
323 size_t i = n;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
324
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
325 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
326 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
327 {
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
328 P = P*lambda/(double)tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
329 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
330 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
331
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
332 while (i-- > 0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
333 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
334 double u = RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
335
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
336 /* 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
337 * 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
338 * 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
339 * 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
340 * problem. */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
341 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
342
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
343 /* 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
344 * 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
345 * 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
346 * the inner loop. */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
347 nextk:
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
348 if (u <= t[k])
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 p[i] = (double) k;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
351 continue;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
352 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
353 if (++k < tableidx) goto nextk;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
354
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
355 /* 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
356 * 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
357 while (tableidx < TABLESIZE)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
358 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
359 P = P*lambda/(double)tableidx;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
360 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
361 /* 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
362 * 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
363 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
364 tableidx++;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
365 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
366 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
367
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
368 /* 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
369 * 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
370 * the range [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
371 p[i] = (double)(tableidx-1);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
372 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
373 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
374
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
375 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
376 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
377 {
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 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
379
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
380 /* 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
381 * We will almost certainly need it. */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
382 int intlambda = (int)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
383 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
384 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
385 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
386
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
387 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
388 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
389 {
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
390 P = P*lambda/(double)tableidx;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
391 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
392 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
393
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
394 while (i-- > 0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
395 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
396 double u = RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
397 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
398 nextk:
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
399 if (u <= t[k])
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
400 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
401 p[i] = (float) k;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
402 continue;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
403 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
404 if (++k < tableidx) goto nextk;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
405
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
406 while (tableidx < TABLESIZE)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
407 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
408 P = P*lambda/(double)tableidx;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
409 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
410 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
411 tableidx++;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
412 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
413 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
414
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
415 p[i] = (float)(tableidx-1);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
416 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
417 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
418
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
419 /* From Press, et al., Numerical Recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
420 static void
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
421 poisson_rejection (double lambda, double *p, size_t n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
422 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
423 double sq = sqrt (2.0*lambda);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
424 double alxm = log (lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
425 double g = lambda*alxm - LGAMMA(lambda+1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
426 size_t i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
427
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
428 for (i = 0; i < n; i++)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
429 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
430 double y, em, t;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
431 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
432 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
433 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
434 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
435 y = tan (M_PI*RUNI);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
436 em = sq * y + lambda;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
437 } while (em < 0.0);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
438 em = floor (em);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
439 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
440 } while (RUNI > t);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
441 p[i] = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
442 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
443 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
444
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
445 /* 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
446 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
447 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
448 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
449 double sq = sqrt (2.0*lambda);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
450 double alxm = 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
451 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
452 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
453
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
454 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
455 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
456 double y, em, t;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
457 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
458 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
459 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
460 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
461 y = tan (M_PI*RUNI);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
462 em = sq * y + lambda;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
463 } while (em < 0.0);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
464 em = floor (em);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
465 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
466 } 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
467 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
468 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
469 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
470
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
471 /* 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
472 * the normal approximation is based on:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
473 * > L=1e8; x=floor(linspace(0,2*L,1000));
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
474 * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
475 * ans = 1.1376e-28
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
476 * 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
477 * 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
478 * 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
479
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
480 /* 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
481 void
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
482 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
483 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
484 octave_idx_type i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
485 if (L < 0.0 || INFINITE(L))
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
486 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
487 for (i=0; i<n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
488 p[i] = NAN;
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 <= 10.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
491 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
492 poisson_cdf_lookup (L, p, n);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
493 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
494 else if (L <= 1e8)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
495 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
496 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
497 p[i] = pprsc (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
498 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
499 else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
500 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
501 /* 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
502 const double sqrtL = sqrt (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
503 for (i = 0; i < n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
504 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
505 p[i] = floor (RNOR*sqrtL + L + 0.5);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
506 if (p[i] < 0.0)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
507 p[i] = 0.0; /* will probably never happen */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
508 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
509 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
510 }
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 /* Generate one poisson variate */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
513 double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
514 oct_randp (double L)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
515 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
516 double ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
517 if (L < 0.0) ret = NAN;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
518 else if (L <= 12.0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
519 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
520 /* 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
521 double g = exp (-L);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
522 int em = -1;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
523 double t = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
524 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
525 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
526 ++em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
527 t *= RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
528 } while (t > g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
529 ret = em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
530 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
531 else if (L <= 1e8)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
532 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
533 /* numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
534 poisson_rejection (L, &ret, 1);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
535 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
536 else if (INFINITE(L))
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 /* FIXME R uses NaN, but the normal approx. suggests that as
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
539 * limit should be inf. Which is correct? */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
540 ret = NAN;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
541 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
542 else
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
543 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
544 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
545 ret = floor (RNOR*sqrt (L) + L + 0.5);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
546 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
547 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
548 return ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
549 }
14655
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 /* 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
552 void
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
553 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
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 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
556 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
557 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
558 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
559 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
560 p[i] = NAN;
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 <= 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
563 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
564 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
565 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
566 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
567 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
568 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
569 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
570 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
571 else
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
572 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
573 /* 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
574 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
575 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
576 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
577 p[i] = 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
578 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
579 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
580 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
581 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
582 }
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 /* 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
585 float
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
586 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
587 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
588 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
589 float ret;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
590 if (L < 0.0) ret = NAN;
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
591 else if (L <= 12.0)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
592 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
593 /* 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
594 double g = exp (-L);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
595 int em = -1;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
596 double t = 1.0;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
597 do
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
598 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
599 ++em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
600 t *= RUNI;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
601 } while (t > g);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
602 ret = em;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
603 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
604 else if (L <= 1e8)
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
605 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
606 /* numerical recipes */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
607 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
608 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
609 else if (INFINITE(L))
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 /* FIXME R uses NaN, but the normal approx. suggests that as
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
612 * limit should be inf. Which is correct? */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
613 ret = NAN;
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
614 }
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
615 else
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
616 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
617 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
618 ret = floor (RNOR*sqrt (L) + L + 0.5);
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
619 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
620 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
621 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
622 }