octave-nkf: liboctave/numeric/randgamma.c annotate

annotate liboctave/numeric/randgamma.c @ 17769:49a5a4be04a1

maint: Use GNU style coding conventions for code in liboctave/ * liboctave/array/Array-C.cc, liboctave/array/Array-b.cc, liboctave/array/Array-ch.cc, liboctave/array/Array-d.cc, liboctave/array/Array-f.cc, liboctave/array/Array-fC.cc, liboctave/array/Array-util.cc, liboctave/array/Array-util.h, liboctave/array/Array.cc, liboctave/array/Array.h, liboctave/array/Array3.h, liboctave/array/CColVector.cc, liboctave/array/CColVector.h, liboctave/array/CDiagMatrix.cc, liboctave/array/CDiagMatrix.h, liboctave/array/CMatrix.cc, liboctave/array/CMatrix.h, liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h, liboctave/array/CRowVector.cc, liboctave/array/CRowVector.h, liboctave/array/CSparse.cc, liboctave/array/CSparse.h, liboctave/array/DiagArray2.h, liboctave/array/MArray.cc, liboctave/array/MArray.h, liboctave/array/MDiagArray2.cc, liboctave/array/MDiagArray2.h, liboctave/array/MSparse.cc, liboctave/array/MSparse.h, liboctave/array/MatrixType.cc, liboctave/array/MatrixType.h, liboctave/array/PermMatrix.h, liboctave/array/Range.cc, liboctave/array/Range.h, liboctave/array/Sparse.cc, liboctave/array/Sparse.h, liboctave/array/boolMatrix.cc, liboctave/array/boolMatrix.h, liboctave/array/boolNDArray.cc, liboctave/array/boolNDArray.h, liboctave/array/boolSparse.cc, liboctave/array/boolSparse.h, liboctave/array/chMatrix.cc, liboctave/array/chMatrix.h, liboctave/array/chNDArray.cc, liboctave/array/chNDArray.h, liboctave/array/dColVector.h, liboctave/array/dDiagMatrix.cc, liboctave/array/dDiagMatrix.h, liboctave/array/dMatrix.cc, liboctave/array/dMatrix.h, liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h, liboctave/array/dRowVector.h, liboctave/array/dSparse.cc, liboctave/array/dSparse.h, liboctave/array/dim-vector.cc, liboctave/array/dim-vector.h, liboctave/array/fCColVector.cc, liboctave/array/fCColVector.h, liboctave/array/fCDiagMatrix.cc, liboctave/array/fCDiagMatrix.h, liboctave/array/fCMatrix.cc, liboctave/array/fCMatrix.h, liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h, liboctave/array/fCRowVector.cc, liboctave/array/fCRowVector.h, liboctave/array/fColVector.h, liboctave/array/fDiagMatrix.cc, liboctave/array/fDiagMatrix.h, liboctave/array/fMatrix.cc, liboctave/array/fMatrix.h, liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h, liboctave/array/fRowVector.h, liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h, liboctave/array/intNDArray.cc, liboctave/array/intNDArray.h, liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/quit.cc, liboctave/numeric/CmplxCHOL.cc, liboctave/numeric/CmplxCHOL.h, liboctave/numeric/CmplxGEPBAL.cc, liboctave/numeric/CmplxGEPBAL.h, liboctave/numeric/CmplxHESS.h, liboctave/numeric/CmplxLU.cc, liboctave/numeric/CmplxLU.h, liboctave/numeric/CmplxQR.cc, liboctave/numeric/CmplxQRP.cc, liboctave/numeric/CmplxQRP.h, liboctave/numeric/CmplxSCHUR.h, liboctave/numeric/CmplxSVD.cc, liboctave/numeric/CmplxSVD.h, liboctave/numeric/CollocWt.h, liboctave/numeric/DAE.h, liboctave/numeric/DAEFunc.h, liboctave/numeric/DAERT.h, liboctave/numeric/DAERTFunc.h, liboctave/numeric/DASPK.cc, liboctave/numeric/DASRT.cc, liboctave/numeric/DASRT.h, liboctave/numeric/DASSL.cc, liboctave/numeric/DET.h, liboctave/numeric/EIG.cc, liboctave/numeric/EIG.h, liboctave/numeric/LSODE.cc, liboctave/numeric/ODE.h, liboctave/numeric/ODEFunc.h, liboctave/numeric/ODES.h, liboctave/numeric/ODESFunc.h, liboctave/numeric/Quad.cc, liboctave/numeric/Quad.h, liboctave/numeric/SparseCmplxCHOL.h, liboctave/numeric/SparseCmplxLU.cc, liboctave/numeric/SparseCmplxLU.h, liboctave/numeric/SparseCmplxQR.cc, liboctave/numeric/SparseCmplxQR.h, liboctave/numeric/SparseQR.cc, liboctave/numeric/SparseQR.h, liboctave/numeric/SparsedbleCHOL.h, liboctave/numeric/SparsedbleLU.cc, liboctave/numeric/SparsedbleLU.h, liboctave/numeric/base-aepbal.h, liboctave/numeric/base-dae.h, liboctave/numeric/base-de.h, liboctave/numeric/base-lu.cc, liboctave/numeric/base-lu.h, liboctave/numeric/base-min.h, liboctave/numeric/base-qr.h, liboctave/numeric/bsxfun.h, liboctave/numeric/dbleCHOL.cc, liboctave/numeric/dbleCHOL.h, liboctave/numeric/dbleGEPBAL.h, liboctave/numeric/dbleHESS.h, liboctave/numeric/dbleLU.cc, liboctave/numeric/dbleLU.h, liboctave/numeric/dbleQR.cc, liboctave/numeric/dbleQRP.cc, liboctave/numeric/dbleQRP.h, liboctave/numeric/dbleSCHUR.cc, liboctave/numeric/dbleSCHUR.h, liboctave/numeric/dbleSVD.cc, liboctave/numeric/dbleSVD.h, liboctave/numeric/eigs-base.cc, liboctave/numeric/fCmplxAEPBAL.cc, liboctave/numeric/fCmplxAEPBAL.h, liboctave/numeric/fCmplxCHOL.cc, liboctave/numeric/fCmplxCHOL.h, liboctave/numeric/fCmplxGEPBAL.cc, liboctave/numeric/fCmplxGEPBAL.h, liboctave/numeric/fCmplxHESS.h, liboctave/numeric/fCmplxLU.cc, liboctave/numeric/fCmplxLU.h, liboctave/numeric/fCmplxQR.cc, liboctave/numeric/fCmplxQR.h, liboctave/numeric/fCmplxQRP.cc, liboctave/numeric/fCmplxQRP.h, liboctave/numeric/fCmplxSCHUR.cc, liboctave/numeric/fCmplxSCHUR.h, liboctave/numeric/fCmplxSVD.h, liboctave/numeric/fEIG.cc, liboctave/numeric/fEIG.h, liboctave/numeric/floatCHOL.cc, liboctave/numeric/floatCHOL.h, liboctave/numeric/floatGEPBAL.cc, liboctave/numeric/floatGEPBAL.h, liboctave/numeric/floatHESS.h, liboctave/numeric/floatLU.cc, liboctave/numeric/floatLU.h, liboctave/numeric/floatQR.cc, liboctave/numeric/floatQRP.cc, liboctave/numeric/floatQRP.h, liboctave/numeric/floatSCHUR.cc, liboctave/numeric/floatSCHUR.h, liboctave/numeric/floatSVD.cc, liboctave/numeric/floatSVD.h, liboctave/numeric/lo-mappers.cc, liboctave/numeric/lo-mappers.h, liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h, liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-fftw.cc, liboctave/numeric/oct-fftw.h, liboctave/numeric/oct-norm.cc, liboctave/numeric/oct-rand.cc, liboctave/numeric/oct-rand.h, liboctave/numeric/randgamma.c, liboctave/numeric/randgamma.h, liboctave/numeric/randmtzig.c, liboctave/numeric/randpoisson.c, liboctave/numeric/randpoisson.h, liboctave/numeric/sparse-base-chol.h, liboctave/numeric/sparse-base-lu.h, liboctave/numeric/sparse-dmsolve.cc, liboctave/operators/Sparse-diag-op-defs.h, liboctave/operators/Sparse-op-defs.h, liboctave/operators/mx-inlines.cc, liboctave/system/dir-ops.h, liboctave/system/file-ops.cc, liboctave/system/file-stat.cc, liboctave/system/file-stat.h, liboctave/system/lo-sysdep.cc, liboctave/system/lo-sysdep.h, liboctave/system/mach-info.cc, liboctave/system/mach-info.h, liboctave/system/oct-env.cc, liboctave/system/oct-group.cc, liboctave/system/oct-syscalls.cc, liboctave/system/oct-syscalls.h, liboctave/system/oct-time.h, liboctave/system/tempname.c, liboctave/util/action-container.h, liboctave/util/base-list.h, liboctave/util/cmd-edit.cc, liboctave/util/cmd-edit.h, liboctave/util/cmd-hist.cc, liboctave/util/cmd-hist.h, liboctave/util/data-conv.cc, liboctave/util/data-conv.h, liboctave/util/kpse.cc, liboctave/util/lo-array-gripes.cc, liboctave/util/lo-cieee.c, liboctave/util/lo-regexp.cc, liboctave/util/lo-utils.cc, liboctave/util/oct-alloc.cc, liboctave/util/oct-base64.cc, liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h, liboctave/util/oct-glob.cc, liboctave/util/oct-inttypes.cc, liboctave/util/oct-inttypes.h, liboctave/util/oct-locbuf.cc, liboctave/util/oct-locbuf.h, liboctave/util/oct-mem.h, liboctave/util/oct-mutex.cc, liboctave/util/oct-refcount.h, liboctave/util/oct-shlib.cc, liboctave/util/oct-shlib.h, liboctave/util/oct-sort.cc, liboctave/util/oct-sort.h, liboctave/util/pathsearch.cc, liboctave/util/pathsearch.h, liboctave/util/sparse-util.cc, liboctave/util/str-vec.cc, liboctave/util/str-vec.h, liboctave/util/unwind-prot.h, liboctave/util/url-transfer.cc, liboctave/util/url-transfer.h: Use GNU style coding conventions.

