From dac7b0aca943ac7db905e3e97ebe8ca27e7b0a71 Mon Sep 17 00:00:00 2001 From: tege Date: Wed, 8 May 2002 21:00:44 +0200 Subject: *** empty log message *** --- ChangeLog | 10 ++++ mpn/generic/rootrem.c | 158 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+) create mode 100644 mpn/generic/rootrem.c diff --git a/ChangeLog b/ChangeLog index f0e203994..7e083ddc9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -19,6 +19,16 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +2002-05-08 Torbjorn Granlund + + * mpn/generic/rootrem.c: New file. + + * mpz/perfpow.c: Amend allocations for mpn_rootrem requirements. + + * gmp-impl.h (mpn_rootrem): Add declaration. + + * mpz/root.c: Rewrite to use mpn_rootrem. + 2002-05-08 Kevin Ryde * gmp-impl.h (MUL_KARATSUBA_THRESHOLD etc): Remove forced nail values. diff --git a/mpn/generic/rootrem.c b/mpn/generic/rootrem.c new file mode 100644 index 000000000..7c2e50b0d --- /dev/null +++ b/mpn/generic/rootrem.c @@ -0,0 +1,158 @@ +/* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and + store the truncated integer part at rootp and the remainder at remp. + + THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE + INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. + IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A + FUTURE GNU MP RELEASE. + + +Copyright 2002 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. */ + +/* + We use Newton's method to compute the root of a: + + n + f(x) := x - a + + + n - 1 + f'(x) := x n + + + n-1 n-1 n-1 + x - a/x a/x - x a/x + (n-1)x + new x = x - f(x)/f'(x) = x - ---------- = x + --------- = -------------- + n n n + +*/ + + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +mp_size_t +mpn_rootrem (mp_ptr rootp, mp_ptr remp, + mp_srcptr up, mp_size_t un, mp_limb_t nth) +{ + mp_ptr pp, qp; + mp_size_t pn, xn, qn; + unsigned long int unb, xnb, bit; + unsigned int cnt; + mp_size_t i; + unsigned long int n_valid_bits, adj; + TMP_DECL (marker); + + TMP_MARK (marker); + pp = TMP_ALLOC_LIMBS (un + 2); + qp = TMP_ALLOC_LIMBS (un); /* should try and reduce this allocation */ + + count_leading_zeros (cnt, up[un - 1]); + unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS; + + xnb = (unb + nth - 1) / nth; + if (xnb == 1) + { + rootp[0] = 1; + if (remp == NULL) + remp = pp; + mpn_sub_1 (remp, up, un, (mp_limb_t) 1); + MPN_NORMALIZE (remp, un); + TMP_FREE (marker); + return un; + } + + xn = (xnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; + + /* Set initial root to only ones. This is an overestimate of the actual root + by less than a factor of 2. */ + for (i = 0; i < xn; i++) + rootp[i] = GMP_NUMB_MAX; + rootp[xnb / GMP_NUMB_BITS] = ((mp_limb_t) 1 << (xnb % GMP_NUMB_BITS)) - 1; + + /* Improve the initial approximation, one bit at a time. Keep the + approximations >= root(U,nth). */ + bit = xnb - 2; + n_valid_bits = 0; + for (i = 0; (nth >> i) != 0; i++) + { + mp_limb_t xl = rootp[bit / GMP_NUMB_BITS]; + rootp[bit / GMP_NUMB_BITS] = xl ^ (mp_limb_t) 1 << bit % GMP_NUMB_BITS; + pn = mpn_pow_1 (pp, rootp, xn, nth, qp); + /* If the new root approximation is too small, restore old value. */ + if (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0))) + rootp[bit / GMP_NUMB_BITS] = xl; /* restore old value */ + n_valid_bits += 1; + if (bit == 0) + goto done; + bit--; + } + + adj = n_valid_bits - 1; + + /* Newton loop. Converges downwards towards root(U,nth). Currently we use + full precision from iteration 1. Clearly, we should use just n_valid_bits + of precision in each step, and thus save most of the computations. */ + while (n_valid_bits <= xnb) + { + mp_limb_t cy; + + pn = mpn_pow_1 (pp, rootp, xn, nth - 1, qp); + mpn_tdiv_qr (qp, pp, (mp_size_t) 0, up, un, pp, pn); /* junk remainder */ + qn = un - pn + 1; + cy = mpn_addmul_1 (qp, rootp, xn, nth - 1); + if (qn == xn + 1) + { + cy += qp[xn]; + if (cy == nth) + { + for (i = xn - 1; i >= 0; i--) + qp[i] = GMP_NUMB_MAX; + cy = nth - 1; + } + } + + qp[xn] = cy; + qn = xn + (cy != 0); + + mpn_divrem_1 (rootp, (mp_size_t) 0, qp, qn, nth); + n_valid_bits = n_valid_bits * 2 - adj; + } + + /* The computed result might be one unit too large. Adjust as necessary. */ + done: + pn = mpn_pow_1 (pp, rootp, xn, nth, qp); + if (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)) + { + mpn_decr_u (rootp, 1); + pn = mpn_pow_1 (pp, rootp, xn, nth, qp); + + ASSERT_ALWAYS (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0))); + } + + if (remp == NULL) + remp = pp; + mpn_sub (remp, up, un, pp, pn); + MPN_NORMALIZE (remp, un); + + TMP_FREE (marker); + return un; +} -- cgit v1.2.1