Annotation of OpenXM_contrib/gmp/mpn/generic/rootrem.c, Revision 1.1
1.1 ! ohara 1: /* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and
! 2: store the truncated integer part at rootp and the remainder at remp.
! 3:
! 4: THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
! 5: INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
! 6: IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
! 7: FUTURE GNU MP RELEASE.
! 8:
! 9:
! 10: Copyright 2002 Free Software Foundation, Inc.
! 11:
! 12: This file is part of the GNU MP Library.
! 13:
! 14: The GNU MP Library is free software; you can redistribute it and/or modify
! 15: it under the terms of the GNU Lesser General Public License as published by
! 16: the Free Software Foundation; either version 2.1 of the License, or (at your
! 17: option) any later version.
! 18:
! 19: The GNU MP Library is distributed in the hope that it will be useful, but
! 20: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 21: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 22: License for more details.
! 23:
! 24: You should have received a copy of the GNU Lesser General Public License
! 25: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 26: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 27: MA 02111-1307, USA. */
! 28:
! 29: /*
! 30: We use Newton's method to compute the root of a:
! 31:
! 32: n
! 33: f(x) := x - a
! 34:
! 35:
! 36: n - 1
! 37: f'(x) := x n
! 38:
! 39:
! 40: n-1 n-1 n-1
! 41: x - a/x a/x - x a/x + (n-1)x
! 42: new x = x - f(x)/f'(x) = x - ---------- = x + --------- = --------------
! 43: n n n
! 44:
! 45: */
! 46:
! 47:
! 48: #include <stdio.h>
! 49: #include <stdlib.h>
! 50:
! 51: #include "gmp.h"
! 52: #include "gmp-impl.h"
! 53: #include "longlong.h"
! 54:
! 55: mp_size_t
! 56: mpn_rootrem (mp_ptr rootp, mp_ptr remp,
! 57: mp_srcptr up, mp_size_t un, mp_limb_t nth)
! 58: {
! 59: mp_ptr pp, qp, xp;
! 60: mp_size_t pn, xn, qn;
! 61: unsigned long int unb, xnb, bit;
! 62: unsigned int cnt;
! 63: mp_size_t i;
! 64: unsigned long int n_valid_bits, adj;
! 65: TMP_DECL (marker);
! 66:
! 67: TMP_MARK (marker);
! 68:
! 69: /* The extra factor 1.585 = log(3)/log(2) here is for the worst case
! 70: overestimate of the root, i.e., when the code rounds a root that is
! 71: 2+epsilon to 3, and the powers this to a potentially huge power. We
! 72: could generalize the code for detecting root=1 a few lines below to deal
! 73: with xnb <= k, for some small k. For example, when xnb <= 2, meaning
! 74: the root should be 1, 2, or 3, we could replace this factor by the much
! 75: smaller log(5)/log(4). */
! 76:
! 77: #define PP_ALLOC (2 + (mp_size_t) (un*1.585))
! 78: pp = TMP_ALLOC_LIMBS (PP_ALLOC);
! 79:
! 80: count_leading_zeros (cnt, up[un - 1]);
! 81: unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS;
! 82:
! 83: xnb = (unb - 1) / nth + 1;
! 84: if (xnb == 1)
! 85: {
! 86: if (remp == NULL)
! 87: remp = pp;
! 88: mpn_sub_1 (remp, up, un, (mp_limb_t) 1);
! 89: MPN_NORMALIZE (remp, un);
! 90: rootp[0] = 1;
! 91: TMP_FREE (marker);
! 92: return un;
! 93: }
! 94:
! 95: xn = (xnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
! 96:
! 97: qp = TMP_ALLOC_LIMBS (un + 1);
! 98: xp = TMP_ALLOC_LIMBS (xn + 1);
! 99:
! 100: /* Set initial root to only ones. This is an overestimate of the actual root
! 101: by less than a factor of 2. */
! 102: for (i = 0; i < xn; i++)
! 103: xp[i] = GMP_NUMB_MAX;
! 104: xp[xnb / GMP_NUMB_BITS] = ((mp_limb_t) 1 << (xnb % GMP_NUMB_BITS)) - 1;
! 105:
! 106: /* Improve the initial approximation, one bit at a time. Keep the
! 107: approximations >= root(U,nth). */
! 108: bit = xnb - 2;
! 109: n_valid_bits = 0;
! 110: for (i = 0; (nth >> i) != 0; i++)
! 111: {
! 112: mp_limb_t xl = xp[bit / GMP_NUMB_BITS];
! 113: xp[bit / GMP_NUMB_BITS] = xl ^ (mp_limb_t) 1 << bit % GMP_NUMB_BITS;
! 114: pn = mpn_pow_1 (pp, xp, xn, nth, qp);
! 115: ASSERT_ALWAYS (pn < PP_ALLOC);
! 116: /* If the new root approximation is too small, restore old value. */
! 117: if (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)))
! 118: xp[bit / GMP_NUMB_BITS] = xl; /* restore old value */
! 119: n_valid_bits += 1;
! 120: if (bit == 0)
! 121: goto done;
! 122: bit--;
! 123: }
! 124:
! 125: adj = n_valid_bits - 1;
! 126:
! 127: /* Newton loop. Converges downwards towards root(U,nth). Currently we use
! 128: full precision from iteration 1. Clearly, we should use just n_valid_bits
! 129: of precision in each step, and thus save most of the computations. */
! 130: while (n_valid_bits <= xnb)
! 131: {
! 132: mp_limb_t cy;
! 133:
! 134: pn = mpn_pow_1 (pp, xp, xn, nth - 1, qp);
! 135: ASSERT_ALWAYS (pn < PP_ALLOC);
! 136: qp[xn - 1] = 0; /* pad quotient to make it always xn limbs */
! 137: mpn_tdiv_qr (qp, pp, (mp_size_t) 0, up, un, pp, pn); /* junk remainder */
! 138: cy = mpn_addmul_1 (qp, xp, xn, nth - 1);
! 139: if (un - pn == xn)
! 140: {
! 141: cy += qp[xn];
! 142: if (cy == nth)
! 143: {
! 144: for (i = xn - 1; i >= 0; i--)
! 145: qp[i] = GMP_NUMB_MAX;
! 146: cy = nth - 1;
! 147: }
! 148: }
! 149:
! 150: qp[xn] = cy;
! 151: qn = xn + (cy != 0);
! 152:
! 153: mpn_divrem_1 (xp, (mp_size_t) 0, qp, qn, nth);
! 154: n_valid_bits = n_valid_bits * 2 - adj;
! 155: }
! 156:
! 157: /* The computed result might be one unit too large. Adjust as necessary. */
! 158: done:
! 159: pn = mpn_pow_1 (pp, xp, xn, nth, qp);
! 160: ASSERT_ALWAYS (pn < PP_ALLOC);
! 161: if (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0))
! 162: {
! 163: mpn_decr_u (xp, 1);
! 164: pn = mpn_pow_1 (pp, xp, xn, nth, qp);
! 165: ASSERT_ALWAYS (pn < PP_ALLOC);
! 166:
! 167: ASSERT_ALWAYS (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)));
! 168: }
! 169:
! 170: if (remp == NULL)
! 171: remp = pp;
! 172: mpn_sub (remp, up, un, pp, pn);
! 173: MPN_NORMALIZE (remp, un);
! 174: MPN_COPY (rootp, xp, xn);
! 175: TMP_FREE (marker);
! 176: return un;
! 177: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rootrem.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic