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

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

version 1.1.1.1, 2000/09/09 14:12:51 version 1.1.1.2, 2003/08/25 16:06:33
Line 1 
Line 1 
 /* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number.  /* mpz_fib_ui -- calculate Fibonacci numbers.
   
 Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.  Copyright 2000, 2001 Free Software Foundation, Inc.
   
 This file is part of the GNU MP Library.  This file is part of the GNU MP Library.
   
Line 19  along with the GNU MP Library; see the file COPYING.LI
Line 19  along with the GNU MP Library; see the file COPYING.LI
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */  MA 02111-1307, USA. */
   
   #include <stdio.h>
 #include "gmp.h"  #include "gmp.h"
 #include "gmp-impl.h"  #include "gmp-impl.h"
   #include "longlong.h"
   
 /* This is fast, but could be made somewhat faster and neater.  
    The timing is somewhat fluctuating for even/odd sizes because  
    of the extra hair used to save variables and operations.  Here  
    are a few things one might want to address:  
      1. Avoid using 4 intermediate variables in mpz_fib_bigcase.  
      2. Call mpn functions directly.  Straightforward for these functions.  
      3. Merge the three functions into one.  
   
 Said by Kevin:  /* change to "#define TRACE(x) x" to get some traces */
    Consider using the Lucas numbers L[n] as an auxiliary sequence, making  #define TRACE(x)
    it possible to do the "doubling" operation in mpz_fib_bigcase with two  
    squares rather than two multiplies.  The formulas are a little more  
    complicated, something like the following (untested).  
   
        F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2  
        L[2n] = 5*F[n]^2 + 2*(-1)^n  
   
        F[2n+1] = (F[2n] + L[2n]) / 2  /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
        L[2n+1] = (5*F[2n] + L[2n]) / 2     limb because the result F[2k+1] is an F[4m+3] and such numbers are always
      == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7.  (This is
      the same as in mpn_fib2_ui.)
   
    The Lucas number that comes for free here could even be returned.     In the F[2k+1] for k even, the +2 won't give a carry out of the low limb
      in normal circumstances.  This is an F[4m+1] and we claim that F[3*2^b+1]
      == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence
      if n < 2^BITS_PER_MP_LIMB then F[n] cannot have a low limb of 0 or 1.  No
      proof for this claim, but it's been verified up to b==32 and has such a
      nice pattern it must be true :-).  Of interest is that F[3*2^b] == 0 mod
      2^(b+1) seems to hold too.
   
    Maybe there's formulas with two squares using just F[n], but I don't     When n >= 2^BITS_PER_MP_LIMB, which can arise in test setups with a small
    know of any.     limb, then the low limb of F[4m+1] can certainly be 1, and an mpn_add_1
 */     must be used.  */
   
 /* Determine the needed storage for Fib(n).  */  
 #define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2)  
   
 static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int));  
 static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int));  
   
   
 #ifndef FIB_THRESHOLD  
 #define FIB_THRESHOLD 60  
 #endif  
   
 void  void
 #if __STDC__  mpz_fib_ui (mpz_ptr fn, unsigned long n)
 mpz_fib_ui (mpz_t r, unsigned long int n)  
 #else  
 mpz_fib_ui (r, n)  
      mpz_t r;  
      unsigned long int n;  
 #endif  
 {  {
   if (n == 0)    mp_ptr         fp, xp, yp;
     mpz_set_ui (r, 0);    mp_size_t      size, xalloc;
   else    unsigned long  n2;
     mp_limb_t      c, c2;
     TMP_DECL (marker);
   
     if (n <= FIB_TABLE_LIMIT)
     {      {
       mpz_t t1;        PTR(fn)[0] = FIB_TABLE (n);
       mpz_init (t1);        SIZ(fn) = (n != 0);      /* F[0]==0, others are !=0 */
       if (n < FIB_THRESHOLD)        return;
         mpz_fib_basecase (t1, r, n);  
       else  
         mpz_fib_bigcase (t1, r, n);  
       mpz_clear (t1);  
     }      }
 }  
   
 static void    n2 = n/2;
 #if __STDC__    xalloc = MPN_FIB2_SIZE (n2) + 1;
 mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n)    MPZ_REALLOC (fn, 2*xalloc+1);
 #else    fp = PTR (fn);
 mpz_fib_basecase (t1, t2, n)  
      mpz_t t1;  
      mpz_t t2;  
      unsigned long int n;  
 #endif  
 {  
   unsigned long int m, i;  
   
   mpz_set_ui (t1, 0);    TMP_MARK (marker);
   mpz_set_ui (t2, 1);    TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
   m = n/2;    size = mpn_fib2_ui (xp, yp, n2);
   for (i = 0; i < m; i++)  
     TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n",
                    n >> 1, size, n&1);
            mpn_trace ("xp", xp, size);
            mpn_trace ("yp", yp, size));
   
     if (n & 1)
     {      {
       mpz_add (t1, t1, t2);        /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k  */
       mpz_add (t2, t1, t2);        mp_size_t  xsize, ysize;
     }  
   if ((n & 1) == 0)  
     {  
       mpz_sub (t1, t2, t1);  
       mpz_sub (t2, t2, t1);     /* trick: recover t1 value just overwritten */  
     }  
 }  
   
 static void  #if HAVE_NATIVE_mpn_addsub_n
 #if __STDC__        xp[size] = mpn_lshift (xp, xp, size, 1);
 mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n)        yp[size] = 0;
         ASSERT_NOCARRY (mpn_addsub_n (xp, yp, xp, yp, size+1));
         xsize = size + (xp[size] != 0);
         ysize = size + (yp[size] != 0);
 #else  #else
 mpz_fib_bigcase (t1, t2, n)        c2 = mpn_lshift (fp, xp, size, 1);
      mpz_t t1;        c = c2 + mpn_add_n (xp, fp, yp, size);
      mpz_t t2;        xp[size] = c;
      unsigned long int n;        xsize = size + (c != 0);
         c2 -= mpn_sub_n (yp, fp, yp, size);
         yp[size] = c2;
         ASSERT (c2 <= 1);
         ysize = size + c2;
 #endif  #endif
 {  
   unsigned long int n2;  
   int ni, i;  
   mpz_t x1, x2, u1, u2;  
   
   ni = 0;        size = xsize + ysize;
   for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2)        c = mpn_mul (fp, xp, xsize, yp, ysize);
     ni++;  
   
   mpz_fib_basecase (t1, t2, n2);  #if BITS_PER_MP_LIMB >= BITS_PER_ULONG
         /* no overflow, see comments above */
   mpz_init (x1);        ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= MP_LIMB_T_MAX-2);
   mpz_init (x2);        fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
   mpz_init (u1);  #else
   mpz_init (u2);        /* this code only for testing with small limbs, limb<ulong is unusual */
         if (n & 2)
   for (i = ni - 1; i >= 0; i--)          {
             ASSERT (fp[0] >= 2);
             fp[0] -= 2;
           }
         else
           {
             ASSERT (c != MP_LIMB_T_MAX); /* because it's the high of a mul */
             c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2));
             fp[size-1] = c;
           }
   #endif
       }
     else
     {      {
       mpz_mul_2exp (x1, t1, 1);        /* F[2k] = F[k]*(F[k]+2F[k-1]) */
       mpz_mul_2exp (x2, t2, 1);  
   
       mpz_add (x1, x1, t2);        mp_size_t  xsize, ysize;
       mpz_sub (x2, x2, t1);        c = mpn_lshift (yp, yp, size, 1);
         c += mpn_add_n (yp, yp, xp, size);
         yp[size] = c;
         xsize = size;
         ysize = size + (c != 0);
         size += ysize;
         c = mpn_mul (fp, yp, ysize, xp, xsize);
       }
   
       mpz_mul (u1, t2, x1);    /* one or two high zeros */
       mpz_mul (u2, t1, x2);    size -= (c == 0);
     size -= (fp[size-1] == 0);
     SIZ(fn) = size;
   
       if (((n >> i) & 1) == 0)    TRACE (printf ("done special, size=%ld\n", size);
         {           mpn_trace ("fp ", fp, size));
           mpz_sub (t1, u1, u2);  
           mpz_set (t2, u1);  
         }  
       else  
         {  
           mpz_set (t1, u1);  
           mpz_mul_2exp (t2, u1, 1);  
           mpz_sub (t2, t2, u2);  
         }  
     }  
   
   mpz_clear (x1);    TMP_FREE (marker);
   mpz_clear (x2);  
   mpz_clear (u1);  
   mpz_clear (u2);  
 }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>