diff options
author | tege <tege@gmplib.org> | 2002-05-08 20:58:23 +0200 |
---|---|---|
committer | tege <tege@gmplib.org> | 2002-05-08 20:58:23 +0200 |
commit | 170a6691239fc423cb66cb46211e1ae30bf817e2 (patch) | |
tree | 2721ffa397de84be546766332f26c8bc1d484330 /mpz/root.c | |
parent | d239457d3ae0e7f590b66ae1340325156f29d0d7 (diff) | |
download | gmp-170a6691239fc423cb66cb46211e1ae30bf817e2.tar.gz |
Rewrite to use mpn_rootrem.
Diffstat (limited to 'mpz/root.c')
-rw-r--r-- | mpz/root.c | 143 |
1 files changed, 31 insertions, 112 deletions
diff --git a/mpz/root.c b/mpz/root.c index 73c4c49ea..677ac5bad 100644 --- a/mpz/root.c +++ b/mpz/root.c @@ -33,19 +33,19 @@ MA 02111-1307, USA. */ #include "longlong.h" int -mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth) +mpz_root (mpz_ptr r, mpz_srcptr u, unsigned long int nth) { - mpz_t x, t0, t1, t2; - __mpz_struct ccs, *cc = &ccs; - unsigned long int nbits; - int bit; + mp_ptr rootp, up; + mp_size_t us, un, rootn; int exact; - int i; - unsigned long int lowz; - unsigned long int rl; + unsigned int cnt; + unsigned long int rootnb, unb; + + up = PTR(u); + us = SIZ(u); /* even roots of negatives provoke an exception */ - if (mpz_sgn (c) < 0 && (nth & 1) == 0) + if (us < 0 && (nth & 1) == 0) SQRT_OF_NEGATIVE; /* root extraction interpreted as c^(1/nth) means a zeroth root should @@ -53,125 +53,44 @@ mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth) if (nth == 0) DIVIDE_BY_ZERO; - if (mpz_sgn (c) == 0) + if (us == 0) { if (r != NULL) - mpz_set_ui (r, 0); + SIZ(r) = 0; return 1; /* exact result */ } - PTR(cc) = PTR(c); - SIZ(cc) = ABSIZ(c); + un = ABS (us); + + count_leading_zeros (cnt, up[un - 1]); + unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS; + rootnb = (unb + nth - 1) / nth; + rootn = (rootnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; - nbits = (mpz_sizeinbase (cc, 2) - 1) / nth; - if (nbits == 0) + if (r != NULL) { - 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; + rootp = MPZ_REALLOC (r, rootn + 1); + up = PTR(u); } - - 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++) + else { - 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; - } + rootp = __GMP_ALLOCATE_FUNC_LIMBS (rootn + 1); } - 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 + if (nth == 1) { - 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); + MPN_COPY (rootp, up, un); + exact = 1; } - 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, (unsigned long) (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) + else { - mpz_set (r, x); - if (mpz_sgn (c) < 0) - SIZ(r) = -SIZ(r); + exact = 0 == mpn_rootrem (rootp, NULL, up, un, nth); } - mpz_clear (t2); - mpz_clear (t1); - mpz_clear (t0); - mpz_clear (x); + if (r != NULL) + SIZ(r) = us >= 0 ? rootn : -rootn; + else + __GMP_FREE_FUNC_LIMBS (rootp, rootn + 1); return exact; } |