Fast Fibonacci and Fomula Solution

fastfib.c

// file: fastfib.c
// compile with gcc -lgmp fastfib.c

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <ctype.h>

// return pair F(n), F(n+1) with mpz_t in argument
void recur_fib(uint64_t n, mpz_t fn, mpz_t fnp1) {
  if (n == 0) {
    mpz_set_ui(fn, 0);
    mpz_set_ui(fnp1, 1);
  } else {
    mpz_t fn_half1, fn_half2, temp_sq;

    mpz_inits(fn_half1, fn_half2, temp_sq, NULL);
    recur_fib(n >> 1, fn_half1, fn_half2);

    if (n % 2 == 0) { // n is even
      mpz_mul(temp_sq, fn_half1, fn_half1);

      // F(2n) = F(n) * (2*F(n+1) - F(n)).
      mpz_mul(fn, fn_half2, fn_half1);
      mpz_mul_ui(fn, fn, 2);
      mpz_sub(fn, fn, temp_sq);

      // F(2n+1) = F(n+1)^2 + F(n)^2.
      mpz_mul(fnp1, fn_half2, fn_half2);
      mpz_add(fnp1, fnp1, temp_sq);
    } else { // n is odd
      mpz_mul(temp_sq, fn_half2, fn_half2);

      // F(2n+1) = F(n+1)^2 + F(n)^2.
      mpz_mul(fn, fn_half1, fn_half1);
      mpz_add(fn, fn, temp_sq);

      // F(2n+2) = F(n+1)^2 + 2*F(n)*F(n+1).
      mpz_mul(fnp1, fn_half1, fn_half2);
      mpz_mul_ui(fnp1, fnp1, 2);
      mpz_add(fnp1, fnp1, temp_sq);
    }

    mpz_clears(fn_half1, fn_half2, temp_sq, NULL);
  }
}

// initialize F(n),F(n+1) and call recur_fib
void fast_fib(uint64_t n) {
  mpz_t fib_n, fib_np1;
  mpz_inits(fib_n, fib_np1, NULL);
  recur_fib(n, fib_n, fib_np1);
  gmp_printf ("fib(%d) = %Zd\n", n, fib_n);
  mpz_clears(fib_n, fib_np1, NULL);
}

int main(int argc, char *argv[])
{
  if (argc != 2) {
    fprintf(stderr, "Usage: <exe> <n>\n");
    return -1;
  }
  uint64_t n = strtoull(argv[1], NULL, 10);
  fast_fib(n);

  return 0;
}

formulafib.c

// file: formulafib.c
// compile with gcc -lmpfr -lm -lgmp formulafib.c

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <mpfr.h>
#include <gmp.h>
#include <ctype.h>

const double LOG2 = 0.3010299956639812L;

long double logfib(uint64_t n) {
  // estimated by F(n) < 2^n
  return LOG2 * n;
}

void calc_fib(uint64_t n) {
  mpfr_t sqrt5, phi1, phi2;
  mpz_t fib;
  uint64_t precision;
  uint64_t digits = ceill(logfib(n));

  precision = ceill(logfib(n)*logl(10)/logl(2));

  // init
  mpz_init(fib);
  mpfr_inits2(precision, sqrt5, phi1, phi2, NULL);

  // sqrt5
  mpfr_sqrt_ui(sqrt5, 5, MPFR_RNDN);

  // phi1
  mpfr_add_ui(phi1, sqrt5, 1, MPFR_RNDN);
  mpfr_div_ui(phi1, phi1,  2, MPFR_RNDN);
  mpfr_pow_ui(phi1, phi1,  n, MPFR_RNDN);

  // phi2
  mpfr_ui_sub(phi2, 1, sqrt5, MPFR_RNDN);
  mpfr_div_ui(phi2, phi2,  2, MPFR_RNDN);
  mpfr_pow_ui(phi2, phi2,  n, MPFR_RNDN);

  // combine
  mpfr_sub(phi1, phi1, phi2, MPFR_RNDN);
  mpfr_div(phi1, phi1, sqrt5, MPFR_RNDN);

  mpfr_get_z(fib, phi1, MPFR_RNDN);
  gmp_printf ("fib(%d) = %Zd\n", n, fib);

  // clear
  mpfr_clears(sqrt5, phi1, phi2, NULL);
  mpz_clear(fib);
}

int main(int argc, char *argv[])
{
  if (argc != 2) {
    fprintf(stderr, "Usage: <exe> <n>\n");
    return -1;
  }
  uint64_t n = strtoull(argv[1], NULL, 10);
  calc_fib(n);

  return 0;
}