Now that all is done—yes, with the exception of the grammar, I know, as a lot of other things, too—that we never did Fibonacci numbers. I’m not full sure about the pronunciation of his name (although two Italian men with different dialects give an example) but quite sure about the fast computation of the individual values.

The normal implementation is by following the definition:

static unsigned long fibonacci(unsigned long n){ unsigned long i = 1, j = 0, k, l; for (k = 1; k <= n; k++) { l = i + j; i = j; j = l; } return j; }

This is good for up to for 32-bit `unsigned long`

.

We want more, of course, at least.

One fast algorithm is by matrix exponentiation which is fast but has a lot of redundancies.

These can be stripped easily and leaves us with the following:

which assumes that and are known which is not the biggest problem, I think. I had to massage it a bit for easier implementation with libtommath but I hope it is still readable.

int mp_fibonacci(unsigned long n, mp_int * r) { unsigned long i = n - 1; mp_int a, b, c, d, t, t1, t2, t3; int e; if ((e = mp_init_multi(&a, &b, &c, &d, &t, &t1, &t2, &t3, NULL)) != MP_OKAY) { return e; } mp_set(&a, (mp_digit) 1); mp_set(&b, (mp_digit) 0); mp_set(&c, (mp_digit) 0); mp_set(&d, (mp_digit) 1); while (i > 0) { if (i & 0x1) { //t = d*(a + b) + c*b; if ((e = mp_mul(&c, &b, &t1)) != MP_OKAY) { return e; } if ((e = mp_add(&a, &b, &t2)) != MP_OKAY) { return e; } if ((e = mp_mul(&d, &t2, &t3)) != MP_OKAY) { return e; } if ((e = mp_add(&t3, &t1, &t)) != MP_OKAY) { return e; } //a = d*b + c*a; if ((e = mp_mul(&d, &b, &t1)) != MP_OKAY) { return e; } if ((e = mp_mul(&c, &a, &t2)) != MP_OKAY) { return e; } if ((e = mp_add(&t1, &t2, &a)) != MP_OKAY) { return e; } //b = t; if ((e = mp_copy(&t, &b)) != MP_OKAY) { return e; } } //t = d*(2*c + d); if ((e = mp_mul_2d(&c, 1, &t1)) != MP_OKAY) { return e; } if ((e = mp_add(&t1, &d, &t2)) != MP_OKAY) { return e; } if ((e = mp_mul(&d, &t2, &t)) != MP_OKAY) { return e; } //c = c*c + d*d; if ((e = mp_sqr(&c, &t1)) != MP_OKAY) { return e; } if ((e = mp_sqr(&d, &t2)) != MP_OKAY) { return e; } if ((e = mp_add(&t1, &t2, &c)) != MP_OKAY) { return e; } //d = t; if ((e = mp_copy(&t, &d)) != MP_OKAY) { return e; } i /= 2; } if ((e = mp_add(&a, &b, r)) != MP_OKAY) { return e; } mp_clear_multi(&a, &b, &c, &d, &t, &t1, &t2, &t3, NULL); return MP_OKAY; }

This function is implemented as `bn_mp_fibonacci.c`

in my fork of libtommath.

Once we got the Fibonacci numbers the Lucas numbers fall out quasi automatically.

We have two ways to solve the problem. One is the obvious and most probably also the faster one and is based on the fact that Lucas numbers lie between two consecutive Fibonacci numbers and can be calculated in this way:

Another way is by way of the closed form for the Fibonacci numbers, known as Binet’s formula

From this follows for Lucas numbers

Which makes it obvious now[1] that we can have

take the square root of both sides, keep the principal values, drop the rest and get

It is one evaluation of a Fibonacci number replaced with a square root and a squaring.

It’s slower but less than expected.

The implementation for the Lucas number evaluation with the very first algorithm is

int mp_lucas(unsigned long n, mp_int * c){ mp_int fnm1, fnp1; int e; if ((e = mp_init_multi(&fnm1, &fnp1, NULL)) != MP_OKAY) { return e; } if ((e = mp_fibonacci(n - 1, &fnm1)) != MP_OKAY) { return e; } if ((e = mp_fibonacci(n + 1, &fnp1)) != MP_OKAY) { return e; } if ((e = mp_add(&fnm1, &fnp1, c)) != MP_OKAY) { return e; } return MP_OKAY; }

The above code has been merged as `bn_mp_lucas.c`

into my fork of libtommath.

[1] took me a couple of sheets of paper and two pencils but once you got it it is obviously—obvious.