-
Notifications
You must be signed in to change notification settings - Fork 3
/
fibonacci.c
63 lines (60 loc) · 1.64 KB
/
fibonacci.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <gmp.h>
#include <stdio.h>
void fibonacci(mpz_t res, long n) {
/*
On intialize le résultat à M⁰, i.e
|0 1|⁰
|1 1|
(i.e., l'identité)
*/
mpz_t a00, a01, a10, a11;
mpz_init_set_ui(a00, 1); mpz_init_set_ui(a01, 0);
mpz_init_set_ui(a10, 0); mpz_init_set_ui(a11, 1);
// des variables auxiliaires
mpz_t tmp1, tmp2;
mpz_init(tmp1); mpz_init(tmp2);
// On calcule M^n par un square and multiply de gauche à droite
long i;
for (i = 1 ; i < n ; i <<= 1);
for ( ; i >= 1 ; i >>= 1) {
/*
Square (générique) :
|a₀₀ a₀₁|² = |a₀₀²+a₀₁a₁₀ a₀₁(a₀₀+a₁₁)|
|a₁₀ a₁₁| |a₁₀(a₀₀+a₁₁) a₀₁a₁₀+a₁₁²|
*/
mpz_add(tmp1, a00, a11);
mpz_mul(tmp2, a01, a10);
mpz_mul(a01, a01, tmp1);
mpz_mul(a10, a10, tmp1);
mpz_mul(tmp1, a00, a00);
mpz_add(a00, tmp1, tmp2);
mpz_mul(tmp1, a11, a11);
mpz_add(a11, tmp1, tmp2);
if (n & i) {
/*
Multiply (spécifique à M) :
|a₀₀ a₀₁| |0 1| = |a₀₁ a₀₀+a₀₁|
|a₁₀ a₁₁| |1 1| |a₁₁ a₁₀+a₁₁|
*/
mpz_add(a01, a00, a01);
mpz_add(a11, a10, a11);
mpz_sub(a00, a01, a00);
mpz_sub(a10, a11, a10);
}
}
// f(n) apparaît sur la diagonale de M
mpz_set(res, a01);
}
void main(int argc, char** argv) {
printf("GNU GMP %d.%d.%d, compiled with %s %s\n", __GNU_MP_VERSION,
__GNU_MP_VERSION_MINOR, __GNU_MP_VERSION_PATCHLEVEL,
__GMP_CC, __GMP_CFLAGS);
if (argc > 1) {
long n = atol(argv[1]);
mpz_t f;
mpz_init(f);
fibonacci(f, n);
gmp_printf("%Zd\n", f);
mpz_clear(f);
}
}