Issue finding last digits of the nth term in the Fibonacci sequence

Question

I wrote some code to determine the nth Fibonacci number using the nice blog post given in the accepted answer to this question: Finding out nth fibonacci number for very large 'n'. I am doing this as a way of practising a more difficult recursion problem given on projecteuler but that is not really relevant. The method relies on changing the problem to a small linear algebra problem of the form

Fn = T^n F1

where F1 = (1 1)^t and Fn contains the nth and (n-1)th Fibonacci number. The term T^n can then be determined in O(log n) time. I implemented this succesfully and it seems to work fine. When I perform the matrix exponentiation I use %10000 so I get the last 4 digits only, which seems to work (I checked against some large Fibonacci numbers). However, I wanted to try to get more last digits by increasing the number 10000. This doesn't seem to work however. I no longer get the correct answer. Here is my code

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

const unsigned long M = 10000;

unsigned long int * matProd(unsigned long int * A, unsigned long int * B){
  unsigned long int * C;
  C = malloc(4*sizeof(unsigned long int));
  C[0] = ((A[0]*B[0]%M) + (A[1]*B[2]%M)) % M;
  C[1] = ((A[0]*B[1]%M) + (A[1]*B[3]%M)) % M;
  C[2] = ((A[2]*B[0]%M) + (A[3]*B[2]%M)) % M;
  C[3] = ((A[2]*B[1]%M) + (A[3]*B[3]%M)) % M;
  return C;
}

unsigned long int * matExp(unsigned long int *A, unsigned long int n){
  if (n==1){
    return A;
  }
  if (n%2==0){
    return matExp(matProd(A,A),n/2);
  }
  return matProd(A,matExp(A,n-1));
}

unsigned long int findFib(unsigned long int n){
  unsigned long int A[4] = {0, 1, 1, 1};
  unsigned long int * C;
  C = malloc(4*sizeof(unsigned long int));
  C = matExp(A,n-2);
  return (C[2]+C[3]);
}

main(){
  unsigned long int n = 300;
  printf("%ld\n",findFib(n));
}

There are probably several issues there with regards to proper coding conventions and things that can be improved. I thought changing to long int might solve the problem but this does not seem to do the trick. So basically the problem is that increasing M to for instance 1000000 does not give me more digits but instead gives me nonsense. What mistake am I making?

P.S. sorry for the poor math formatting, I am used to math.stackexchange.

score 1 · Accepted Answer · answered Jul 22 '15 at 13:35

1

The issue is probably that you are running on a system where long is 32-bits in size, as I believe is the case for Windows. You can check this by compiling and running printf("%d\n", sizeof(long)) which should output 4.

Since with M=1000000=10^6, the product of two numbers smaller than M can go up to 10^12, you get overflow issues when you are computing your matrix entries since unsigned long can hold up to at most 2^32-1 or roughly 4 * 10^9.

To fix this simply using unsigned long long as your type instead of unsigned long. Or better yet, uint64_t, which is guaranteed to be 64-bits in all platforms (and which will require #include <stdint.h>). This should make your code work for M up to sqrt(2^64)~10^9. If you need bigger than that you'll need to use a big integer library.

answered Jul 22 '15 at 13:35

toth

2,519
1
15
23

Thanks for the answer, this indeed works. I was considering perhaps defining my own recursive function for multiplying two numbers which finds a*b by adding b to itself a times, each time performing the modulo operation. This will be horridly slow probably but should alleviate the problems with overflow. Anyway, thanks for the helpful answer! – Slugger Jul 22 '15 at 13:55
Unclear as why `uint64_t` is _better_ than `unsigned long long` here. `unsigned long long` is at least 64-bit and `uint64_t` need not exist (albeit on rare machines) to be compliant. If anything, code could use `uint_least64_t`. – chux - Reinstate Monica Jul 22 '15 at 15:51

score 1 · Answer 2 · answered Jul 22 '15 at 13:36

If the program works for M == 10000 but fails for M == 1000000 (or even for M == 100000) then that probably means that your C implementation's unsigned long int type is 32 bits wide.

If your matrix elements are drawn exclusively from Z¹⁰⁰⁰⁰, then they require at most 14 significant binary digits. The products you compute in your matrix multiplication function, before reducing modulo M, may therefore require up to 28 binary digits. If you increase M even to 100000, however, then the matrix elements require up to 17 binary digits, and the intermediate products require up to 34. The reduction modulo M is too late to prevent that overflowing a 32-bit integer and therefore giving you garbage results.

You could consider declaring the element type as uint64_t instead. If it's an overflow problem then that should give you enough extra digits to handle M == 1000000.

Issue finding last digits of the nth term in the Fibonacci sequence

2 Answers2