Mercurial > octave
changeset 21746:ac8b8bdae98f
use gnulib:: tag in rand functions as needed
* randgamma.cc, randmtzig.cc, randpoisson.cc: Tag functions with
gnulib:: namespace id as needed.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 19 May 2016 14:40:35 -0400 |
| parents | bf1121302404 |
| children | 61f3575250e4 |
| files | liboctave/numeric/randgamma.cc liboctave/numeric/randmtzig.cc liboctave/numeric/randpoisson.cc |
| diffstat | 3 files changed, 41 insertions(+), 41 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/numeric/randgamma.cc Thu May 19 13:59:28 2016 -0400 +++ b/liboctave/numeric/randgamma.cc Thu May 19 14:40:35 2016 -0400 @@ -118,7 +118,7 @@ goto restart; /* rare, so don't bother moving up */ u = RUNI; xsq = x*x; - if (u >= 1.-0.0331*xsq*xsq && log (u) >= 0.5*xsq + d*(1-v+log (v))) + if (u >= 1.-0.0331*xsq*xsq && gnulib::log (u) >= 0.5*xsq + d*(1-v+gnulib::log (v))) goto restart; r[i] = d*v; } @@ -173,7 +173,7 @@ goto frestart; /* rare, so don't bother moving up */ u = RUNI; xsq = x*x; - if (u >= 1.-0.0331*xsq*xsq && log (u) >= 0.5*xsq + d*(1-v+log (v))) + if (u >= 1.-0.0331*xsq*xsq && gnulib::log (u) >= 0.5*xsq + d*(1-v+gnulib::log (v))) goto frestart; r[i] = d*v; }
--- a/liboctave/numeric/randmtzig.cc Thu May 19 13:59:28 2016 -0400 +++ b/liboctave/numeric/randmtzig.cc Thu May 19 14:40:35 2016 -0400 @@ -261,17 +261,17 @@ int n = 0; /* Look for entropy in /dev/urandom */ - FILE* urandom = fopen ("/dev/urandom", "rb"); + FILE* urandom = gnulib::fopen ("/dev/urandom", "rb"); if (urandom) { while (n < MT_N) { unsigned char word[4]; - if (fread (word, 4, 1, urandom) != 1) + if (gnulib::fread (word, 4, 1, urandom) != 1) break; entropy[n++] = word[0]+(word[1]<<8)+(word[2]<<16)+(static_cast<uint32_t>(word[3])<<24); } - fclose (urandom); + gnulib::fclose (urandom); } /* If there isn't enough entropy, gather some from various sources */ @@ -283,7 +283,7 @@ if (n < MT_N) { struct timeval tv; - if (gettimeofday (&tv, 0) != -1) + if (gnulib::gettimeofday (&tv, 0) != -1) entropy[n++] = tv.tv_usec; /* Fractional part of current time */ } #endif @@ -507,7 +507,7 @@ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) */ - x = sqrt (-2. * log (NOR_SECTION_AREA / x1 + fi[i+1])); + x = sqrt (-2. * gnulib::log (NOR_SECTION_AREA / x1 + fi[i+1])); ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA); wi[i] = x / NMANTISSA; fi[i] = exp (-0.5 * x * x); @@ -535,7 +535,7 @@ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus * need inverse operator of y = exp(-x) -> x = -ln(y) */ - x = - log (EXP_SECTION_AREA / x1 + fe[i+1]); + x = - gnulib::log (EXP_SECTION_AREA / x1 + fe[i+1]); ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA); we[i] = x / EMANTISSA; fe[i] = exp (-x); @@ -615,8 +615,8 @@ double xx, yy; do { - xx = - ZIGGURAT_NOR_INV_R * log (RANDU); - yy = - log (RANDU); + xx = - ZIGGURAT_NOR_INV_R * gnulib::log (RANDU); + yy = - gnulib::log (RANDU); } while ( yy+yy <= xx*xx); return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); @@ -646,7 +646,7 @@ * For the exponential tail, the method of Marsaglia[5] provides: * x = r - ln(U); */ - return ZIGGURAT_EXP_R - log (RANDU); + return ZIGGURAT_EXP_R - gnulib::log (RANDU); } else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x)) return x; @@ -698,7 +698,7 @@ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) */ - x = sqrt (-2. * log (NOR_SECTION_AREA / x1 + ffi[i+1])); + x = sqrt (-2. * gnulib::log (NOR_SECTION_AREA / x1 + ffi[i+1])); fki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA); fwi[i] = x / NMANTISSA; ffi[i] = exp (-0.5 * x * x); @@ -726,7 +726,7 @@ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus * need inverse operator of y = exp(-x) -> x = -ln(y) */ - x = - log (EXP_SECTION_AREA / x1 + ffe[i+1]); + x = - gnulib::log (EXP_SECTION_AREA / x1 + ffe[i+1]); fke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA); fwe[i] = x / EMANTISSA; ffe[i] = exp (-x); @@ -782,8 +782,8 @@ float xx, yy; do { - xx = - ZIGGURAT_NOR_INV_R * log (RANDU); - yy = - log (RANDU); + xx = - ZIGGURAT_NOR_INV_R * gnulib::log (RANDU); + yy = - gnulib::log (RANDU); } while ( yy+yy <= xx*xx); return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); @@ -813,7 +813,7 @@ * For the exponential tail, the method of Marsaglia[5] provides: * x = r - ln(U); */ - return ZIGGURAT_EXP_R - log (RANDU); + return ZIGGURAT_EXP_R - gnulib::log (RANDU); } else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x)) return x;
--- a/liboctave/numeric/randpoisson.cc Thu May 19 13:59:28 2016 -0400 +++ b/liboctave/numeric/randpoisson.cc Thu May 19 14:40:35 2016 -0400 @@ -102,7 +102,7 @@ { r = 1.0 / k; rr = r * r; - return ((k + 0.5)*log (k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7)))); + return ((k + 0.5)*gnulib::log (k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7)))); } else return (logfak[static_cast<int> (k)]); @@ -164,9 +164,9 @@ /* mode m, reflection points k2 and k4, and points k1 and k5, */ /* which delimit the centre region of h(x) */ - m = floor (my); + m = gnulib::floor (my); k2 = ceil (my - 0.5 - Ds); - k4 = floor (my - 0.5 + Ds); + k4 = gnulib::floor (my - 0.5 + Ds); k1 = k2 + k2 - m + 1L; k5 = k4 + k4 - m; @@ -181,11 +181,11 @@ r5 = my / (k5 + 1.0); /* reciprocal values of the scale parameters of exp. tail envelope */ - ll = log (r1); /* expon. tail left */ - lr = -log (r5); /* expon. tail right*/ + ll = gnulib::log (r1); /* expon. tail left */ + lr = -gnulib::log (r5); /* expon. tail right*/ /* Poisson constants, necessary for computing function values f(k) */ - l_my = log (my); + l_my = gnulib::log (my); c_pm = m * l_my - flogfak (m); /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */ @@ -213,14 +213,14 @@ /* immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1 */ - if ((V = U - p1) < 0.0) return (k2 + floor (U/f2)); + if ((V = U - p1) < 0.0) return (k2 + gnulib::floor (U/f2)); /* immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */ - if ((W = V / dl) < f1 ) return (k1 + floor (V/f1)); + if ((W = V / dl) < f1 ) return (k1 + gnulib::floor (V/f1)); /* computation of candidate X < k2, and its counterpart Y > k2 */ /* either squeeze-acceptance of X or acceptance-rejection of Y */ - Dk = floor (dl * RUNI) + 1.0; + Dk = gnulib::floor (dl * RUNI) + 1.0; if (W <= f2 - Dk * (f2 - f2/r2)) { /* quick accept of */ return (k2 - Dk); /* X = k2 - Dk */ @@ -240,14 +240,14 @@ { /* centre right */ /* immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4 */ - if ((V = U - p3) < 0.0) return (k4 - floor ((U - p2)/f4)); + if ((V = U - p3) < 0.0) return (k4 - gnulib::floor ((U - p2)/f4)); /* immediate acceptance region R4 = [k4+1, k5+1)*[0, f5) */ - if ((W = V / dr) < f5 ) return (k5 - floor (V/f5)); + if ((W = V / dr) < f5 ) return (k5 - gnulib::floor (V/f5)); /* computation of candidate X > k4, and its counterpart Y < k4 */ /* either squeeze-acceptance of X or acceptance-rejection of Y */ - Dk = floor (dr * RUNI) + 1.0; + Dk = gnulib::floor (dr * RUNI) + 1.0; if (W <= f4 - Dk * (f4 - f4*r4)) { /* quick accept of */ return (k4 + Dk); /* X = k4 + Dk */ @@ -268,7 +268,7 @@ W = RUNI; if (U < p5) { /* expon. tail left */ - Dk = floor (1.0 - log (W)/ll); + Dk = gnulib::floor (1.0 - gnulib::log (W)/ll); if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */ W *= (U - p4) * ll; /* W -- U(0, h(x)) */ if (W <= f1 - Dk * (f1 - f1/r1)) @@ -276,7 +276,7 @@ } else { /* expon. tail right*/ - Dk = floor (1.0 - log (W)/lr); + Dk = gnulib::floor (1.0 - gnulib::log (W)/lr); X = k5 + Dk; /* X >= k5 + 1 */ W *= (U - p5) * lr; /* W -- U(0, h(x)) */ if (W <= f5 - Dk * (f5 - f5*r5)) @@ -287,7 +287,7 @@ /* acceptance-rejection test of candidate X from the original area */ /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/ /* log f(X) = (X - m)*log(my) - log X! + log m! */ - if (log (W) <= X * l_my - flogfak (X) - c_pm) return (X); + if (gnulib::log (W) <= X * l_my - flogfak (X) - c_pm) return (X); } } /* ---- pprsc.c end ------ */ @@ -313,7 +313,7 @@ /* Precompute the table for the u up to and including 0.458. * We will almost certainly need it. */ - int intlambda = static_cast<int> (floor (lambda)); + int intlambda = static_cast<int> (gnulib::floor (lambda)); double P; int tableidx; size_t i = n; @@ -376,7 +376,7 @@ /* Precompute the table for the u up to and including 0.458. * We will almost certainly need it. */ - int intlambda = static_cast<int> (floor (lambda)); + int intlambda = static_cast<int> (gnulib::floor (lambda)); double P; int tableidx; size_t i = n; @@ -419,7 +419,7 @@ poisson_rejection (double lambda, double *p, size_t n) { double sq = sqrt (2.0*lambda); - double alxm = log (lambda); + double alxm = gnulib::log (lambda); double g = lambda*alxm - LGAMMA(lambda+1.0); size_t i; @@ -433,7 +433,7 @@ y = tan (M_PI*RUNI); em = sq * y + lambda; } while (em < 0.0); - em = floor (em); + em = gnulib::floor (em); t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); } while (RUNI > t); p[i] = em; @@ -445,7 +445,7 @@ poisson_rejection_float (double lambda, float *p, size_t n) { double sq = sqrt (2.0*lambda); - double alxm = log (lambda); + double alxm = gnulib::log (lambda); double g = lambda*alxm - LGAMMA(lambda+1.0); size_t i; @@ -459,7 +459,7 @@ y = tan (M_PI*RUNI); em = sq * y + lambda; } while (em < 0.0); - em = floor (em); + em = gnulib::floor (em); t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); } while (RUNI > t); p[i] = em; @@ -500,7 +500,7 @@ const double sqrtL = sqrt (L); for (i = 0; i < n; i++) { - p[i] = floor (RNOR*sqrtL + L + 0.5); + p[i] = gnulib::floor (RNOR*sqrtL + L + 0.5); if (p[i] < 0.0) p[i] = 0.0; /* will probably never happen */ } @@ -540,7 +540,7 @@ else { /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = floor (RNOR*sqrt (L) + L + 0.5); + ret = gnulib::floor (RNOR*sqrt (L) + L + 0.5); if (ret < 0.0) ret = 0.0; /* will probably never happen */ } return ret; @@ -572,7 +572,7 @@ const double sqrtL = sqrt (L); for (i = 0; i < n; i++) { - p[i] = floor (RNOR*sqrtL + L + 0.5); + p[i] = gnulib::floor (RNOR*sqrtL + L + 0.5); if (p[i] < 0.0) p[i] = 0.0; /* will probably never happen */ } @@ -613,7 +613,7 @@ else { /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = floor (RNOR*sqrt (L) + L + 0.5); + ret = gnulib::floor (RNOR*sqrt (L) + L + 0.5); if (ret < 0.0) ret = 0.0; /* will probably never happen */ } return ret;