#include <stdio.h>
    #include <time.h>
    
    int fib_mod_9837(int n) {
        int a = 0;
        int b = 1;
        int k = 1;
        while (k < n) {
            k <<= 1;
        }
        while (k) {
            int a2 = a*a;
            int b2 = b*b;
            int ab2 = a*b << 1;
            a = a2 + ab2;
            b = a2 + b2;
            if (n & k) {
                int t = a;
                a += b;
                b = t;
            }
            a %= 9837;
            b %= 9837;
            k >>= 1;
        }
        return a;
    }
    
    int main() {
        struct timespec begin, end;
        clock_gettime(CLOCK_MONOTONIC_RAW, &begin);
        int result = fib_mod_9837(99999999);
        clock_gettime(CLOCK_MONOTONIC_RAW, &end);
        printf("%lld\n", result);
        printf("elapsed nanos: %lld\n", (end.tv_nsec - begin.tv_nsec) + (end.tv_sec  - begin.tv_sec)*1000000000);
    }

Unintentionally, the article is a good cautionary tale of how algorithmic knowledge can beat raw hardware power, if you know what you're doing.

It wouldn’t surprise me at all if doing those 5 matrix multiplications on a CPU were faster than doing 31 on the GPU.

Also, FTA: “To calculate extremly large fibonacci numbers we need to avoid integer overflow. A simple way to avoid that is to just take everything modulo a given number in the matrix multiplication“

⇒ This also is NOT calculating the Fibonacci numbers, copping out when things get difficult on a GPU.

If you want to compute Fn mod some constant for 1 ≤ n ≤ some large constant fast on a GPU, I think I would compute quite a few of the intermediate matrices using the fast way, then spin up some ‘threads’ to compute the intermediate versions by copying those with the ‘standard’ matrix.

In my experience where I didn't modulo the numbers, the compute time was dominated by the last couple of multiplications with Megabyte-sized integers.

Or recursive hashing?