version 1.1.1.1, 2000/09/09 14:12:56 |
version 1.1.1.2, 2003/08/25 16:06:33 |
|
|
/* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)). |
/* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)). |
Return an indication if the result is exact. |
Return an indication if the result is exact. |
|
|
Copyright (C) 1999, 2000 Free Software Foundation, Inc. |
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
Line 26 MA 02111-1307, USA. */ |
|
Line 26 MA 02111-1307, USA. */ |
|
to some extent. It would be natural to avoid representing the low zero |
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. */ |
bits mpz_scan1 is counting, and at the same time call mpn directly. */ |
|
|
#include <stdio.h> /* for NULL */ |
#include <stdio.h> /* for NULL */ |
|
#include <stdlib.h> /* for abort */ |
#include "gmp.h" |
#include "gmp.h" |
#include "gmp-impl.h" |
#include "gmp-impl.h" |
#include "longlong.h" |
#include "longlong.h" |
|
|
int |
int |
#if __STDC__ |
mpz_root (mpz_ptr r, mpz_srcptr u, unsigned long int nth) |
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; |
mp_ptr rootp, up; |
__mpz_struct ccs, *cc = &ccs; |
mp_size_t us, un, rootn; |
unsigned long int nbits; |
|
int bit; |
|
int exact; |
int exact; |
int i; |
unsigned int cnt; |
unsigned long int lowz; |
unsigned long int rootnb, unb; |
unsigned long int rl; |
|
|
|
|
up = PTR(u); |
|
us = SIZ(u); |
|
|
/* even roots of negatives provoke an exception */ |
/* even roots of negatives provoke an exception */ |
if (mpz_sgn (c) < 0 && (nth & 1) == 0) |
if (us < 0 && (nth & 1) == 0) |
SQRT_OF_NEGATIVE; |
SQRT_OF_NEGATIVE; |
|
|
/* root extraction interpreted as c^(1/nth) means a zeroth root should |
/* root extraction interpreted as c^(1/nth) means a zeroth root should |
Line 59 mpz_root (r, c, nth) |
|
Line 53 mpz_root (r, c, nth) |
|
if (nth == 0) |
if (nth == 0) |
DIVIDE_BY_ZERO; |
DIVIDE_BY_ZERO; |
|
|
if (mpz_sgn (c) == 0) |
if (us == 0) |
{ |
{ |
if (r != NULL) |
if (r != NULL) |
mpz_set_ui (r, 0); |
SIZ(r) = 0; |
return 1; /* exact result */ |
return 1; /* exact result */ |
} |
} |
|
|
PTR(cc) = PTR(c); |
un = ABS (us); |
SIZ(cc) = ABSIZ(c); |
|
|
|
nbits = (mpz_sizeinbase (cc, 2) - 1) / nth; |
count_leading_zeros (cnt, up[un - 1]); |
if (nbits == 0) |
unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS; |
|
rootnb = (unb - 1) / nth + 1; |
|
rootn = (rootnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; |
|
|
|
if (r != NULL) |
{ |
{ |
if (r != NULL) |
rootp = MPZ_REALLOC (r, rootn); |
mpz_set_ui (r, 1); |
up = PTR(u); |
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; |
|
} |
} |
|
else |
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); |
rootp = __GMP_ALLOCATE_FUNC_LIMBS (rootn); |
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 |
if (nth == 1) |
/* 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); |
MPN_COPY (rootp, up, un); |
mpz_tdiv_q_2exp (t0, x, lowz); |
exact = 1; |
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); |
else |
|
|
/* 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); |
exact = 0 == mpn_rootrem (rootp, NULL, up, un, nth); |
if (mpz_sgn (c) < 0) |
|
SIZ(r) = -SIZ(r); |
|
} |
} |
|
|
mpz_clear (t2); |
if (r != NULL) |
mpz_clear (t1); |
SIZ(r) = us >= 0 ? rootn : -rootn; |
mpz_clear (t0); |
else |
mpz_clear (x); |
__GMP_FREE_FUNC_LIMBS (rootp, rootn); |
|
|
return exact; |
return exact; |
} |
} |