Adventures of a Programmer: Parser Writing Peril XI

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 F_{47} for 32-bit unsigned long.
We want more, of course,F_{10^6} at least.
One fast algorithm is by matrix exponentiation which is fast but has a lot of redundancies.

{\displaystyle{ \left[\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}\right]^n = \left[\begin{matrix} F_{n+1} & F_{n} \\ F_{n} & F_{n-1} \end{matrix} \right] }}

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

{\displaystyle{ \begin{aligned} F_{2n} &= F_{n} \bigl( 2F_{n+1} - F_{n}) \bigr) \\ F_{2n+1} &= F_{n+1}^2 + F_{n}^2 \end{aligned} }}

which assumes that F_{n} and F_{n+1} 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:

{\displaystyle{ L_n = F_{n-1}+F_{n+1}=F_n+2F_{n-1} }}

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

{\displaystyle{  \frac{\varphi^n-(-\varphi)^{-n})}{\sqrt{5}} }}

From this follows for Lucas numbers

{\displaystyle{\varphi^n = {{L_n + F_n \sqrt{5}} \over 2} }}

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

{\displaystyle{L_n^2 = 5 F_n^2 + 4 (-1)^n }}

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

{\displaystyle{L_n = \left|\left(5 F_n^2 + 4 (-1)^n \right)^\frac{1}{2}\right|}}

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.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s