Mercurial > hg > octave-nkf
diff liboctave/numeric/randpoisson.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 |
line wrap: on
line diff
--- a/liboctave/numeric/randpoisson.c +++ b/liboctave/numeric/randpoisson.c @@ -80,13 +80,14 @@ static double flogfak (double k) { -#define C0 9.18938533204672742e-01 -#define C1 8.33333333333333333e-02 -#define C3 -2.77777777777777778e-03 -#define C5 7.93650793650793651e-04 -#define C7 -5.95238095238095238e-04 +#define C0 9.18938533204672742e-01 +#define C1 8.33333333333333333e-02 +#define C3 -2.77777777777777778e-03 +#define C5 7.93650793650793651e-04 +#define C7 -5.95238095238095238e-04 - static double logfak[30L] = { + static double logfak[30L] = + { 0.00000000000000000, 0.00000000000000000, 0.69314718055994531, 1.79175946922805500, 3.17805383034794562, 4.78749174278204599, 6.57925121201010100, 8.52516136106541430, 10.60460290274525023, @@ -152,12 +153,12 @@ static double pprsc (double my) { - static double my_last = -1.0; - static double m, k2, k4, k1, k5; - static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm, - f1, f2, f4, f5, p1, p2, p3, p4, p5, p6; - double Dk, X, Y; - double Ds, U, V, W; + static double my_last = -1.0; + static double m, k2, k4, k1, k5; + static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm, + f1, f2, f4, f5, p1, p2, p3, p4, p5, p6; + double Dk, X, Y; + double Ds, U, V, W; if (my != my_last) { /* set-up */ @@ -322,50 +323,53 @@ size_t i = n; t[0] = P = exp (-lambda); - for (tableidx = 1; tableidx <= intlambda; tableidx++) { - P = P*lambda/(double)tableidx; - t[tableidx] = t[tableidx-1] + P; - } - - while (i-- > 0) { - double u = RUNI; - - /* If u > 0.458 we know we can jump to floor(lambda) before - * comparing (this observation is based on Stadlober's winrand - * code). For lambda >= 1, this will be a win. Lambda < 1 - * is already fast, so adding an extra comparison is not a - * problem. */ - int k = (u > 0.458 ? intlambda : 0); - - /* We aren't using a for loop here because when we find the - * right k we want to jump to the next iteration of the - * outer loop, and the continue statement will only work for - * the inner loop. */ - nextk: - if ( u <= t[k] ) { - p[i] = (double) k; - continue; - } - if (++k < tableidx) goto nextk; - - /* We only need high values of the table very rarely so we - * don't automatically compute the entire table. */ - while (tableidx < TABLESIZE) { + for (tableidx = 1; tableidx <= intlambda; tableidx++) + { P = P*lambda/(double)tableidx; t[tableidx] = t[tableidx-1] + P; - /* Make sure we converge to 1.0 just in case u is uniform - * on [0,1] rather than [0,1). */ - if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; - tableidx++; - if (u <= t[tableidx-1]) break; } - /* We are assuming that the table size is big enough here. - * This should be true even if RUNI is returning values in - * the range [0,1] rather than [0,1). - */ - p[i] = (double)(tableidx-1); - } + while (i-- > 0) + { + double u = RUNI; + + /* If u > 0.458 we know we can jump to floor(lambda) before + * comparing (this observation is based on Stadlober's winrand + * code). For lambda >= 1, this will be a win. Lambda < 1 + * is already fast, so adding an extra comparison is not a + * problem. */ + int k = (u > 0.458 ? intlambda : 0); + + /* We aren't using a for loop here because when we find the + * right k we want to jump to the next iteration of the + * outer loop, and the continue statement will only work for + * the inner loop. */ + nextk: + if (u <= t[k]) + { + p[i] = (double) k; + continue; + } + if (++k < tableidx) goto nextk; + + /* We only need high values of the table very rarely so we + * don't automatically compute the entire table. */ + while (tableidx < TABLESIZE) + { + P = P*lambda/(double)tableidx; + t[tableidx] = t[tableidx-1] + P; + /* Make sure we converge to 1.0 just in case u is uniform + * on [0,1] rather than [0,1). */ + if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; + tableidx++; + if (u <= t[tableidx-1]) break; + } + + /* We are assuming that the table size is big enough here. + * This should be true even if RUNI is returning values in + * the range [0,1] rather than [0,1). */ + p[i] = (double)(tableidx-1); + } } static void @@ -381,31 +385,35 @@ size_t i = n; t[0] = P = exp (-lambda); - for (tableidx = 1; tableidx <= intlambda; tableidx++) { - P = P*lambda/(double)tableidx; - t[tableidx] = t[tableidx-1] + P; - } - - while (i-- > 0) { - double u = RUNI; - int k = (u > 0.458 ? intlambda : 0); - nextk: - if ( u <= t[k] ) { - p[i] = (float) k; - continue; - } - if (++k < tableidx) goto nextk; - - while (tableidx < TABLESIZE) { + for (tableidx = 1; tableidx <= intlambda; tableidx++) + { P = P*lambda/(double)tableidx; t[tableidx] = t[tableidx-1] + P; - if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; - tableidx++; - if (u <= t[tableidx-1]) break; } - p[i] = (float)(tableidx-1); - } + while (i-- > 0) + { + double u = RUNI; + int k = (u > 0.458 ? intlambda : 0); + nextk: + if (u <= t[k]) + { + p[i] = (float) k; + continue; + } + if (++k < tableidx) goto nextk; + + while (tableidx < TABLESIZE) + { + P = P*lambda/(double)tableidx; + t[tableidx] = t[tableidx-1] + P; + if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; + tableidx++; + if (u <= t[tableidx-1]) break; + } + + p[i] = (float)(tableidx-1); + } } /* From Press, et al., Numerical Recipes */ @@ -420,14 +428,16 @@ for (i = 0; i < n; i++) { double y, em, t; - do { - do { - y = tan (M_PI*RUNI); - em = sq * y + lambda; - } while (em < 0.0); - em = floor (em); - t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); - } while (RUNI > t); + do + { + do + { + y = tan (M_PI*RUNI); + em = sq * y + lambda; + } while (em < 0.0); + em = floor (em); + t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); + } while (RUNI > t); p[i] = em; } } @@ -444,14 +454,16 @@ for (i = 0; i < n; i++) { double y, em, t; - do { - do { - y = tan (M_PI*RUNI); - em = sq * y + lambda; - } while (em < 0.0); - em = floor (em); - t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); - } while (RUNI > t); + do + { + do + { + y = tan (M_PI*RUNI); + em = sq * y + lambda; + } while (em < 0.0); + em = floor (em); + t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); + } while (RUNI > t); p[i] = em; } } @@ -503,28 +515,36 @@ { double ret; if (L < 0.0) ret = NAN; - else if (L <= 12.0) { - /* From Press, et al. Numerical recipes */ - double g = exp (-L); - int em = -1; - double t = 1.0; - do { - ++em; - t *= RUNI; - } while (t > g); - ret = em; - } else if (L <= 1e8) { - /* numerical recipes */ - poisson_rejection (L, &ret, 1); - } else if (INFINITE(L)) { - /* FIXME R uses NaN, but the normal approx. suggests that as - * limit should be inf. Which is correct? */ - ret = NAN; - } else { - /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = floor (RNOR*sqrt (L) + L + 0.5); - if (ret < 0.0) ret = 0.0; /* will probably never happen */ - } + else if (L <= 12.0) + { + /* From Press, et al. Numerical recipes */ + double g = exp (-L); + int em = -1; + double t = 1.0; + do + { + ++em; + t *= RUNI; + } while (t > g); + ret = em; + } + else if (L <= 1e8) + { + /* numerical recipes */ + poisson_rejection (L, &ret, 1); + } + else if (INFINITE(L)) + { + /* FIXME R uses NaN, but the normal approx. suggests that as + * limit should be inf. Which is correct? */ + ret = NAN; + } + else + { + /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ + ret = floor (RNOR*sqrt (L) + L + 0.5); + if (ret < 0.0) ret = 0.0; /* will probably never happen */ + } return ret; } @@ -568,27 +588,35 @@ double L = FL; float ret; if (L < 0.0) ret = NAN; - else if (L <= 12.0) { - /* From Press, et al. Numerical recipes */ - double g = exp (-L); - int em = -1; - double t = 1.0; - do { - ++em; - t *= RUNI; - } while (t > g); - ret = em; - } else if (L <= 1e8) { - /* numerical recipes */ - poisson_rejection_float (L, &ret, 1); - } else if (INFINITE(L)) { - /* FIXME R uses NaN, but the normal approx. suggests that as - * limit should be inf. Which is correct? */ - ret = NAN; - } else { - /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = floor (RNOR*sqrt (L) + L + 0.5); - if (ret < 0.0) ret = 0.0; /* will probably never happen */ - } + else if (L <= 12.0) + { + /* From Press, et al. Numerical recipes */ + double g = exp (-L); + int em = -1; + double t = 1.0; + do + { + ++em; + t *= RUNI; + } while (t > g); + ret = em; + } + else if (L <= 1e8) + { + /* numerical recipes */ + poisson_rejection_float (L, &ret, 1); + } + else if (INFINITE(L)) + { + /* FIXME R uses NaN, but the normal approx. suggests that as + * limit should be inf. Which is correct? */ + ret = NAN; + } + else + { + /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ + ret = floor (RNOR*sqrt (L) + L + 0.5); + if (ret < 0.0) ret = 0.0; /* will probably never happen */ + } return ret; }