Created
February 11, 2019 17:25
-
-
Save yodalee/4e221b081be4b367e9c7ef328ada7db5 to your computer and use it in GitHub Desktop.
Fast Fibonacci and Fomula Solution
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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; | |
} | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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