Skip to content

Instantly share code, notes, and snippets.

@yodalee
Created February 11, 2019 17:25
Show Gist options
  • Save yodalee/4e221b081be4b367e9c7ef328ada7db5 to your computer and use it in GitHub Desktop.
Save yodalee/4e221b081be4b367e9c7ef328ada7db5 to your computer and use it in GitHub Desktop.
Fast Fibonacci and Fomula Solution
// 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;
}
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment