Annotation of OpenXM_contrib/gmp/mpz/root.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)).
2: Return an indication if the result is exact.
3:
4: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
5:
6: This file is part of the GNU MP Library.
7:
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.
12:
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.
17:
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. */
22:
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. */
28:
29: #include <stdio.h> /* for NULL */
30: #include "gmp.h"
31: #include "gmp-impl.h"
32: #include "longlong.h"
33:
34: int
35: #if __STDC__
36: mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth)
37: #else
38: mpz_root (r, c, nth)
39: mpz_ptr r;
40: mpz_srcptr c;
41: unsigned long int nth;
42: #endif
43: {
44: mpz_t x, t0, t1, t2;
45: __mpz_struct ccs, *cc = &ccs;
46: unsigned long int nbits;
47: int bit;
48: int exact;
49: int i;
50: unsigned long int lowz;
51: unsigned long int rl;
52:
53: /* even roots of negatives provoke an exception */
54: if (mpz_sgn (c) < 0 && (nth & 1) == 0)
55: SQRT_OF_NEGATIVE;
56:
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 */
59: if (nth == 0)
60: DIVIDE_BY_ZERO;
61:
62: if (mpz_sgn (c) == 0)
63: {
64: if (r != NULL)
65: mpz_set_ui (r, 0);
66: return 1; /* exact result */
67: }
68:
69: PTR(cc) = PTR(c);
70: SIZ(cc) = ABSIZ(c);
71:
72: nbits = (mpz_sizeinbase (cc, 2) - 1) / nth;
73: if (nbits == 0)
74: {
75: if (r != NULL)
76: mpz_set_ui (r, 1);
77: if (mpz_sgn (c) < 0)
78: {
79: if (r != NULL)
80: SIZ(r) = -SIZ(r);
81: return mpz_cmp_si (c, -1L) == 0;
82: }
83: return mpz_cmp_ui (c, 1L) == 0;
84: }
85:
86: mpz_init (x);
87: mpz_init (t0);
88: mpz_init (t1);
89: mpz_init (t2);
90:
91: /* Create a one-bit approximation. */
92: mpz_set_ui (x, 0);
93: mpz_setbit (x, nbits);
94:
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. */
98: bit = nbits - 1;
99: for (i = 1; (nth >> i) != 0; i++)
100: {
101: mpz_setbit (x, bit);
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)
106: mpz_clrbit (x, bit);
107:
108: bit--; /* check/set next bit */
109: if (bit < 0)
110: {
111: /* We're done. */
112: mpz_pow_ui (t1, x, nth);
113: goto done;
114: }
115: }
116: mpz_setbit (x, bit);
117: mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2);
118:
119: #if DEBUG
120: /* Check that the starting approximation is >= than the root. */
121: mpz_pow_ui (t1, x, nth);
122: if (mpz_cmp (cc, t1) >= 0)
123: abort ();
124: #endif
125:
126: mpz_add_ui (x, x, 1);
127:
128: /* Main loop */
129: do
130: {
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);
136: mpz_sub (t2, x, t2);
137: rl = mpz_tdiv_q_ui (t2, t2, nth);
138: mpz_sub (x, x, t2);
139: }
140: while (mpz_sgn (t2) != 0);
141:
142: /* If we got a non-zero remainder in the last division, we know our root
143: is too large. */
144: mpz_sub_ui (x, x, (mp_limb_t) (rl != 0));
145:
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. */
148: {
149: int bad = 0;
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)
155: {
156: bad++;
157: if (bad > 2)
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);
164: }
165: }
166:
167: done:
168: exact = mpz_cmp (t1, cc) == 0;
169:
170: if (r != NULL)
171: {
172: mpz_set (r, x);
173: if (mpz_sgn (c) < 0)
174: SIZ(r) = -SIZ(r);
175: }
176:
177: mpz_clear (t2);
178: mpz_clear (t1);
179: mpz_clear (t0);
180: mpz_clear (x);
181:
182: return exact;
183: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to root.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz