Annotation of OpenXM_contrib/gmp/mpz/root.c, Revision 1.1
1.1 ! maekawa 1: /* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)).
! 2: Return an indication if the result is exact.
! 3:
! 4: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
! 5:
! 6: This file is part of the GNU MP Library.
! 7:
! 8: The GNU MP Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The GNU MP Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA. */
! 22:
! 23: /* Naive implementation of nth root extraction. It would probably be a
! 24: better idea to use a division-free Newton iteration. It is insane
! 25: to use full precision from iteration 1. The mpz_scan1 trick compensates
! 26: to some extent. It would be natural to avoid representing the low zero
! 27: bits mpz_scan1 is counting, and at the same time call mpn directly. */
! 28:
! 29: #include <stdio.h> /* for NULL */
! 30: #include "gmp.h"
! 31: #include "gmp-impl.h"
! 32: #include "longlong.h"
! 33:
! 34: int
! 35: #if __STDC__
! 36: mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth)
! 37: #else
! 38: mpz_root (r, c, nth)
! 39: mpz_ptr r;
! 40: mpz_srcptr c;
! 41: unsigned long int nth;
! 42: #endif
! 43: {
! 44: mpz_t x, t0, t1, t2;
! 45: __mpz_struct ccs, *cc = &ccs;
! 46: unsigned long int nbits;
! 47: int bit;
! 48: int exact;
! 49: int i;
! 50: unsigned long int lowz;
! 51: unsigned long int rl;
! 52:
! 53: /* even roots of negatives provoke an exception */
! 54: if (mpz_sgn (c) < 0 && (nth & 1) == 0)
! 55: SQRT_OF_NEGATIVE;
! 56:
! 57: /* root extraction interpreted as c^(1/nth) means a zeroth root should
! 58: provoke a divide by zero, do this even if c==0 */
! 59: if (nth == 0)
! 60: DIVIDE_BY_ZERO;
! 61:
! 62: if (mpz_sgn (c) == 0)
! 63: {
! 64: if (r != NULL)
! 65: mpz_set_ui (r, 0);
! 66: return 1; /* exact result */
! 67: }
! 68:
! 69: PTR(cc) = PTR(c);
! 70: SIZ(cc) = ABSIZ(c);
! 71:
! 72: nbits = (mpz_sizeinbase (cc, 2) - 1) / nth;
! 73: if (nbits == 0)
! 74: {
! 75: if (r != NULL)
! 76: mpz_set_ui (r, 1);
! 77: if (mpz_sgn (c) < 0)
! 78: {
! 79: if (r != NULL)
! 80: SIZ(r) = -SIZ(r);
! 81: return mpz_cmp_si (c, -1L) == 0;
! 82: }
! 83: return mpz_cmp_ui (c, 1L) == 0;
! 84: }
! 85:
! 86: mpz_init (x);
! 87: mpz_init (t0);
! 88: mpz_init (t1);
! 89: mpz_init (t2);
! 90:
! 91: /* Create a one-bit approximation. */
! 92: mpz_set_ui (x, 0);
! 93: mpz_setbit (x, nbits);
! 94:
! 95: /* Make the approximation better, one bit at a time. This odd-looking
! 96: termination criteria makes large nth get better initial approximation,
! 97: which avoids slow convergence for such values. */
! 98: bit = nbits - 1;
! 99: for (i = 1; (nth >> i) != 0; i++)
! 100: {
! 101: mpz_setbit (x, bit);
! 102: mpz_tdiv_q_2exp (t0, x, bit);
! 103: mpz_pow_ui (t1, t0, nth);
! 104: mpz_mul_2exp (t1, t1, bit * nth);
! 105: if (mpz_cmp (cc, t1) < 0)
! 106: mpz_clrbit (x, bit);
! 107:
! 108: bit--; /* check/set next bit */
! 109: if (bit < 0)
! 110: {
! 111: /* We're done. */
! 112: mpz_pow_ui (t1, x, nth);
! 113: goto done;
! 114: }
! 115: }
! 116: mpz_setbit (x, bit);
! 117: mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2);
! 118:
! 119: #if DEBUG
! 120: /* Check that the starting approximation is >= than the root. */
! 121: mpz_pow_ui (t1, x, nth);
! 122: if (mpz_cmp (cc, t1) >= 0)
! 123: abort ();
! 124: #endif
! 125:
! 126: mpz_add_ui (x, x, 1);
! 127:
! 128: /* Main loop */
! 129: do
! 130: {
! 131: lowz = mpz_scan1 (x, 0);
! 132: mpz_tdiv_q_2exp (t0, x, lowz);
! 133: mpz_pow_ui (t1, t0, nth - 1);
! 134: mpz_mul_2exp (t1, t1, lowz * (nth - 1));
! 135: mpz_tdiv_q (t2, cc, t1);
! 136: mpz_sub (t2, x, t2);
! 137: rl = mpz_tdiv_q_ui (t2, t2, nth);
! 138: mpz_sub (x, x, t2);
! 139: }
! 140: while (mpz_sgn (t2) != 0);
! 141:
! 142: /* If we got a non-zero remainder in the last division, we know our root
! 143: is too large. */
! 144: mpz_sub_ui (x, x, (mp_limb_t) (rl != 0));
! 145:
! 146: /* Adjustment loop. If we spend more care on rounding in the loop above,
! 147: we could probably get rid of this, or greatly simplify it. */
! 148: {
! 149: int bad = 0;
! 150: lowz = mpz_scan1 (x, 0);
! 151: mpz_tdiv_q_2exp (t0, x, lowz);
! 152: mpz_pow_ui (t1, t0, nth);
! 153: mpz_mul_2exp (t1, t1, lowz * nth);
! 154: while (mpz_cmp (cc, t1) < 0)
! 155: {
! 156: bad++;
! 157: if (bad > 2)
! 158: abort (); /* abort if our root is far off */
! 159: mpz_sub_ui (x, x, 1);
! 160: lowz = mpz_scan1 (x, 0);
! 161: mpz_tdiv_q_2exp (t0, x, lowz);
! 162: mpz_pow_ui (t1, t0, nth);
! 163: mpz_mul_2exp (t1, t1, lowz * nth);
! 164: }
! 165: }
! 166:
! 167: done:
! 168: exact = mpz_cmp (t1, cc) == 0;
! 169:
! 170: if (r != NULL)
! 171: {
! 172: mpz_set (r, x);
! 173: if (mpz_sgn (c) < 0)
! 174: SIZ(r) = -SIZ(r);
! 175: }
! 176:
! 177: mpz_clear (t2);
! 178: mpz_clear (t1);
! 179: mpz_clear (t0);
! 180: mpz_clear (x);
! 181:
! 182: return exact;
! 183: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to root.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz