1 /* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)).
2 Return an indication if the result is exact.
4 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
6 This file is part of the GNU MP Library.
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.
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.
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. */
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. */
29 #include <stdio.h> /* for NULL */
36 mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth)
41 unsigned long int nth;
45 __mpz_struct ccs, *cc = &ccs;
46 unsigned long int nbits;
50 unsigned long int lowz;
53 /* even roots of negatives provoke an exception */
54 if (mpz_sgn (c) < 0 && (nth & 1) == 0)
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 */
66 return 1; /* exact result */
72 nbits = (mpz_sizeinbase (cc, 2) - 1) / nth;
81 return mpz_cmp_si (c, -1L) == 0;
83 return mpz_cmp_ui (c, 1L) == 0;
91 /* Create a one-bit approximation. */
93 mpz_setbit (x, nbits);
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. */
99 for (i = 1; (nth >> i) != 0; i++)
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)
108 bit--; /* check/set next bit */
112 mpz_pow_ui (t1, x, nth);
117 mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2);
120 /* Check that the starting approximation is >= than the root. */
121 mpz_pow_ui (t1, x, nth);
122 if (mpz_cmp (cc, t1) >= 0)
126 mpz_add_ui (x, x, 1);
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);
137 rl = mpz_tdiv_q_ui (t2, t2, nth);
140 while (mpz_sgn (t2) != 0);
142 /* If we got a non-zero remainder in the last division, we know our root
144 mpz_sub_ui (x, x, (mp_limb_t) (rl != 0));
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. */
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)
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);
168 exact = mpz_cmp (t1, cc) == 0;