diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/grandom.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/grandom.c | 40 |
1 files changed, 22 insertions, 18 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/grandom.c b/Build/source/libs/mpfr/mpfr-src/src/grandom.c index 37f58cdc85a..03e1b49edaf 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/grandom.c +++ b/Build/source/libs/mpfr/mpfr-src/src/grandom.c @@ -3,7 +3,7 @@ distribution and round it to the precision of rop1, rop2 according to the given rounding mode. -Copyright 2011-2019 Free Software Foundation, Inc. +Copyright 2011-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -40,13 +40,9 @@ mpfr_grandom (mpfr_ptr rop1, mpfr_ptr rop2, gmp_randstate_t rstate, inex2 = inex1 = 0; if (rop2 == NULL) /* only one output requested. */ - { - tprec0 = MPFR_PREC (rop1); - } + tprec0 = MPFR_PREC (rop1); else - { - tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2)); - } + tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2)); tprec0 += 11; @@ -87,21 +83,29 @@ mpfr_grandom (mpfr_ptr rop1, mpfr_ptr rop2, gmp_randstate_t rstate, mpz_mul (b, yp, yp); mpz_add (s, a, b); } - while (mpz_sizeinbase (s, 2) > tprec * 2); /* x^2 + y^2 <= 2^{2tprec} */ + while (mpz_sizeinbase (s, 2) > tprec * 2); + + /* now s = x^2 + y^2 < 2^{2tprec} */ for (;;) { - /* FIXME: compute s as s += 2x + 2y + 2 */ - mpz_add_ui (a, xp, 1); - mpz_add_ui (b, yp, 1); - mpz_mul (a, a, a); - mpz_mul (b, b, b); - mpz_add (s, a, b); - if ((mpz_sizeinbase (s, 2) <= 2 * tprec) || - ((mpz_sizeinbase (s, 2) == 2 * tprec + 1) && - (mpz_scan1 (s, 0) == 2 * tprec))) + /* we compute (xp+1)^2 + (yp+1)^2 as s + 2xp + 2yp + 2 */ + mpz_addmul_ui (s, xp, 2); + mpz_addmul_ui (s, yp, 2); + mpz_add_ui (s, s, 2); + /* The case s = 2^(2*tprec) is not possible: + (a) if xp and yp have different parities, s is odd + (b) if xp and yp are even, (xp+1)^2 and (yp+1)^2 are 1 mod 4, + thus s = 2 mod 4 (and tprec >= 1); + (c) if xp and yp are odd, if we note x = xp+1, y = yp+1 and + p = tprec, we would have x^2 + y^2 = 2^(2p) with x and y even + 0 < x, y <= 2^p, thus if x' = x/2, y' = y/2 and p'=p-1, + we would have x'^2 + y'^2 = 2^(2p') with + 0 < x', y' <= 2^p', and we conclude by induction. */ + if (mpz_sizeinbase (s, 2) <= 2 * tprec) goto yeepee; - /* Extend by 32 bits */ + /* Extend by 32 bits: for tprec=12, the probability we get here + is 8191/13180825, i.e., about 0.000621 */ mpz_mul_2exp (xp, xp, 32); mpz_mul_2exp (yp, yp, 32); mpz_urandomb (x, rstate, 32); |