1
/* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)).
2
Return an indication if the result is exact.
4
Copyright 1999, 2000, 2001 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 */
35
mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth)
38
__mpz_struct ccs, *cc = &ccs;
39
unsigned long int nbits;
43
unsigned long int lowz;
46
/* even roots of negatives provoke an exception */
47
if (mpz_sgn (c) < 0 && (nth & 1) == 0)
50
/* root extraction interpreted as c^(1/nth) means a zeroth root should
51
provoke a divide by zero, do this even if c==0 */
59
return 1; /* exact result */
65
nbits = (mpz_sizeinbase (cc, 2) - 1) / nth;
74
return mpz_cmp_si (c, -1L) == 0;
76
return mpz_cmp_ui (c, 1L) == 0;
84
/* Create a one-bit approximation. */
86
mpz_setbit (x, nbits);
88
/* Make the approximation better, one bit at a time. This odd-looking
89
termination criteria makes large nth get better initial approximation,
90
which avoids slow convergence for such values. */
92
for (i = 1; (nth >> i) != 0; i++)
95
mpz_tdiv_q_2exp (t0, x, bit);
96
mpz_pow_ui (t1, t0, nth);
97
mpz_mul_2exp (t1, t1, bit * nth);
98
if (mpz_cmp (cc, t1) < 0)
101
bit--; /* check/set next bit */
105
mpz_pow_ui (t1, x, nth);
110
mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2);
113
/* Check that the starting approximation is >= than the root. */
114
mpz_pow_ui (t1, x, nth);
115
if (mpz_cmp (cc, t1) >= 0)
119
mpz_add_ui (x, x, 1);
124
lowz = mpz_scan1 (x, 0);
125
mpz_tdiv_q_2exp (t0, x, lowz);
126
mpz_pow_ui (t1, t0, nth - 1);
127
mpz_mul_2exp (t1, t1, lowz * (nth - 1));
128
mpz_tdiv_q (t2, cc, t1);
130
rl = mpz_tdiv_q_ui (t2, t2, nth);
133
while (mpz_sgn (t2) != 0);
135
/* If we got a non-zero remainder in the last division, we know our root
137
mpz_sub_ui (x, x, (mp_limb_t) (rl != 0));
139
/* Adjustment loop. If we spend more care on rounding in the loop above,
140
we could probably get rid of this, or greatly simplify it. */
143
lowz = mpz_scan1 (x, 0);
144
mpz_tdiv_q_2exp (t0, x, lowz);
145
mpz_pow_ui (t1, t0, nth);
146
mpz_mul_2exp (t1, t1, lowz * nth);
147
while (mpz_cmp (cc, t1) < 0)
151
abort (); /* abort if our root is far off */
152
mpz_sub_ui (x, x, 1);
153
lowz = mpz_scan1 (x, 0);
154
mpz_tdiv_q_2exp (t0, x, lowz);
155
mpz_pow_ui (t1, t0, nth);
156
mpz_mul_2exp (t1, t1, lowz * nth);
161
exact = mpz_cmp (t1, cc) == 0;