'Is this code use the Fast Doubling method for Fibonacci number calculation?

I need to understand how the fib_ui.c function in gmp repo work :

in their documentation, they state that :

Beyond the table, values are generated with a binary powering algorithm, calculating a pair F[n] and F[n-1] working from high to low across the bits of n. The formulas used are

F[2k+1] = 4F[k]^2 - F[k-1]^2 + 2(-1)^k F[2k-1] = F[k]^2 + F[k-1]^2

F[2k] = F[2k+1] - F[2k-1]

  1. What they mean by " binary powering algorithm"

  2. Is this the same as the "Fast Doubling method"

  3. Do they change the number to binary? and why?

    
    #include <stdio.h>
    #include "gmp-impl.h"
    #include "longlong.h"
    
    
    /* change to "#define TRACE(x) x" to get some traces */
    #define TRACE(x)
    
    
    /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
       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.)
    
       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^GMP_NUMB_BITS 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.
    
       When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low
       limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used.  */
    
    void
    mpz_fib_ui (mpz_ptr fn, unsigned long n)
    {
      mp_ptr         fp, xp, yp;
      mp_size_t      size, xalloc;
      unsigned long  n2;
      mp_limb_t      c;
      TMP_DECL;
    
      if (n <= FIB_TABLE_LIMIT)
        {
          MPZ_NEWALLOC (fn, 1)[0] = FIB_TABLE (n);
          SIZ(fn) = (n != 0);      /* F[0]==0, others are !=0 */
          return;
        }
    
      n2 = n/2;
      xalloc = MPN_FIB2_SIZE (n2) + 1;
      fp = MPZ_NEWALLOC (fn, 2 * xalloc);
    
      TMP_MARK;
      TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
      size = mpn_fib2_ui (xp, yp, n2);
    
      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)
        {
          /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k  */
          mp_size_t  xsize, ysize;
    
    #if HAVE_NATIVE_mpn_add_n_sub_n
          xp[size] = mpn_lshift (xp, xp, size, 1);
          yp[size] = 0;
          ASSERT_NOCARRY (mpn_add_n_sub_n (xp, yp, xp, yp, size+1));
          xsize = size + (xp[size] != 0);
          ASSERT (yp[size] <= 1);
          ysize = size + yp[size];
    #else
          mp_limb_t  c2;
    
          c2 = mpn_lshift (fp, xp, size, 1);
          c = c2 + mpn_add_n (xp, fp, yp, size);
          xp[size] = c;
          xsize = size + (c != 0);
          c2 -= mpn_sub_n (yp, fp, yp, size);
          yp[size] = c2;
          ASSERT (c2 <= 1);
          ysize = size + c2;
    #endif
    
          size = xsize + ysize;
          c = mpn_mul (fp, xp, xsize, yp, ysize);
    
    #if GMP_NUMB_BITS >= BITS_PER_ULONG
          /* no overflow, see comments above */
          ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2);
          fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
    #else
          if (n & 2)
        {
          ASSERT (fp[0] >= 2);
          fp[0] -= 2;
        }
          else
        {
          ASSERT (c != GMP_NUMB_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
        {
          /* F[2k] = F[k]*(F[k]+2F[k-1]) */
    
          mp_size_t  xsize, ysize;
    #if HAVE_NATIVE_mpn_addlsh1_n
          c = mpn_addlsh1_n (yp, xp, yp, size);
    #else
          c = mpn_lshift (yp, yp, size, 1);
          c += mpn_add_n (yp, yp, xp, size);
    #endif
          yp[size] = c;
          xsize = size;
          ysize = size + (c != 0);
          size += ysize;
          c = mpn_mul (fp, yp, ysize, xp, xsize);
        }
    
      /* one or two high zeros */
    +
    −  size -= (c == 0);
      size -= (fp[size-1] == 0);
      SIZ(fn) = size;
    
      TRACE (printf ("done special, size=%ld\n", size);
         mpn_trace ("fp ", fp, size));
    
      TMP_FREE;
    }


Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source