Mercurial > hg > octave-nkf

annotate liboctave/numeric/randpoisson.c @ 16057:c3057d80cf91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Created common octave_dock_widget class * octave_dock_widget.h : New class octave_dock_widget with common dock widget methods, slots and signals * file_editor_interface.h: Now inherit from octave_dock_widget, removed common dock widget code * file_editor.h: Now inherit from octave_dock_widget, removed common dock widget code * file_editor.cc: Now inherit from octave_dock_widget, removed common dock widget code * files_dockwidget.h: Now inherit from octave_dock_widget, removed common dock widget code * files_dockwidget.cc: Now inherit from octave_dock_widget, removed common dock widget code * history_dockwidget.h: Now inherit from octave_dock_widget, removed common dock widget code * history_dockwidget.cc: Now inherit from octave_dock_widget, removed common dock widget methods * terminal_dockwidget.h: Now inherit from octave_dock_widget, removed common dock widget code * terminal_dockwidget.cc: Now inherit from octave_dock_widget, removed common dock widget code * module.mk: Added octave_dock_widget
author Richard Crozier <richard.crozier@yahoo.co.uk>
date Wed, 06 Feb 2013 21:45:04 +0000
parents 4bfe605f5ecf
children d63878346099
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
14138
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 11586
diff changeset
3 Copyright (C) 2006-2012 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 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
83 #define C0 9.18938533204672742e-01
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
84 #define C1 8.33333333333333333e-02
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
85 #define C3 -2.77777777777777778e-03
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
86 #define C5 7.93650793650793651e-04
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
87 #define C7 -5.95238095238095238e-04
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
89 static double logfak[30L] = {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
90 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
91 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
92 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
93 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
94 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
95 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
96 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
97 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
98 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
99 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
100 };
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
101
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
102 double r, rr;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
103
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
104 if (k >= 30.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
105 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
106 r = 1.0 / k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
107 rr = r * r;
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
108 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
109 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
110 else
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
111 return (logfak[(int)k]);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112 }
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 * Poisson Distribution - Patchwork Rejection/Inversion *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
118 * *
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 * For parameter my < 10 Tabulated Inversion is applied. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 * For my >= 10 Patchwork Rejection is employed: *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
123 * 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
124 * its body by certain point reflections. Within a large center *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
125 * interval variates are sampled efficiently by rejection from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
126 * uniform hats. Rectangular immediate acceptance regions speed *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
127 * up the generation. The remaining tails are covered by *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
128 * exponential functions. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
129 * *
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 * FUNCTION : - pprsc samples a random number from the Poisson *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
133 * distribution with parameter my > 0. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
134 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
135 * continuous and discrete unimodal distributions, *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
136 * Doctoral Dissertation, 156 pp., Technical *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
137 * University Graz, Austria. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
138 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139 * unsigned long integer *seed. *
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 * Implemented by H. Zechner, January 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
142 * Revised by F. Niederl, July 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
143 * *
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
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
146 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
147 f (double k, double l_nu, double c_pm)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
148 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
149 return exp (k * l_nu - flogfak (k) - c_pm);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
150 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
151
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
152 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
153 pprsc (double my)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
154 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
155 static double my_last = -1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
156 static double m, k2, k4, k1, k5;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
157 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
158 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
159 double Dk, X, Y;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
160 double Ds, U, V, W;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
161
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
162 if (my != my_last)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
163 { /* set-up */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
164 my_last = my;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
165 /* 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
166 Ds = sqrt (my + 0.25);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
167
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
168 /* 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
169 /* 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
170 m = floor (my);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
171 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
172 k4 = floor (my - 0.5 + Ds);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
173 k1 = k2 + k2 - m + 1L;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
174 k5 = k4 + k4 - m;
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 /* range width of the critical left and right centre region */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
177 dl = (k2 - k1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
178 dr = (k5 - k4);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
179
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
180 /* 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
181 r1 = my / k1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
182 r2 = my / k2;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
183 r4 = my / (k4 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184 r5 = my / (k5 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 /* 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
187 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
188 lr = -log (r5); /* expon. tail right*/
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 /* 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
191 l_my = log (my);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
192 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
193
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
194 /* 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
195 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
196 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
197 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
198 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
199
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
200 /* 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
201 /* area of the two immediate acceptance regions between k2, k4 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
202 p1 = f2 * (dl + 1.0); /* immed. left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
203 p2 = f2 * dl + p1; /* centre left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
204 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
205 p4 = f4 * dr + p3; /* centre right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
206 p5 = f1 / ll + p4; /* exp. tail left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
207 p6 = f5 / lr + p5; /* exp. tail right*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
208 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
209
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
210 for (;;)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
211 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
212 /* generate uniform number U -- U(0, p6) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
213 /* case distinction corresponding to U */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
214 if ((U = RUNI * p6) < p2)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
215 { /* centre left */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
216
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
217 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
218 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
219 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
220 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
221 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
222 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
223
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
224 /* 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
225 /* 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
226 Dk = floor (dl * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
227 if (W <= f2 - Dk * (f2 - f2/r2))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
228 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
229 return (k2 - Dk); /* X = k2 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
230 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
231 if ((V = f2 + f2 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
232 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
233 Y = k2 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
234 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
235 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
236 return (Y); /* Y = k2 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
237 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
238 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
239 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
240 X = k2 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
241 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
242 else if (U < p4)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
243 { /* centre right */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
244 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
245 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
246 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
247 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
248 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
249 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
250
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
251 /* 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
252 /* 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
253 Dk = floor (dr * RUNI) + 1.0;
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
254 if (W <= f4 - Dk * (f4 - f4*r4))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
255 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
256 return (k4 + Dk); /* X = k4 + Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
257 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
258 if ((V = f4 + f4 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
259 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
260 Y = k4 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
261 if (V <= f4 + Dk * (1.0 - f4)/ dr)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
262 { /* quick accept of */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
263 return (Y); /* Y = k4 - Dk */
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
264 }
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
265 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
266 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
267 X = k4 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
268 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
269 else
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
270 {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
271 W = RUNI;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
272 if (U < p5)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
273 { /* expon. tail left */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
274 Dk = floor (1.0 - log (W)/ll);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
275 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
276 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
277 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
278 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
279 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
280 else
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
281 { /* expon. tail right*/
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
282 Dk = floor (1.0 - log (W)/lr);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
283 X = k5 + Dk; /* X >= k5 + 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
284 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
285 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
286 return (X); /* quick accept of X*/
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
287 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
288 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
289
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
290 /* acceptance-rejection test of candidate X from the original area */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
291 /* 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
292 /* 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
293 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
294 }
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 /* ---- pprsc.c end ------ */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
297
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 /* The remainder of the file is by Paul Kienzle */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
300
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
301 /* 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
302 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
303 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
304 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
305 /* Table size is predicated on the maximum value of lambda
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
306 * 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
307 * returned by the uniform random number generator on [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
308 * 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
309 * 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
310 * 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
311 * 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
312 * 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
313 * will need to be longer still. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
314 #define TABLESIZE 46
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
315 double t[TABLESIZE];
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
316
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
317 /* 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
318 * 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
319 int intlambda = (int)floor (lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
320 double P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
321 int tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
322 size_t i = n;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
323
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
324 t[0] = P = exp (-lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
325 for (tableidx = 1; tableidx <= intlambda; tableidx++) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
326 P = P*lambda/(double)tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
327 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
328 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
329
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
330 while (i-- > 0) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
331 double u = RUNI;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
332
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
333 /* If u > 0.458 we know we can jump to floor(lambda) before
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
334 * comparing (this observation is based on Stadlober's winrand
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
335 * code). For lambda >= 1, this will be a win. Lambda < 1
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
336 * is already fast, so adding an extra comparison is not a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
337 * problem. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
338 int k = (u > 0.458 ? intlambda : 0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
339
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
340 /* We aren't using a for loop here because when we find the
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
341 * right k we want to jump to the next iteration of the
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
342 * outer loop, and the continue statement will only work for
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
343 * the inner loop. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
344 nextk:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
345 if ( u <= t[k] ) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
346 p[i] = (double) k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
347 continue;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
348 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
349 if (++k < tableidx) goto nextk;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
350
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
351 /* We only need high values of the table very rarely so we
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
352 * don't automatically compute the entire table. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
353 while (tableidx < TABLESIZE) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
354 P = P*lambda/(double)tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
355 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
356 /* Make sure we converge to 1.0 just in case u is uniform
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
357 * on [0,1] rather than [0,1). */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
358 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
359 tableidx++;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
360 if (u <= t[tableidx-1]) break;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
361 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
362
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
363 /* We are assuming that the table size is big enough here.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
364 * This should be true even if RUNI is returning values in
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
365 * the range [0,1] rather than [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
366 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
367 p[i] = (double)(tableidx-1);
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 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
370
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
371 static void
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
372 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
373 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
374 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
375
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
376 /* 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
377 * 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
378 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
379 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
380 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
381 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
382
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
383 t[0] = P = exp (-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
384 for (tableidx = 1; tableidx <= intlambda; 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 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
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
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
389 while (i-- > 0) {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
390 double u = RUNI;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
391 int k = (u > 0.458 ? intlambda : 0);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
392 nextk:
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
393 if ( u <= t[k] ) {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
394 p[i] = (float) k;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
395 continue;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
396 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
397 if (++k < tableidx) goto nextk;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
398
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
399 while (tableidx < TABLESIZE) {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
400 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
401 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
402 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 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
403 tableidx++;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
404 if (u <= t[tableidx-1]) break;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
405 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
406
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
407 p[i] = (float)(tableidx-1);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
408 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
409 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
410
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
411 /* From Press, et al., Numerical Recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
412 static void
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
413 poisson_rejection (double lambda, double *p, size_t n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
414 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
415 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
416 double alxm = log (lambda);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
417 double g = lambda*alxm - LGAMMA(lambda+1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
418 size_t i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
419
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
420 for (i = 0; i < n; i++)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
421 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
422 double y, em, t;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
423 do {
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
424 do {
15851
4bfe605f5ecf randpoisson.c: Use Octave spacing convention on function calls.
Rik <rik@octave.org>
parents: 15271
diff changeset
425 y = tan (M_PI*RUNI);
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
426 em = sq * y + lambda;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
427 } while (em < 0.0);
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
428 em = floor (em);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
429 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
430 } while (RUNI > t);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
431 p[i] = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
432 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
433 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
434
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
435 /* 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
436 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
437 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
438 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
439 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
440 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
441 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
442 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
443
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
444 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
445 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
446 double y, em, t;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
447 do {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
448 do {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
449 y = tan (M_PI*RUNI);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
450 em = sq * y + lambda;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
451 } while (em < 0.0);
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
452 em = floor (em);
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
453 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
454 } while (RUNI > t);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
455 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
456 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
457 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
458
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
459 /* 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
460 * the normal approximation is based on:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
461 * > L=1e8; x=floor(linspace(0,2*L,1000));
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
462 * > 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
463 * ans = 1.1376e-28
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
464 * 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
465 * 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
466 * 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
467
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
468 /* 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
469 void
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
470 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
471 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
472 octave_idx_type i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
473 if (L < 0.0 || INFINITE(L))
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
474 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
475 for (i=0; i<n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
476 p[i] = NAN;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
477 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
478 else if (L <= 10.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
479 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
480 poisson_cdf_lookup (L, p, n);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
481 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
482 else if (L <= 1e8)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
483 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
484 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
485 p[i] = pprsc (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
486 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
487 else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
488 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
489 /* 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
490 const double sqrtL = sqrt (L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
491 for (i = 0; i < n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
492 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
493 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
494 if (p[i] < 0.0)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
495 p[i] = 0.0; /* will probably never happen */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
496 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
497 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
498 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
499
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
500 /* Generate one poisson variate */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
501 double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
502 oct_randp (double L)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
503 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
504 double ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
505 if (L < 0.0) ret = NAN;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
506 else if (L <= 12.0) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
507 /* From Press, et al. Numerical recipes */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
508 double g = exp (-L);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
509 int em = -1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
510 double t = 1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
511 do {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
512 ++em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
513 t *= RUNI;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
514 } while (t > g);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
515 ret = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
516 } else if (L <= 1e8) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
517 /* numerical recipes */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
518 poisson_rejection (L, &ret, 1);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
519 } else if (INFINITE(L)) {
5775
ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe]
jwe
parents: 5742
diff changeset
520 /* FIXME R uses NaN, but the normal approx. suggests that as
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
521 * limit should be inf. Which is correct? */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
522 ret = NAN;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
523 } else {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
524 /* 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
525 ret = floor (RNOR*sqrt (L) + L + 0.5);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
526 if (ret < 0.0) ret = 0.0; /* will probably never happen */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
527 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
528 return ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
529 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
530
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
531 /* 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
532 void
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
533 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
534 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
535 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
536 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
537 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
538 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
539 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
540 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
541 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
542 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
543 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
544 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
545 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
546 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
547 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
548 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
549 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
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 else
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
552 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
553 /* 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
554 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
555 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
556 {
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
557 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
558 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
559 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
560 }
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 }
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 /* 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
565 float
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
566 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
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 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
569 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
570 if (L < 0.0) ret = NAN;
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 if (L <= 12.0) {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
572 /* From Press, et al. Numerical recipes */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
573 double g = exp (-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
574 int em = -1;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
575 double t = 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
576 do {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
577 ++em;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
578 t *= RUNI;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
579 } while (t > g);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
580 ret = em;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
581 } 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
582 /* numerical recipes */
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14655
diff changeset
583 poisson_rejection_float (L, &ret, 1);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
584 } else if (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
585 /* FIXME R uses NaN, but the normal approx. suggests that as
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
586 * limit should be inf. Which is correct? */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
587 ret = NAN;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
588 } else {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
589 /* 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
590 ret = floor (RNOR*sqrt (L) + 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
591 if (ret < 0.0) ret = 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
592 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
593 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
594 }