I want to determine the nth Fibonacci term in the sequence using large integer values

Question

The code below is able to determine the correct sequence up to a point namely 70 using the data type unsigned long long. I know the sequence can become large thus I mod 10,000 the results. I want to determine the nth term for 10,000 using the best data type or improve the algo to calculate the nth term.

#define MOD %10000

unsigned long long calc(long nth) {
    return (pow( 1 + sqrt(5), nth ) - pow( 1 - sqrt(5), nth )) / (pow(2.0, nth)*(sqrt(5)));
}

int main() {
    long t, nth;
    for (std::cin>>t;  t-- && std::cin>>nth; ) {
        std::cout<<calc(nth-2)MOD<<" "<<calc(nth-1)MOD<<" "<<calc(nth)MOD<<std::endl;
    }   
    return 0;
}

Answer 1

Your algorithm will not compute the correct result for large N's, due to the floating point errors of sqrn(5).

In order to speed up your algorithm you can use fast doubling Fibonacci:

   F(2k) = F(k)[2F(k+1) - F(k)]
   F(2k+1) = F(k+1)^2 + F(k)^2

Applying modulo arithmetics, your final fastest algorithm would be:

   F(2k) = F(k)[2F(k+1) - F(k)] % 10000
   F(2k+1) = (F(k+1)^2 + F(k)^2) % 10000

Using this approach, your function never exceeds 10000, thus an int type suffices.

EDIT: Okay I had some free time on a Friday night (not a good thing I guess) and implemented the algorithm. I implemented two versions, first one with O(1) memory and O(lg n) time complexity and second one using a cache, with memory and worst-case runtime of O(lg n), but with a best case runtime of O(1).

#include <iostream>
#include <unordered_map>

using namespace std;

const int P = 10000;

/* Fast Fibonacci with O(1) memory and O(lg n) time complexity. No cache. */

int fib_uncached (int n)
{
    /* find MSB position */
    int msb_position = 31;
    while (!((1 << (msb_position-1) & n)) && msb_position >= 0)
        msb_position--;

    int a=0, b=1; 

    for (int i=msb_position; i>=0;i--)
    {       
        int d = (a%P) * ((b%P)*2 - (a%P) + P),
            e = (a%P) * (a%P) + (b%P)*(b%P);
        a=d%P;
        b=e%P;

        if (((n >> i) & 1) != 0)
        {
            int c = (a + b) % P;
            a = b;
            b = c;
        }
    }
    return a;
}  

/* Fast Fibonacci using cache */
int fib (int n)
{
    static std::unordered_map<int,int> cache;

    if (cache.find(n) == cache.end()) 
    {
        int f;
        if (n==0)
            f = 0;
        else if (n < 3)
            f = 1;
        else if (n % 2 == 0)
        {
            int k = n/2;
            f = (fib(k) * (2*fib(k+1) - fib(k))) % P;
        } 
        else
        {
            int k = (n-1)/2;
            f = (fib(k+1)*fib(k+1)+ fib(k) * fib(k)) % P;
        }
        if (f<0)
            f += P;

        cache[n] = f;
    }
    return cache.at(n);
}

int main ()
{
    int i ;
    cin >> i;
    cout << i << " : " << fib(i) << endl;
return 0;
}

Reference for cache-less implementations: https://www.nayuki.io/page/fast-fibonacci-algorithms

Answer 2

Calculate the terms successively, taking the mod at each step. Since each term only depends on the previous two, your computational space is just a 3-element array.

#include <iostream>
using namespace std;

typedef unsigned long long numtype;

const numtype MOD = 10000;

// Assume n is 1-based

numtype fib(int n)
{
  numtype seq[3] = {1,1,2};
  if( --n < 3 ) return seq[n]; // make n 0-based
  for( int i=3 ; i<=n ; ++i )
  {
    seq[i%3] = (seq[(i-1)%3] + seq[(i-2)%3]) % MOD;
    cout << seq[i%3] << ' '; // comment out for large n
  }
  return seq[n%3];
}

int main()
{
//numtype answer = fib(10000000); // 6875
  numtype answer = fib(70); // 9135
  cout << endl;
  cout << "answer = " << answer << endl;
}

Output:

3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 946 7711 8657 6368 5025 1393 6418 7811 4229 2040 6269 8309 4578 2887 7465 352 7817 8169 5986 4155 141 4296 4437 8733 3170 1903 5073 6976 2049 9025 1074 99 1173 1272 2445 3717 6162 9879 6041 5920 1961 7881 9842 7723 7565 5288 2853 8141 994 9135

(The first 3 terms 1,1,2 are intentionally missing.) The 70th term is 9135.

Time is O(n) and memory is O(1).