Return to root.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

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: }