[BACK]Return to root.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpz

Diff for /OpenXM_contrib/gmp/mpz/Attic/root.c between version 1.1.1.1 and 1.1.1.2

version 1.1.1.1, 2000/09/09 14:12:56 version 1.1.1.2, 2003/08/25 16:06:33
Line 1 
Line 1 
 /* mpz_root(root, u, nth) --  Set ROOT to floor(U^(1/nth)).  /* mpz_root(root, u, nth) --  Set ROOT to floor(U^(1/nth)).
    Return an indication if the result is exact.     Return an indication if the result is exact.
   
 Copyright (C) 1999, 2000 Free Software Foundation, Inc.  Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
   
 This file is part of the GNU MP Library.  This file is part of the GNU MP Library.
   
Line 26  MA 02111-1307, USA. */
Line 26  MA 02111-1307, USA. */
    to some extent.  It would be natural to avoid representing the low zero     to some extent.  It would be natural to avoid representing the low zero
    bits mpz_scan1 is counting, and at the same time call mpn directly.  */     bits mpz_scan1 is counting, and at the same time call mpn directly.  */
   
 #include <stdio.h> /* for NULL */  #include <stdio.h>  /* for NULL */
   #include <stdlib.h> /* for abort */
 #include "gmp.h"  #include "gmp.h"
 #include "gmp-impl.h"  #include "gmp-impl.h"
 #include "longlong.h"  #include "longlong.h"
   
 int  int
 #if __STDC__  mpz_root (mpz_ptr r, mpz_srcptr u, unsigned long int nth)
 mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth)  
 #else  
 mpz_root (r, c, nth)  
      mpz_ptr r;  
      mpz_srcptr c;  
      unsigned long int nth;  
 #endif  
 {  {
   mpz_t x, t0, t1, t2;    mp_ptr rootp, up;
   __mpz_struct ccs, *cc = &ccs;    mp_size_t us, un, rootn;
   unsigned long int nbits;  
   int bit;  
   int exact;    int exact;
   int i;    unsigned int cnt;
   unsigned long int lowz;    unsigned long int rootnb, unb;
   unsigned long int rl;  
   
     up = PTR(u);
     us = SIZ(u);
   
   /* even roots of negatives provoke an exception */    /* even roots of negatives provoke an exception */
   if (mpz_sgn (c) < 0 && (nth & 1) == 0)    if (us < 0 && (nth & 1) == 0)
     SQRT_OF_NEGATIVE;      SQRT_OF_NEGATIVE;
   
   /* root extraction interpreted as c^(1/nth) means a zeroth root should    /* root extraction interpreted as c^(1/nth) means a zeroth root should
Line 59  mpz_root (r, c, nth)
Line 53  mpz_root (r, c, nth)
   if (nth == 0)    if (nth == 0)
     DIVIDE_BY_ZERO;      DIVIDE_BY_ZERO;
   
   if (mpz_sgn (c) == 0)    if (us == 0)
     {      {
       if (r != NULL)        if (r != NULL)
         mpz_set_ui (r, 0);          SIZ(r) = 0;
       return 1;                 /* exact result */        return 1;                 /* exact result */
     }      }
   
   PTR(cc) = PTR(c);    un = ABS (us);
   SIZ(cc) = ABSIZ(c);  
   
   nbits = (mpz_sizeinbase (cc, 2) - 1) / nth;    count_leading_zeros (cnt, up[un - 1]);
   if (nbits == 0)    unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS;
     rootnb = (unb - 1) / nth + 1;
     rootn = (rootnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
   
     if (r != NULL)
     {      {
       if (r != NULL)        rootp = MPZ_REALLOC (r, rootn);
         mpz_set_ui (r, 1);        up = PTR(u);
       if (mpz_sgn (c) < 0)  
         {  
           if (r != NULL)  
             SIZ(r) = -SIZ(r);  
           return mpz_cmp_si (c, -1L) == 0;  
         }  
       return mpz_cmp_ui (c, 1L) == 0;  
     }      }
     else
   mpz_init (x);  
   mpz_init (t0);  
   mpz_init (t1);  
   mpz_init (t2);  
   
   /* Create a one-bit approximation.  */  
   mpz_set_ui (x, 0);  
   mpz_setbit (x, nbits);  
   
   /* Make the approximation better, one bit at a time.  This odd-looking  
      termination criteria makes large nth get better initial approximation,  
      which avoids slow convergence for such values.  */  
   bit = nbits - 1;  
   for (i = 1; (nth >> i) != 0; i++)  
     {      {
       mpz_setbit (x, bit);        rootp = __GMP_ALLOCATE_FUNC_LIMBS (rootn);
       mpz_tdiv_q_2exp (t0, x, bit);  
       mpz_pow_ui (t1, t0, nth);  
       mpz_mul_2exp (t1, t1, bit * nth);  
       if (mpz_cmp (cc, t1) < 0)  
         mpz_clrbit (x, bit);  
   
       bit--;                    /* check/set next bit */  
       if (bit < 0)  
         {  
           /* We're done.  */  
           mpz_pow_ui (t1, x, nth);  
           goto done;  
         }  
     }      }
   mpz_setbit (x, bit);  
   mpz_set_ui (t2, 0); mpz_setbit (t2, bit);  mpz_add (x, x, t2);  
   
 #if DEBUG    if (nth == 1)
   /* Check that the starting approximation is >= than the root.  */  
   mpz_pow_ui (t1, x, nth);  
   if (mpz_cmp (cc, t1) >= 0)  
     abort ();  
 #endif  
   
   mpz_add_ui (x, x, 1);  
   
   /* Main loop */  
   do  
     {      {
       lowz = mpz_scan1 (x, 0);        MPN_COPY (rootp, up, un);
       mpz_tdiv_q_2exp (t0, x, lowz);        exact = 1;
       mpz_pow_ui (t1, t0, nth - 1);  
       mpz_mul_2exp (t1, t1, lowz * (nth - 1));  
       mpz_tdiv_q (t2, cc, t1);  
       mpz_sub (t2, x, t2);  
       rl = mpz_tdiv_q_ui (t2, t2, nth);  
       mpz_sub (x, x, t2);  
     }      }
   while (mpz_sgn (t2) != 0);    else
   
   /* If we got a non-zero remainder in the last division, we know our root  
      is too large.  */  
   mpz_sub_ui (x, x, (mp_limb_t) (rl != 0));  
   
   /* Adjustment loop.  If we spend more care on rounding in the loop above,  
      we could probably get rid of this, or greatly simplify it.  */  
   {  
     int bad = 0;  
     lowz = mpz_scan1 (x, 0);  
     mpz_tdiv_q_2exp (t0, x, lowz);  
     mpz_pow_ui (t1, t0, nth);  
     mpz_mul_2exp (t1, t1, lowz * nth);  
     while (mpz_cmp (cc, t1) < 0)  
       {  
         bad++;  
         if (bad > 2)  
           abort ();                     /* abort if our root is far off */  
         mpz_sub_ui (x, x, 1);  
         lowz = mpz_scan1 (x, 0);  
         mpz_tdiv_q_2exp (t0, x, lowz);  
         mpz_pow_ui (t1, t0, nth);  
         mpz_mul_2exp (t1, t1, lowz * nth);  
       }  
   }  
   
  done:  
   exact = mpz_cmp (t1, cc) == 0;  
   
   if (r != NULL)  
     {      {
       mpz_set (r, x);        exact = 0 == mpn_rootrem (rootp, NULL, up, un, nth);
       if (mpz_sgn (c) < 0)  
         SIZ(r) = -SIZ(r);  
     }      }
   
   mpz_clear (t2);    if (r != NULL)
   mpz_clear (t1);      SIZ(r) = us >= 0 ? rootn : -rootn;
   mpz_clear (t0);    else
   mpz_clear (x);      __GMP_FREE_FUNC_LIMBS (rootp, rootn);
   
   return exact;    return exact;
 }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>