octave: liboctave/numeric/randpoisson.cc annotate

annotate liboctave/numeric/randpoisson.cc @ 30564:796f54d4ddbf stable

update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2021. In all .txi and .texi files except gpl.txi and gpl.texi in the doc/liboctave and doc/interpreter directories, change the copyright to "Octave Project Developers", the same as used for other source files. Update copyright notices for 2022 (not done since 2019). For gpl.txi and gpl.texi, change the copyright notice to be "Free Software Foundation, Inc." and leave the date at 2007 only because this file only contains the text of the GPL, not anything created by the Octave Project Developers. Add Paul Thomas to contributors.in.

author	John W. Eaton <jwe@octave.org>
date	Tue, 28 Dec 2021 18:22:40 -0500
parents	32f4357ac8d9
children	e88a07dec498

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ////////////////////////////////////////////////////////////////////////
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 //
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 29655 diff changeset	3 // Copyright (C) 2006-2022 The Octave Project Developers
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	4 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	5 // See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	6 // distribution or <https://octave.org/copyright/>.
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	7 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	8 // This file is part of Octave.
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	9 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	10 // Octave is free software: you can redistribute it and/or modify it
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	11 // under the terms of the GNU General Public License as published by
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	12 // the Free Software Foundation, either version 3 of the License, or
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	13 // (at your option) any later version.
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	14 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	15 // Octave is distributed in the hope that it will be useful, but
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	16 // WITHOUT ANY WARRANTY; without even the implied warranty of
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	18 // GNU General Public License for more details.
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	19 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	20 // You should have received a copy of the GNU General Public License
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	21 // along with Octave; see the file COPYING. If not, see
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	22 // <https://www.gnu.org/licenses/>.
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	23 //
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	24 ////////////////////////////////////////////////////////////////////////
7019 4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 6198 diff changeset	25
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 6198 diff changeset	26 /* Original version written by Paul Kienzle distributed as free
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 6198 diff changeset	27 software in the in the public domain. */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	28
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	29 #if defined (HAVE_CONFIG_H)
21301 40de9f8f23a6 Use '#include "config.h"' rather than <config.h>. Rik <rik@octave.org> parents: 21202 diff changeset	30 # include "config.h"
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	31 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	32
23662 bd77ab816e43 eliminate obsolete file lo-math.h John W. Eaton <jwe@octave.org> parents: 23649 diff changeset	33 #include <cmath>
23475 d691ed308237 maint: Clean up #includes in liboctave/numeric directory. Rik <rik@octave.org> parents: 23220 diff changeset	34 #include <cstddef>
d691ed308237 maint: Clean up #includes in liboctave/numeric directory. Rik <rik@octave.org> parents: 23220 diff changeset	35
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	36 #include "f77-fcn.h"
7231 2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe] jwe parents: 7019 diff changeset	37 #include "lo-error.h"
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	38 #include "lo-ieee.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	39 #include "randmtzig.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	40 #include "randpoisson.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	41
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	42 namespace octave
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	43 {
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	44 static double xlgamma (double x)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	45 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	46 return std::lgamma (x);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	47 }
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	48
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	49 /* ---- pprsc.c from Stadloeber's winrand --- */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	50
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	51 /* flogfak(k) = ln(k!) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	52 static double flogfak (double k)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	53 {
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	54 #define C0 9.18938533204672742e-01
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	55 #define C1 8.33333333333333333e-02
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	56 #define C3 -2.77777777777777778e-03
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	57 #define C5 7.93650793650793651e-04
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	58 #define C7 -5.95238095238095238e-04
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	59
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	60 static double logfak[30L] =
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	61 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	62 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	63 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	64 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	65 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	66 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	67 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	68 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	69 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	70 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	71 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	72 };
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	73
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	74 double r, rr;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	75
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	76 if (k >= 30.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	77 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	78 r = 1.0 / k;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	79 rr = r * r;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	80 return ((k + 0.5)*std::log (k) - k + C0
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	81 + r(C1 + rr(C3 + rr(C5 + rrC7))));
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	82 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	83 else
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	84 return (logfak[static_cast<int> (k)]);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	85 }
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	86
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	87 /******************************************************************
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	88 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	89 * Poisson Distribution - Patchwork Rejection/Inversion *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	90 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	91 ******************************************************************
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	92 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	93 * For parameter my < 10, Tabulated Inversion is applied. *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	94 * For my >= 10, Patchwork Rejection is employed: *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	95 * The area below the histogram function f(x) is rearranged in *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	96 * its body by certain point reflections. Within a large center *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	97 * interval variates are sampled efficiently by rejection from *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	98 * uniform hats. Rectangular immediate acceptance regions speed *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	99 * up the generation. The remaining tails are covered by *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	100 * exponential functions. *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	101 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	102 ******************************************************************
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	103 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	104 * FUNCTION : - pprsc samples a random number from the Poisson *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	105 * distribution with parameter my > 0. *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	106 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	107 * continuous and discrete unimodal distributions, *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	108 * Doctoral Dissertation, 156 pp., Technical *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	109 * University Graz, Austria. *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	110 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	111 * unsigned long integer seed.
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	112 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	113 * Implemented by H. Zechner, January 1994 *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	114 * Revised by F. Niederl, July 1994 *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	115 * *
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	116 ******************************************************************/
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	117
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	118 static double f (double k, double l_nu, double c_pm)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	119 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	120 return exp (k * l_nu - flogfak (k) - c_pm);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	121 }
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	122
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	123 static double pprsc (double my)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	124 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	125 static double my_last = -1.0;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	126 static double m, k2, k4, k1, k5;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	127 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	128 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	129 double Dk, X, Y;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	130 double Ds, U, V, W;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	131
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	132 if (my != my_last)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	133 { /* set-up */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	134 my_last = my;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	135 /* approximate deviation of reflection points k2, k4 from my - 1/2 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	136 Ds = std::sqrt (my + 0.25);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	137
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	138 /* mode m, reflection points k2 and k4, and points k1 and k5, */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	139 /* which delimit the centre region of h(x) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	140 m = std::floor (my);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	141 k2 = ceil (my - 0.5 - Ds);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	142 k4 = std::floor (my - 0.5 + Ds);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	143 k1 = k2 + k2 - m + 1L;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	144 k5 = k4 + k4 - m;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	145
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	146 /* range width of the critical left and right centre region */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	147 dl = (k2 - k1);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	148 dr = (k5 - k4);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	149
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	150 /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	151 r1 = my / k1;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	152 r2 = my / k2;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	153 r4 = my / (k4 + 1.0);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	154 r5 = my / (k5 + 1.0);
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	155
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	156 /* reciprocal values of the scale parameters of exp. tail envelope */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	157 ll = std::log (r1); /* expon. tail left */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	158 lr = -std::log (r5); /* expon. tail right*/
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	159
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	160 /* Poisson constants, necessary for computing function values f(k) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	161 l_my = std::log (my);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	162 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	163
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	164 /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	165 f2 = f (k2, l_my, c_pm);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	166 f4 = f (k4, l_my, c_pm);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	167 f1 = f (k1, l_my, c_pm);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	168 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	169
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	170 /* area of the two centre and the two exponential tail regions */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	171 /* area of the two immediate acceptance regions between k2, k4 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	172 p1 = f2 * (dl + 1.0); /* immed. left */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	173 p2 = f2 * dl + p1; /* centre left */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	174 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	175 p4 = f4 * dr + p3; /* centre right */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	176 p5 = f1 / ll + p4; /* exp. tail left */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	177 p6 = f5 / lr + p5; /* exp. tail right*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	178 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	179
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	180 for (;;)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	181 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	182 /* generate uniform number U -- U(0, p6) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	183 /* case distinction corresponding to U */
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	184 if ((U = rand_uniform<double> () * p6) < p2)
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	185 { /* centre left */
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	186
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	187 /* immediate acceptance region
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	188 R2 = [k2, m) [0, f2), X = k2, ... m -1 /
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	189 if ((V = U - p1) < 0.0) return (k2 + std::floor (U/f2));
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	190 /* immediate acceptance region
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	191 R1 = [k1, k2)[0, f1), X = k1, ... k2-1 /
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	192 if ((W = V / dl) < f1 ) return (k1 + std::floor (V/f1));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	193
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	194 /* computation of candidate X < k2, and its counterpart Y > k2 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	195 /* either squeeze-acceptance of X or acceptance-rejection of Y */
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	196 Dk = std::floor (dl * rand_uniform<double> ()) + 1.0;
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	197 if (W <= f2 - Dk * (f2 - f2/r2))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	198 { /* quick accept of */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	199 return (k2 - Dk); /* X = k2 - Dk */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	200 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	201 if ((V = f2 + f2 - W) < 1.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	202 { /* quick reject of Y*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	203 Y = k2 + Dk;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	204 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	205 { /* quick accept of */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	206 return (Y); /* Y = k2 + Dk */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	207 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	208 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	209 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	210 X = k2 - Dk;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	211 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	212 else if (U < p4)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	213 { /* centre right */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	214 /* immediate acceptance region
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	215 R3 = [m, k4+1)[0, f4), X = m, ... k4 /
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	216 if ((V = U - p3) < 0.0) return (k4 - std::floor ((U - p2)/f4));
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	217 /* immediate acceptance region
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	218 R4 = [k4+1, k5+1)[0, f5) /
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	219 if ((W = V / dr) < f5 ) return (k5 - std::floor (V/f5));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	220
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	221 /* computation of candidate X > k4, and its counterpart Y < k4 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	222 /* either squeeze-acceptance of X or acceptance-rejection of Y */
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	223 Dk = std::floor (dr * rand_uniform<double> ()) + 1.0;
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	224 if (W <= f4 - Dk * (f4 - f4*r4))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	225 { /* quick accept of */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	226 return (k4 + Dk); /* X = k4 + Dk */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	227 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	228 if ((V = f4 + f4 - W) < 1.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	229 { /* quick reject of Y*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	230 Y = k4 - Dk;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	231 if (V <= f4 + Dk * (1.0 - f4)/ dr)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	232 { /* quick accept of */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	233 return (Y); /* Y = k4 - Dk */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	234 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	235 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	236 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	237 X = k4 + Dk;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	238 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	239 else
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	240 {
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	241 W = rand_uniform<double> ();
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	242 if (U < p5)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	243 { /* expon. tail left */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	244 Dk = std::floor (1.0 - std::log (W)/ll);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	245 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	246 W = (U - p4) ll; /* W -- U(0, h(x)) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	247 if (W <= f1 - Dk * (f1 - f1/r1))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	248 return (X); /* quick accept of X*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	249 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	250 else
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	251 { /* expon. tail right*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	252 Dk = std::floor (1.0 - std::log (W)/lr);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	253 X = k5 + Dk; /* X >= k5 + 1 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	254 W = (U - p5) lr; /* W -- U(0, h(x)) */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	255 if (W <= f5 - Dk * (f5 - f5*r5))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	256 return (X); /* quick accept of X*/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	257 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	258 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	259
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	260 /* acceptance-rejection test of candidate X from the original area */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	261 /* test, whether W <= f(k), with W = Uh(x) and U -- U(0, 1)/
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	262 /* log f(X) = (X - m)log(my) - log X! + log m! /
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	263 if (std::log (W) <= X * l_my - flogfak (X) - c_pm) return (X);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	264 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	265 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	266 /* ---- pprsc.c end ------ */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	267
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	268 /* The remainder of the file is by Paul Kienzle */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	269
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	270 /* Table size is predicated on the maximum value of lambda
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	271 * we want to store in the table, and the maximum value of
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	272 * returned by the uniform random number generator on [0,1).
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	273 * With lambda==10 and u_max = 1 - 1/(2^32+1), we
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	274 * have poisson_pdf(lambda,36) < 1-u_max. If instead our
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	275 * generator uses more bits of mantissa or returns a value
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	276 * in the range [0,1], then for lambda==10 we need a table
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	277 * size of 46 instead. For long doubles, the table size
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	278 * will need to be longer still. */
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	279 #define TABLESIZE 46
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	280
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	281 /* Given uniform u, find x such that CDF(L,x)==u. Return x. */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	282
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	283 template <typename T>
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	284 static void
29654 d13d090cb03a use std::size_t and std::ptrdiff_t in C++ code (bug #60471) John W. Eaton <jwe@octave.org> parents: 29358 diff changeset	285 poisson_cdf_lookup (double lambda, T *p, std::size_t n)
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	286 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	287 double t[TABLESIZE];
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	288
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	289 /* Precompute the table for the u up to and including 0.458.
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	290 * We will almost certainly need it. */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	291 int intlambda = static_cast<int> (std::floor (lambda));
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	292 double P;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	293 int tableidx;
29654 d13d090cb03a use std::size_t and std::ptrdiff_t in C++ code (bug #60471) John W. Eaton <jwe@octave.org> parents: 29358 diff changeset	294 std::size_t i = n;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	295
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	296 t[0] = P = exp (-lambda);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	297 for (tableidx = 1; tableidx <= intlambda; tableidx++)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	298 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	299 P = P*lambda/static_cast<double> (tableidx);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	300 t[tableidx] = t[tableidx-1] + P;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	301 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	302
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	303 while (i-- > 0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	304 {
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	305 double u = rand_uniform<double> ();
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	306
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	307 /* If u > 0.458 we know we can jump to floor(lambda) before
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	308 * comparing (this observation is based on Stadlober's winrand
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	309 * code). For lambda >= 1, this will be a win. Lambda < 1
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	310 * is already fast, so adding an extra comparison is not a
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	311 * problem. */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	312 int k = (u > 0.458 ? intlambda : 0);
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	313
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	314 /* We aren't using a for loop here because when we find the
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	315 * right k we want to jump to the next iteration of the
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	316 * outer loop, and the continue statement will only work for
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	317 * the inner loop. */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	318 nextk:
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	319 if (u <= t[k])
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	320 {
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	321 p[i] = static_cast<T> (k);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	322 continue;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	323 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	324 if (++k < tableidx)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	325 goto nextk;
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	326
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	327 /* We only need high values of the table very rarely so we
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	328 * don't automatically compute the entire table. */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	329 while (tableidx < TABLESIZE)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	330 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	331 P = P*lambda/static_cast<double> (tableidx);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	332 t[tableidx] = t[tableidx-1] + P;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	333 /* Make sure we converge to 1.0 just in case u is uniform
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	334 * on [0,1] rather than [0,1). */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	335 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	336 tableidx++;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	337 if (u <= t[tableidx-1]) break;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	338 }
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	339
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	340 /* We are assuming that the table size is big enough here.
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	341 * This should be true even if rand_uniform is returning values in
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	342 * the range [0,1] rather than [0,1). */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	343 p[i] = static_cast<T> (tableidx-1);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	344 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	345 }
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	346
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	347 /* From Press, et al., Numerical Recipes */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	348 template <typename T>
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	349 static void
29654 d13d090cb03a use std::size_t and std::ptrdiff_t in C++ code (bug #60471) John W. Eaton <jwe@octave.org> parents: 29358 diff changeset	350 poisson_rejection (double lambda, T *p, std::size_t n)
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	351 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	352 double sq = std::sqrt (2.0*lambda);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	353 double alxm = std::log (lambda);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	354 double g = lambda*alxm - xlgamma (lambda+1.0);
29654 d13d090cb03a use std::size_t and std::ptrdiff_t in C++ code (bug #60471) John W. Eaton <jwe@octave.org> parents: 29358 diff changeset	355 std::size_t i;
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	356
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	357 for (i = 0; i < n; i++)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	358 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	359 double y, em, t;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	360 do
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	361 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	362 do
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	363 {
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	364 y = tan (M_PI*rand_uniform<double> ());
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	365 em = sq * y + lambda;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	366 } while (em < 0.0);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	367 em = std::floor (em);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	368 t = 0.9(1.0+yy)exp (emalxm-flogfak (em)-g);
27102 84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	369 } while (rand_uniform<double> () > t);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	370 p[i] = em;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	371 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	372 }
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
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	374 /* The cutoff of L <= 1e8 in the following two functions before using
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	375 * the normal approximation is based on:
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	376 * > L=1e8; x=floor(linspace(0,2*L,1000));
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	377 * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	378 * ans = 1.1376e-28
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	379 * For L=1e7, the max is around 1e-9, which is within the step size of
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	380 * rand_uniform. For L>1e10 the pprsc function breaks down, as I saw
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	381 * from the histogram of a large sample, so 1e8 is both small enough
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	382 * and large enough. */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	383
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	384 /* Generate a set of poisson numbers with the same distribution */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	385 template <typename T> void rand_poisson (T L_arg, octave_idx_type n, T *p)
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	386 {
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	387 double L = L_arg;
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	388 octave_idx_type i;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	389 if (L < 0.0 \|\| lo_ieee_isinf (L))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	390 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	391 for (i=0; i<n; i++)
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	392 p[i] = numeric_limits<T>::NaN ();
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	393 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	394 else if (L <= 10.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	395 {
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	396 poisson_cdf_lookup<T> (L, p, n);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	397 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	398 else if (L <= 1e8)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	399 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	400 for (i=0; i<n; i++)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	401 p[i] = pprsc (L);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	402 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	403 else
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	404 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	405 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	406 const double sqrtL = std::sqrt (L);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	407 for (i = 0; i < n; i++)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	408 {
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	409 p[i] = std::floor (rand_normal<T> () * sqrtL + L + 0.5);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	410 if (p[i] < 0.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	411 p[i] = 0.0; /* will probably never happen */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	412 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	413 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	414 }
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	415
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	416 template void rand_poisson<double> (double, octave_idx_type, double *);
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	417 template void rand_poisson<float> (float, octave_idx_type, float *);
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	418
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	419 /* Generate one poisson variate */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	420 template <typename T> T rand_poisson (T L_arg)
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	421 {
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	422 double L = L_arg;
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	423 T ret;
859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	424 if (L < 0.0) ret = numeric_limits<T>::NaN ();
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	425 else if (L <= 12.0)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	426 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	427 /* From Press, et al. Numerical recipes */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	428 double g = exp (-L);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	429 int em = -1;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	430 double t = 1.0;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	431 do
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	432 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	433 ++em;
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	434 t *= rand_uniform<T> ();
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	435 } while (t > g);
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	436 ret = em;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	437 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	438 else if (L <= 1e8)
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	439 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	440 /* numerical recipes */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	441 poisson_rejection<T> (L, &ret, 1);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	442 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	443 else if (lo_ieee_isinf (L))
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	444 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	445 /* FIXME: R uses NaN, but the normal approximation suggests that
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	446 * limit should be Inf. Which is correct? */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	447 ret = numeric_limits<T>::NaN ();
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	448 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	449 else
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	450 {
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	451 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
25434 859ad1f0b85e use templates to eliminate duplicate code in random number generators John W. Eaton <jwe@octave.org> parents: 25433 diff changeset	452 ret = std::floor (rand_normal<T> () * std::sqrt (L) + L + 0.5);
25433 49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	453 if (ret < 0.0) ret = 0.0; /* will probably never happen */
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	454 }
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	455 return ret;
49e0447413ad use templates and move rand functions inside octave namespace John W. Eaton <jwe@octave.org> parents: 25054 diff changeset	456 }
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	457
29228 5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919). Markus Mützel <markus.muetzel@gmx.de> parents: 27923 diff changeset	458 template OCTAVE_API double rand_poisson<double> (double);
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919). Markus Mützel <markus.muetzel@gmx.de> parents: 27923 diff changeset	459 template OCTAVE_API float rand_poisson<float> (float);
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	460 }

Mercurial > octave

annotate liboctave/numeric/randpoisson.cc @ 30564:796f54d4ddbf stable