author	Rik <rik@octave.org>
date	Sat, 26 Oct 2013 18:57:05 -0700
parents	d63878346099
children	4197fc428c7d

rev	line source
7019 4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	1 /*
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	2
17744 d63878346099 maint: Update copyright notices for release. John W. Eaton <jwe@octave.org> parents: 15271 diff changeset	3 Copyright (C) 2006-2013 John W. Eaton
7019 4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	4
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	5 This file is part of Octave.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	6
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 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: 5742 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: 5742 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: 5742 diff changeset	10 option) any later version.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	11
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 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: 5742 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: 5742 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: 5742 diff changeset	15 for more details.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	16
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 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: 5742 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: 5742 diff changeset	19 <http://www.gnu.org/licenses/>.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	20
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	21 */
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	22
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	23 /* Original version written by Paul Kienzle distributed as free
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 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
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	26 /*
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	27
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	28 double randg (a)
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	29 void fill_randg (a,n,x)
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	30
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	31 Generate a series of standard gamma distributions.
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 See: Marsaglia G and Tsang W (2000), "A simple method for generating
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	34 gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	35
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	36 Needs the following defines:
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	37 * NAN: value to return for Not-A-Number
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	38 * RUNI: uniform generator on (0,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	39 * RNOR: normal generator
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	40 * REXP: exponential generator, or -log(RUNI) if one isn't available
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	41 * INFINITE: function to test whether a value is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	42
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	43 Test using:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	44 mean = a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	45 variance = a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	46 skewness = 2/sqrt(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	47 kurtosis = 3 + 6/sqrt(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	48
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	49 Note that randg can be used to generate many distributions:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	50
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	51 gamma(a,b) for a>0, b>0 (from R)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	52 r = b*randg(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	53 beta(a,b) for a>0, b>0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	54 r1 = randg(a,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	55 r = r1 / (r1 + randg(b,1))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	56 Erlang(a,n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	57 r = a*randg(n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	58 chisq(df) for df>0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	59 r = 2*randg(df/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	60 t(df) for 0<df<inf (use randn if df is infinite)
14846 460a3c6d8bf1 maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists. Rik <octave@nomad.inbox5.com> parents: 14655 diff changeset	61 r = randn () / sqrt(2*randg(df/2)/df)
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	62 F(n1,n2) for 0<n1, 0<n2
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	63 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	64 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	65 r = r1 / r2
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	66 negative binonial (n, p) for n>0, 0<p<=1
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	67 r = randp((1-p)/p * randg(n))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	68 (from R, citing Devroye(1986), Non-Uniform Random Variate Generation)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	69 non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	70 r = randp(L/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	71 r(r>0) = 2*randg(r(r>0))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	72 r(df>0) += 2*randg(df(df>0)/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	73 (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	74 Dirichlet(a1,...,ak) for ai > 0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	75 r = (randg(a1),...,randg(ak))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	76 r = r / sum(r)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	77 (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis)
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
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	80 #if defined (HAVE_CONFIG_H)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	81 #include <config.h>
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	82 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	83
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	84 #include <stdio.h>
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	85
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	86 #include "lo-ieee.h"
7231 2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe] jwe parents: 7019 diff changeset	87 #include "lo-math.h"
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	88 #include "randmtzig.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	89 #include "randgamma.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	90
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	91 #undef NAN
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	92 #define NAN octave_NaN
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	93 #define INFINITE lo_ieee_isinf
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	94 #define RUNI oct_randu()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	95 #define RNOR oct_randn()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	96 #define REXP oct_rande()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	97
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	98 void
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	99 oct_fill_randg (double a, octave_idx_type n, double *r)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	100 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	101 octave_idx_type i;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	102 /* If a < 1, start by generating gamma (1+a) */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	103 const double d = (a < 1. ? 1.+a : a) - 1./3.;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	104 const double c = 1./sqrt (9.*d);
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 /* Handle invalid cases */
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	107 if (a <= 0 \|\| INFINITE(a))
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	108 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	109 for (i=0; i < n; i++)
10317 42d098307c30 untabify additional source files John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	110 r[i] = NAN;
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	111 return;
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
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	114 for (i=0; i < n; i++)
5742 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 double x, xsq, v, u;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	117 restart:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	118 x = RNOR;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	119 v = (1+c*x);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	120 v = vv;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	121 if (v <= 0)
10317 42d098307c30 untabify additional source files John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	122 goto restart; /* rare, so don't bother moving up */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	123 u = RUNI;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	124 xsq = x*x;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	125 if (u >= 1.-0.0331xsqxsq && log (u) >= 0.5xsq + d(1-v+log (v)))
10317 42d098307c30 untabify additional source files John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	126 goto restart;
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	127 r[i] = d*v;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	128 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	129 if (a < 1)
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	130 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	131 /* Use gamma(a) = gamma(1+a)U^(1/a) /
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	132 /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	133 for (i = 0; i < n; i++)
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	134 r[i] *= exp (-REXP/a);
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	135 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	136 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	137
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	138 double
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	139 oct_randg (double a)
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 double ret;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	142 oct_fill_randg (a,1,&ret);
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	143 return ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	144 }
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	145
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	146 #undef NAN
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	147 #undef RUNI
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	148 #undef RNOR
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	149 #undef REXP
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	150 #define NAN octave_Float_NaN
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	151 #define RUNI oct_float_randu()
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	152 #define RNOR oct_float_randn()
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	153 #define REXP oct_float_rande()
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	154
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	155 void
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	156 oct_fill_float_randg (float a, octave_idx_type n, float *r)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	157 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	158 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	159 /* If a < 1, start by generating gamma(1+a) */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	160 const float d = (a < 1. ? 1.+a : a) - 1./3.;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	161 const float c = 1./sqrt (9.*d);
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	162
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	163 /* Handle invalid cases */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	164 if (a <= 0 \|\| INFINITE(a))
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	165 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	166 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	167 r[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	168 return;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	169 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	170
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	171 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	172 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	173 float x, xsq, v, u;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	174 frestart:
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	175 x = RNOR;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	176 v = (1+c*x);
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	177 v = vv;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	178 if (v <= 0)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	179 goto frestart; /* rare, so don't bother moving up */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	180 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	181 xsq = x*x;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	182 if (u >= 1.-0.0331xsqxsq && log (u) >= 0.5xsq + d(1-v+log (v)))
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	183 goto frestart;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	184 r[i] = d*v;
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	185 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	186 if (a < 1)
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	187 {
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	188 /* Use gamma(a) = gamma(1+a)U^(1/a) /
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	189 /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	190 for (i = 0; i < n; i++)
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	191 r[i] *= exp (-REXP/a);
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	192 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	193 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	194
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	195 float
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	196 oct_float_randg (float a)
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	197 {
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	198 float ret;
15018 3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/. Rik <rik@octave.org> parents: 14846 diff changeset	199 oct_fill_float_randg (a,1,&ret);
14655 43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) David Bateman <dbateman@free.fr> parents: 14138 diff changeset	200 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	201 }

Mercurial > hg > octave-nkf

annotate liboctave/numeric/randgamma.c @ 17769:49a5a4be04a1