diff options
Diffstat (limited to 'rts/gmp/mpz/root.c')
| -rw-r--r-- | rts/gmp/mpz/root.c | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/rts/gmp/mpz/root.c b/rts/gmp/mpz/root.c new file mode 100644 index 0000000000..0920bf22d3 --- /dev/null +++ b/rts/gmp/mpz/root.c @@ -0,0 +1,183 @@ +/* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)). + Return an indication if the result is exact. + +Copyright (C) 1999, 2000 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +/* Naive implementation of nth root extraction. It would probably be a + better idea to use a division-free Newton iteration. It is insane + to use full precision from iteration 1. The mpz_scan1 trick compensates + to some extent. It would be natural to avoid representing the low zero + bits mpz_scan1 is counting, and at the same time call mpn directly. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +int +#if __STDC__ +mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth) +#else +mpz_root (r, c, nth) + mpz_ptr r; + mpz_srcptr c; + unsigned long int nth; +#endif +{ + mpz_t x, t0, t1, t2; + __mpz_struct ccs, *cc = &ccs; + unsigned long int nbits; + int bit; + int exact; + int i; + unsigned long int lowz; + unsigned long int rl; + + /* even roots of negatives provoke an exception */ + if (mpz_sgn (c) < 0 && (nth & 1) == 0) + SQRT_OF_NEGATIVE; + + /* root extraction interpreted as c^(1/nth) means a zeroth root should + provoke a divide by zero, do this even if c==0 */ + if (nth == 0) + DIVIDE_BY_ZERO; + + if (mpz_sgn (c) == 0) + { + if (r != NULL) + mpz_set_ui (r, 0); + return 1; /* exact result */ + } + + PTR(cc) = PTR(c); + SIZ(cc) = ABSIZ(c); + + nbits = (mpz_sizeinbase (cc, 2) - 1) / nth; + if (nbits == 0) + { + if (r != NULL) + mpz_set_ui (r, 1); + if (mpz_sgn (c) < 0) + { + if (r != NULL) + SIZ(r) = -SIZ(r); + return mpz_cmp_si (c, -1L) == 0; + } + return mpz_cmp_ui (c, 1L) == 0; + } + + mpz_init (x); + mpz_init (t0); + mpz_init (t1); + mpz_init (t2); + + /* Create a one-bit approximation. */ + mpz_set_ui (x, 0); + mpz_setbit (x, nbits); + + /* Make the approximation better, one bit at a time. This odd-looking + termination criteria makes large nth get better initial approximation, + which avoids slow convergence for such values. */ + bit = nbits - 1; + for (i = 1; (nth >> i) != 0; i++) + { + mpz_setbit (x, bit); + mpz_tdiv_q_2exp (t0, x, bit); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, bit * nth); + if (mpz_cmp (cc, t1) < 0) + mpz_clrbit (x, bit); + + bit--; /* check/set next bit */ + if (bit < 0) + { + /* We're done. */ + mpz_pow_ui (t1, x, nth); + goto done; + } + } + mpz_setbit (x, bit); + mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2); + +#if DEBUG + /* Check that the starting approximation is >= than the root. */ + mpz_pow_ui (t1, x, nth); + if (mpz_cmp (cc, t1) >= 0) + abort (); +#endif + + mpz_add_ui (x, x, 1); + + /* Main loop */ + do + { + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth - 1); + mpz_mul_2exp (t1, t1, lowz * (nth - 1)); + mpz_tdiv_q (t2, cc, t1); + mpz_sub (t2, x, t2); + rl = mpz_tdiv_q_ui (t2, t2, nth); + mpz_sub (x, x, t2); + } + while (mpz_sgn (t2) != 0); + + /* If we got a non-zero remainder in the last division, we know our root + is too large. */ + mpz_sub_ui (x, x, (mp_limb_t) (rl != 0)); + + /* Adjustment loop. If we spend more care on rounding in the loop above, + we could probably get rid of this, or greatly simplify it. */ + { + int bad = 0; + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, lowz * nth); + while (mpz_cmp (cc, t1) < 0) + { + bad++; + if (bad > 2) + abort (); /* abort if our root is far off */ + mpz_sub_ui (x, x, 1); + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, lowz * nth); + } + } + + done: + exact = mpz_cmp (t1, cc) == 0; + + if (r != NULL) + { + mpz_set (r, x); + if (mpz_sgn (c) < 0) + SIZ(r) = -SIZ(r); + } + + mpz_clear (t2); + mpz_clear (t1); + mpz_clear (t0); + mpz_clear (x); + + return exact; +} |