How can I calculate Fibonacci(A) mod B with A,B<=10^18 in C++ (g++)

Question

For example, A=10^17, B=10^17 (<64bits)
Typically, in the algorithm above, the calculs to compute F(2n) and those to compute F(2n+1) exceeds long long int types and we can't use modular computation in it.
The best algorithm to compute it I talk is fibonacci fast doubling:

F(0) = 0, F(1) = 1.
F(2n) = F(n)(2*F(n+1) – F(n)).
F(2n + 1) = F(n)² + F(n+1)².

Do you know some types in new C++14 (g++8.3.0 or llvm-clang C++) to use to avoid overflow.
I tried __float128 that is best than long double with no success. (see the g++ code above)
I have heard of the existence of __int128 and __int256 with no printf possibilities but I haven't try it.
Are they availailable in g++ 8.3.0 or are there other fast means to handle 128bits ints to do intermediate calculs you can think of? (time perfs are important)

#include <bits/stdc++.h>
using namespace std;

__float128 a,b,c,d;
long long mod;

void fast_fib(long long n,long long ans[]){
    if(n == 0){
        ans[0] = 0;
        ans[1] = 1;
        return;
    }
    fast_fib((n/2),ans);
    a = ans[0];             /* F(n) */
    b = ans[1];             /* F(n+1) */
    c = 2*b - a;
    if(c < 0) c += mod;
    c = (a * c);      /* F(2n) */
    while(c>=mod)c-=mod;
    d = (a*a + b*b);  /* F(2n + 1) */
    while(d>=mod)d-=mod;
    if(n%2 == 0){
        ans[0] = c;
        ans[1] = d;
    }
    else{
        ans[0] = d;
        ans[1] = c+d;
    }
}

int main(){
  int T=1000;
  long long n;
  while(T--){
    scanf("%lld %lld",&n,&mod);
    long long ans[2]={0};
    fast_fib(n,ans);
    printf("%lld\n", ans[0]);
  }
  return 0;
}

with __float128 I can't implement the modulo efficiently and a,b,c,d must store 128 bits data.

Answer 1

You don't need any floating point type for the calculations. You can use long long type only. First, you need a function that multiplicates two long long numbers (that are less than 10^18) modulo B. This can be done with the similar to exponentiation by squaring method:

long long multiply(long long a, long long b, long long M) {
    long long res = 0;
    long long d = a;

    while (b > 0) {
        if (b & 1) {
            res = (res + d) % M;
        }
        b /= 2;
        d = (d + d) % M;
    }
    return res;
}

Second, you need to add modulo operation to almost all of your arithmetic operations. And you definitely need to replace these loops while(c>=mod)c-=mod (they could be very slow) with the addition of % mod to the corresponding operations.

Your code with __float_128 replaced with long long and with proper modular arithmetic: https://ideone.com/t6R7Tf

---

Another thing you can do is to use (as was mentioned in the comments) Boost.Multiprecision or non-standard __int128 type (if supported) instead of long long type with complicated multiplication.

---

Also, you could use a slightly different (but using the same math actually) approach that seems more obvious to me - the Fibonacci numbers matrix formula

To calculate the Nth power of matrix you can use exponentiation by squaring doing all operations modulo B.