Problem 921: Golden Recurrence

Problem 921 source code

C++

#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kMod = 398874989ULL;
constexpr u64 kGroup = kMod * kMod - 1ULL;

struct Elem {
    u64 a;
    u64 b;
};

u64 add_mod(u64 x, u64 y, u64 mod) {
    x += y;
    if (x >= mod) {
        x -= mod;
    }
    return x;
}

u64 mul_mod(u64 x, u64 y, u64 mod) {
    return static_cast<u64>((static_cast<u128>(x) * y) % mod);
}

Elem mul_elem(const Elem& x, const Elem& y) {
    const u64 aa = mul_mod(x.a, y.a, kMod);
    const u64 bb = mul_mod(x.b, y.b, kMod);
    const u64 real = add_mod(aa, (5ULL * bb) % kMod, kMod);
    const u64 imag = add_mod(mul_mod(x.a, y.b, kMod), mul_mod(x.b, y.a, kMod), kMod);
    return {real, imag};
}

std::array<Elem, 64> precompute_powers(const Elem& base) {
    std::array<Elem, 64> powers{};
    powers[0] = base;
    for (int i = 1; i < 64; ++i) {
        powers[i] = mul_elem(powers[i - 1], powers[i - 1]);
    }
    return powers;
}

Elem pow_from_powers(u64 exp, const std::array<Elem, 64>& powers) {
    Elem result{1, 0};
    while (exp != 0) {
        const int bit = __builtin_ctzll(exp);
        result = mul_elem(result, powers[bit]);
        exp &= (exp - 1);
    }
    return result;
}

Elem pow_elem(Elem base, u64 exp) {
    Elem result{1, 0};
    while (exp != 0) {
        if (exp & 1ULL) {
            result = mul_elem(result, base);
        }
        exp >>= 1U;
        if (exp != 0) {
            base = mul_elem(base, base);
        }
    }
    return result;
}

u64 pow_mod(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 value = base % mod;
    while (exp != 0) {
        if (exp & 1ULL) {
            result = mul_mod(result, value, mod);
        }
        exp >>= 1U;
        if (exp != 0) {
            value = mul_mod(value, value, mod);
        }
    }
    return result;
}

u64 pow5_mod(u64 x) {
    const u64 x2 = mul_mod(x, x, kMod);
    const u64 x4 = mul_mod(x2, x2, kMod);
    return mul_mod(x4, x, kMod);
}

u64 s_from_exp(u64 exp, const std::array<Elem, 64>& base_powers) {
    const Elem value = pow_from_powers(exp, base_powers);
    const u64 p = value.b;
    const u64 q = (value.a == 0 ? 0 : kMod - value.a);
    return add_mod(pow5_mod(p), pow5_mod(q), kMod);
}

u64 solve() {
    constexpr int kM = 1'618'034;
    const Elem base{(kMod + kMod - 2) % kMod, 1};  // -2 + sqrt(5)
    const auto base_powers = precompute_powers(base);

    u64 exp_prev2 = 5;  // 5^{F_1}
    u64 exp_prev1 = 5;  // 5^{F_2}

    u64 answer = s_from_exp(exp_prev1, base_powers);  // i=2
    for (int i = 3; i <= kM; ++i) {
        const u64 current = mul_mod(exp_prev1, exp_prev2, kGroup);
        answer = add_mod(answer, s_from_exp(current, base_powers), kMod);
        exp_prev2 = exp_prev1;
        exp_prev1 = current;
    }
    return answer;
}

void validate() {
    const Elem base{(kMod + kMod - 2) % kMod, 1};
    const auto base_powers = precompute_powers(base);

    assert(s_from_exp(1, base_powers) == 33);
    assert(s_from_exp(5, base_powers) == 257933744);
    assert(s_from_exp(25, base_powers) == 26500067);

    Elem stepped = base;
    u64 exp = 1;
    for (int n = 0; n <= 20; ++n) {
        const Elem by_exp = pow_from_powers(exp, base_powers);
        assert(by_exp.a == stepped.a && by_exp.b == stepped.b);
        stepped = pow_elem(stepped, 5);
        exp = mul_mod(exp, 5, kGroup);
    }

    u64 fib_prev2 = 1;
    u64 fib_prev1 = 1;
    u64 fexp_prev2 = 5;
    u64 fexp_prev1 = 5;
    for (int i = 3; i <= 25; ++i) {
        const u64 fib_cur = fib_prev1 + fib_prev2;
        const u64 fexp_cur = mul_mod(fexp_prev1, fexp_prev2, kGroup);
        assert(fexp_cur == pow_mod(5, fib_cur, kGroup));
        fib_prev2 = fib_prev1;
        fib_prev1 = fib_cur;
        fexp_prev2 = fexp_prev1;
        fexp_prev1 = fexp_cur;
    }
}

}  // namespace

int main() {
    validate();
    std::cout << solve() << '\n';
    return 0;
}

Python

def solve():
    MOD = 398874989
    GROUP = MOD * MOD - 1
    kM = 1618034

    def add_mod(x, y): return (x + y) % MOD
    def mul_mod(x, y): return x * y % MOD

    def mul_elem(x, y):
        aa = x[0]*y[0]%MOD; bb = x[1]*y[1]%MOD
        real = (aa + 5*bb) % MOD
        imag = (x[0]*y[1] + x[1]*y[0]) % MOD
        return (real, imag)

    def pow_elem(base, exp):
        r = (1, 0)
        while exp > 0:
            if exp & 1: r = mul_elem(r, base)
            exp >>= 1
            if exp > 0: base = mul_elem(base, base)
        return r

    def pow5(x):
        x2 = x*x%MOD; x4 = x2*x2%MOD; return x4*x%MOD

    base = ((MOD + MOD - 2) % MOD, 1)
    # Precompute powers of base for fast exponentiation
    powers = [base]
    for _ in range(63): powers.append(mul_elem(powers[-1], powers[-1]))

    def pow_base(exp):
        r = (1, 0)
        while exp:
            b = (exp & -exp).bit_length() - 1
            r = mul_elem(r, powers[b]); exp &= exp - 1
        return r

    def s_from_exp(exp):
        v = pow_base(exp)
        p = v[1]; q = (MOD - v[0]) % MOD
        return (pow5(p) + pow5(q)) % MOD

    ep2 = 5; ep1 = 5  # 5^{F_1}, 5^{F_2}
    answer = s_from_exp(ep1)
    for i in range(3, kM + 1):
        cur = ep1 * ep2 % GROUP
        answer = (answer + s_from_exp(cur)) % MOD
        ep2, ep1 = ep1, cur

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler921 {

    static final long MOD = 398874989L;
    static final long GROUP = MOD * MOD - 1L;

    static class Elem {
        long a, b;

        Elem(long a, long b) {
            this.a = a;
            this.b = b;
        }
    }

    static long addMod(long x, long y, long mod) {
        long res = x + y;
        if (res >= mod)
            res -= mod;
        return res;
    }

    static long mulMod(long x, long y, long mod) {
        long res = 0;
        long a = x;
        long b = y;
        while (b > 0) {
            if ((b & 1L) != 0) {
                res = addMod(res, a, mod);
            }
            a = addMod(a, a, mod);
            b >>= 1L;
        }
        return res;
    }

    static Elem mulElem(Elem x, Elem y) {
        long aa = mulMod(x.a, y.a, MOD);
        long bb = mulMod(x.b, y.b, MOD);
        long real = addMod(aa, mulMod(5L, bb, MOD), MOD);
        long imag = addMod(mulMod(x.a, y.b, MOD), mulMod(x.b, y.a, MOD), MOD);
        return new Elem(real, imag);
    }

    static Elem[] precomputePowers(Elem base) {
        Elem[] powers = new Elem[64];
        powers[0] = base;
        for (int i = 1; i < 64; ++i) {
            powers[i] = mulElem(powers[i - 1], powers[i - 1]);
        }
        return powers;
    }

    static Elem powFromPowers(long exp, Elem[] powers) {
        Elem result = new Elem(1, 0);
        int bit = 0;
        while (exp != 0) {
            if ((exp & 1L) != 0) {
                result = mulElem(result, powers[bit]);
            }
            exp >>= 1L;
            bit++;
        }
        return result;
    }

    static long pow5Mod(long x) {
        long x2 = mulMod(x, x, MOD);
        long x4 = mulMod(x2, x2, MOD);
        return mulMod(x4, x, MOD);
    }

    static long sFromExp(long exp, Elem[] powers) {
        Elem value = powFromPowers(exp, powers);
        long p = value.b;
        long q = value.a == 0 ? 0 : MOD - value.a;
        return addMod(pow5Mod(p), pow5Mod(q), MOD);
    }

    public static String solve() {
        int kM = 1618034;
        Elem base = new Elem((MOD + MOD - 2) % MOD, 1);
        Elem[] basePowers = precomputePowers(base);

        long expPrev2 = 5;
        long expPrev1 = 5;

        long answer = sFromExp(expPrev1, basePowers);
        for (int i = 3; i <= kM; ++i) {
            long current = mulMod(expPrev1, expPrev2, GROUP);
            answer = addMod(answer, sFromExp(current, basePowers), MOD);
            expPrev2 = expPrev1;
            expPrev1 = current;
        }

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}