Annotation of OpenXM_contrib/gmp/mpz/root.c, Revision 1.1.1.2
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:
1.1.1.2 ! ohara 4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 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:
1.1.1.2 ! ohara 29: #include <stdio.h> /* for NULL */
! 30: #include <stdlib.h> /* for abort */
1.1 maekawa 31: #include "gmp.h"
32: #include "gmp-impl.h"
33: #include "longlong.h"
34:
35: int
1.1.1.2 ! ohara 36: mpz_root (mpz_ptr r, mpz_srcptr u, unsigned long int nth)
1.1 maekawa 37: {
1.1.1.2 ! ohara 38: mp_ptr rootp, up;
! 39: mp_size_t us, un, rootn;
1.1 maekawa 40: int exact;
1.1.1.2 ! ohara 41: unsigned int cnt;
! 42: unsigned long int rootnb, unb;
! 43:
! 44: up = PTR(u);
! 45: us = SIZ(u);
1.1 maekawa 46:
47: /* even roots of negatives provoke an exception */
1.1.1.2 ! ohara 48: if (us < 0 && (nth & 1) == 0)
1.1 maekawa 49: SQRT_OF_NEGATIVE;
50:
51: /* root extraction interpreted as c^(1/nth) means a zeroth root should
52: provoke a divide by zero, do this even if c==0 */
53: if (nth == 0)
54: DIVIDE_BY_ZERO;
55:
1.1.1.2 ! ohara 56: if (us == 0)
1.1 maekawa 57: {
58: if (r != NULL)
1.1.1.2 ! ohara 59: SIZ(r) = 0;
1.1 maekawa 60: return 1; /* exact result */
61: }
62:
1.1.1.2 ! ohara 63: un = ABS (us);
! 64:
! 65: count_leading_zeros (cnt, up[un - 1]);
! 66: unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS;
! 67: rootnb = (unb - 1) / nth + 1;
! 68: rootn = (rootnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
1.1 maekawa 69:
1.1.1.2 ! ohara 70: if (r != NULL)
1.1 maekawa 71: {
1.1.1.2 ! ohara 72: rootp = MPZ_REALLOC (r, rootn);
! 73: up = PTR(u);
1.1 maekawa 74: }
1.1.1.2 ! ohara 75: else
1.1 maekawa 76: {
1.1.1.2 ! ohara 77: rootp = __GMP_ALLOCATE_FUNC_LIMBS (rootn);
1.1 maekawa 78: }
79:
1.1.1.2 ! ohara 80: if (nth == 1)
1.1 maekawa 81: {
1.1.1.2 ! ohara 82: MPN_COPY (rootp, up, un);
! 83: exact = 1;
1.1 maekawa 84: }
1.1.1.2 ! ohara 85: else
1.1 maekawa 86: {
1.1.1.2 ! ohara 87: exact = 0 == mpn_rootrem (rootp, NULL, up, un, nth);
1.1 maekawa 88: }
89:
1.1.1.2 ! ohara 90: if (r != NULL)
! 91: SIZ(r) = us >= 0 ? rootn : -rootn;
! 92: else
! 93: __GMP_FREE_FUNC_LIMBS (rootp, rootn);
1.1 maekawa 94:
95: return exact;
96: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>