Created
February 14, 2017 12:40
-
-
Save Randl/fbe86e02fd096dc48c6af99976f2c7e8 to your computer and use it in GitHub Desktop.
Calculate nth Fibonacci number
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
#include <math.h> | |
#include <stdint.h> | |
#include <stdio.h> | |
#include <gmp.h> | |
/* | |
* from Fibo(n) = phi^n / sqrt(5) approximate number of bits in nth fibo | |
* is log_2(phi) * n - log_2(5)/2 For n >> 1, it's a bit less than 0.695n | |
*/ | |
#define FIBO_BITCOUNT(x) ((mp_bitcnt_t)ceil(0.695 * (x))) | |
// for string length same calculation applies with base 10 | |
#define FIBO_STRING_LENGHT(x) ((mp_bitcnt_t)ceil(0.21 * (x))) | |
void fibo(mpz_t *a, mp_bitcnt_t n) { | |
if (n < 3) { | |
mpz_init2(*a, 2); | |
mpz_set_ui(*a, 1); | |
return; | |
} | |
const mp_bitcnt_t num_of_bits = FIBO_BITCOUNT(n); | |
mpz_t b, tmp1, tmp2; | |
mpz_init2(*a, num_of_bits); | |
mpz_init2(b, num_of_bits); | |
mpz_init2(tmp1, num_of_bits); | |
mpz_init2(tmp2, num_of_bits); | |
mpz_set_ui(*a, 1); // F(1) | |
mpz_set_ui(b, 1); // F(2) | |
printf("Num of bits: %d\n", num_of_bits); | |
// going from MSB to LSB, i is a number of current bit, counting from 1 | |
for (int_fast16_t i = sizeof(n) * 8 - 1 - __builtin_clzll(n); i > 0; --i) { | |
// A is F(k) | |
mpz_mul_ui(tmp1, b, 2); // 2F(k+1) | |
mpz_sub(tmp1, tmp1, *a); // 2F(k+1) - F(k) | |
mpz_mul(tmp1, tmp1, *a); // F(2k) = F(k) * (2F(k+1) - F(k)) | |
mpz_pow_ui(tmp2, *a, 2); // F(k)^2 | |
mpz_pow_ui(b, b, 2); // F(k+1)^2 | |
mpz_add(b, tmp2, b); // F(2k+1) = F(k+1)^2+F(k)^2 | |
mpz_swap(tmp1, *a); | |
if ((n >> (i - 1)) % 2) { // If current bit is 1, advance one more | |
mpz_add(tmp1, *a, b); // tmp1 = F(2k+2) | |
mpz_swap(*a, b); // a = F(2k+1) | |
mpz_swap(tmp1, b); // b = F(2k+2) | |
} | |
} | |
mpz_clears(b, tmp1, tmp2, NULL); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment