rts/gmp/mpz/root.c

   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 }