+++ /dev/null
-/* 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;
-}
projects / ghc-hetmet.git / blobdiff