Problem 356: Largest Roots of Cubic Polynomials

View on Project Euler

Project Euler Problem 356 Solution

EulerSolve provides an optimized solution for Project Euler Problem 356, Largest Roots of Cubic Polynomials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each \(n\in\{1,\dots,30\}\), let \(a_n\) be the largest real root of $$x^3-2^n x^2+n=0.$$ The task is to compute $$\sum_{n=1}^{30}\left\lfloor a_n^{987654321}\right\rfloor \pmod{10^8}.$$ A direct floating-point evaluation of such enormous powers would be fragile and unnecessary. The local C++, Python, and Java solutions instead turn the problem into an exact linear recurrence modulo \(10^8\). Mathematical Approach Step 1: Locate the Three Real Roots Write \(c=2^n\) and \(f(x)=x^3-cx^2+n\). The implementations rely on the fact that the cubic has three real roots, with one large positive root and two small roots near the origin. Indeed, $$f(-1)=n-1-c<0,\qquad f(0)=n>0,$$ so there is a root \(\gamma\in(-1,0)\). Also, $$f(1)=n+1-c\le 0,\qquad f(0)=n>0,$$ so there is a root \(\beta\in(0,1]\), with \(\beta=1\) only when \(n=1\). Finally, $$f(c)=n>0,\qquad f(c-1)=n-(c-1)^2\le 0,$$ so the largest root \(\alpha=a_n\) lies in \((c-1,c)\)....

Detailed mathematical approach

Problem Summary

For each \(n\in\{1,\dots,30\}\), let \(a_n\) be the largest real root of

$$x^3-2^n x^2+n=0.$$

The task is to compute

$$\sum_{n=1}^{30}\left\lfloor a_n^{987654321}\right\rfloor \pmod{10^8}.$$

A direct floating-point evaluation of such enormous powers would be fragile and unnecessary. The local C++, Python, and Java solutions instead turn the problem into an exact linear recurrence modulo \(10^8\).

Mathematical Approach

Step 1: Locate the Three Real Roots

Write \(c=2^n\) and \(f(x)=x^3-cx^2+n\). The implementations rely on the fact that the cubic has three real roots, with one large positive root and two small roots near the origin.

Indeed,

$$f(-1)=n-1-c<0,\qquad f(0)=n>0,$$

so there is a root \(\gamma\in(-1,0)\). Also,

$$f(1)=n+1-c\le 0,\qquad f(0)=n>0,$$

so there is a root \(\beta\in(0,1]\), with \(\beta=1\) only when \(n=1\). Finally,

$$f(c)=n>0,\qquad f(c-1)=n-(c-1)^2\le 0,$$

so the largest root \(\alpha=a_n\) lies in \((c-1,c)\). Thus the three roots satisfy

$$\gamma\in(-1,0),\qquad \beta\in(0,1],\qquad \alpha\in(2^n-1,2^n).$$

By Vieta's formulas,

$$\alpha+\beta+\gamma=c,\qquad \alpha\beta+\alpha\gamma+\beta\gamma=0,\qquad \alpha\beta\gamma=-n.$$

Step 2: Define the Power Sums

Let

$$S_k=\alpha^k+\beta^k+\gamma^k.$$

Because every root \(r\in\{\alpha,\beta,\gamma\}\) satisfies

$$r^3=cr^2-n,$$

multiplying by \(r^{k-3}\) gives

$$r^k=cr^{k-1}-nr^{k-3}\qquad (k\ge 3).$$

Summing over the three roots yields the recurrence used in every implementation:

$$\boxed{S_k=cS_{k-1}-nS_{k-3}\qquad (k\ge 3).}$$

The initial values are immediate:

$$S_0=3,$$

$$S_1=\alpha+\beta+\gamma=c,$$

$$S_2=(\alpha+\beta+\gamma)^2-2(\alpha\beta+\alpha\gamma+\beta\gamma)=c^2.$$

So the entire problem is reduced to a third-order linear recurrence.

Step 3: Why \(\left\lfloor a_n^m\right\rfloor=S_m-1\) for Odd \(m\)

Let \(m\) be odd and define the small-root contribution

$$\rho_m=\beta^m+\gamma^m.$$

Since \(\gamma<0\) and \(m\) is odd, we have \(\gamma^m<0\). We also know

$$\beta+\gamma=c-\alpha,$$

and because \(\alpha\in(c-1,c)\), this implies

$$0<\beta+\gamma<1.$$

In particular, \(\beta>|\gamma|\). For odd \(m\), monotonicity of \(x\mapsto x^m\) on positive numbers gives

$$\beta^m>|\gamma|^m,$$

hence

$$0<\rho_m=\beta^m-|\gamma|^m.$$

At the same time, \(\beta\le 1\) and \(|\gamma|<1\), so

$$\rho_m<1.$$

Therefore

$$0<\rho_m<1,$$

and since \(\alpha^m=S_m-\rho_m\), we obtain the crucial identity

$$\boxed{\left\lfloor \alpha^m\right\rfloor=\left\lfloor S_m-\rho_m\right\rfloor=S_m-1.}$$

This is exactly why the production code never needs to approximate \(a_n^{987654321}\) directly.

Step 4: Convert the Recurrence to Matrix Exponentiation

The recurrence has fixed order 3, so the state vector

$$\begin{bmatrix}S_k\\S_{k-1}\\S_{k-2}\end{bmatrix}$$

evolves by

$$\begin{bmatrix}S_k\\S_{k-1}\\S_{k-2}\end{bmatrix} = \begin{bmatrix} c & 0 & -n\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix}S_{k-1}\\S_{k-2}\\S_{k-3}\end{bmatrix}.$$

Hence, for \(m\ge 2\),

$$\begin{bmatrix}S_m\\S_{m-1}\\S_{m-2}\end{bmatrix} = M^{m-2} \begin{bmatrix}S_2\\S_1\\S_0\end{bmatrix}.$$

The solutions compute this power by binary exponentiation, reducing every multiplication modulo \(10^8\). In modular arithmetic the entry \(-n\) is stored as

$$10^8-(n\bmod 10^8),$$

which is why the source code writes the transition matrix with a non-negative last coefficient.

Step 5: A Small Checkpoint

Take \(n=2\), so \(c=4\). Then

$$S_0=3,\qquad S_1=4,\qquad S_2=16.$$

For the odd exponent \(m=3\), the recurrence gives

$$S_3=4\cdot 16-2\cdot 3=58.$$

Since \(0<\beta^3+\gamma^3<1\), we conclude

$$\left\lfloor a_2^3\right\rfloor=S_3-1=57.$$

The C++ implementation performs this kind of spot check numerically for small odd exponents before relying on the modular fast path for the full input.

How the Code Works

The three language versions all implement the same pipeline. First, powMod computes \(c=2^n\bmod 10^8\). Then multiplyMatrix and powerMatrix raise the \(3\times 3\) transition matrix to the power \(m-2\). The function computeNewtonSumMod returns \(S_m\bmod 10^8\), handling the base cases \(m=0,1,2\) directly and using matrix exponentiation otherwise.

After that, computeTermMod applies the identity \(\lfloor a_n^m\rfloor\equiv S_m-1\pmod{10^8}\), and solve sums the terms for \(n=1\) through \(30\). The longer C++ file also contains validation helpers: it checks the statement sample \(\lfloor a_2^2\rfloor=14\), verifies the matrix recurrence against a direct iterative recurrence, compares the floor identity against floating-point root calculations on small cases, and tests thread consistency. Those checks justify the method, but the final answer path itself is entirely integer-based.

Complexity Analysis

For one fixed \(n\), the dominant cost is exponentiating a \(3\times 3\) matrix to power \(m-2\), which takes \(O(\log m)\) matrix multiplications. The problem has only 30 values of \(n\), so the total running time is

$$O(30\log m)=O(\log m),$$

with \(O(1)\) auxiliary memory. The C++ version can parallelize the 30 independent \(n\)-values, but that changes only wall-clock time, not asymptotic complexity.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=356
  2. Vieta's formulas: Wikipedia - Vieta's formulas
  3. Newton sums / power sums: Wikipedia - Newton's identities
  4. Matrix exponentiation for linear recurrences: cp-algorithms - binary exponentiation of recurrence transitions

Problem 356 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr u64 kMod = 100'000'000ULL;
constexpr u32 kDefaultNMin = 1U;
constexpr u32 kDefaultNMax = 30U;
constexpr u64 kDefaultExponent = 987'654'321ULL;
constexpr std::size_t kAutoThreadThreshold = 96ULL;

struct Options {
    u32 n_min = kDefaultNMin;
    u32 n_max = kDefaultNMax;
    u64 exponent = kDefaultExponent;
    bool run_checks = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
        return false;
    }
    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg, const char* prefix, unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checks") {
            options.run_checks = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--n-min=", parsed_u32)) {
            options.n_min = parsed_u32;
            continue;
        }
        if (parse_u32_after_prefix(arg, "--n-max=", parsed_u32)) {
            options.n_max = parsed_u32;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--exp=", parsed_u64)) {
            options.exponent = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.n_min == 0U || options.n_max == 0U) {
        std::cerr << "n bounds must be >= 1.\n";
        return false;
    }
    if (options.n_min > options.n_max) {
        std::cerr << "Invalid range: --n-min must be <= --n-max.\n";
        return false;
    }
    if (options.exponent == 0ULL || (options.exponent & 1ULL) == 0ULL) {
        std::cerr << "Exponent must be a positive odd integer.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    const unsigned max_threads = static_cast<unsigned>(
        std::min<std::size_t>(workload_units, static_cast<std::size_t>(std::numeric_limits<unsigned>::max())));

    if (requested_threads > 0U) {
        return std::max(1U, std::min(requested_threads, max_threads));
    }

    if (workload_units < kAutoThreadThreshold) {
        return 1U;
    }

    unsigned hardware_threads = std::thread::hardware_concurrency();
    if (hardware_threads == 0U) {
        hardware_threads = 1U;
    }
    return std::max(1U, std::min(hardware_threads, max_threads));
}

u64 pow_mod_u64(u64 base, u64 exponent, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }

    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    u64 exp = exponent;

    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % mod);
        exp >>= 1ULL;
    }

    return result;
}

struct Matrix3 {
    u64 v[3][3] = {{0ULL, 0ULL, 0ULL}, {0ULL, 0ULL, 0ULL}, {0ULL, 0ULL, 0ULL}};
};

Matrix3 make_identity_matrix3() {
    Matrix3 out;
    for (int i = 0; i < 3; ++i) {
        out.v[i][i] = 1ULL;
    }
    return out;
}

Matrix3 multiply_matrix3(const Matrix3& a, const Matrix3& b, const u64 mod) {
    Matrix3 out;
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            u128 sum = 0;
            for (int k = 0; k < 3; ++k) {
                sum += static_cast<u128>(a.v[i][k]) * static_cast<u128>(b.v[k][j]);
            }
            out.v[i][j] = static_cast<u64>(sum % mod);
        }
    }
    return out;
}

Matrix3 power_matrix3(Matrix3 base, u64 exponent, const u64 mod) {
    Matrix3 result = make_identity_matrix3();
    u64 exp = exponent;

    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = multiply_matrix3(result, base, mod);
        }
        base = multiply_matrix3(base, base, mod);
        exp >>= 1ULL;
    }

    return result;
}

u64 compute_newton_sum_mod(const u32 n, const u64 exponent, const u64 mod) {
    const u64 c = pow_mod_u64(2ULL, static_cast<u64>(n), mod);
    const u64 s0 = 3ULL % mod;
    const u64 s1 = c;
    const u64 s2 = static_cast<u64>((static_cast<u128>(c) * c) % mod);

    if (exponent == 0ULL) {
        return s0;
    }
    if (exponent == 1ULL) {
        return s1;
    }
    if (exponent == 2ULL) {
        return s2;
    }

    Matrix3 transition;
    transition.v[0][0] = c;
    transition.v[0][1] = 0ULL;
    transition.v[0][2] = (mod - (static_cast<u64>(n) % mod)) % mod;
    transition.v[1][0] = 1ULL % mod;
    transition.v[1][1] = 0ULL;
    transition.v[1][2] = 0ULL;
    transition.v[2][0] = 0ULL;
    transition.v[2][1] = 1ULL % mod;
    transition.v[2][2] = 0ULL;

    const Matrix3 p = power_matrix3(transition, exponent - 2ULL, mod);
    const u128 top = static_cast<u128>(p.v[0][0]) * s2 +
                     static_cast<u128>(p.v[0][1]) * s1 +
                     static_cast<u128>(p.v[0][2]) * s0;
    return static_cast<u64>(top % mod);
}

u64 compute_newton_sum_mod_iterative(const u32 n, const u64 exponent, const u64 mod) {
    const u64 c = pow_mod_u64(2ULL, static_cast<u64>(n), mod);
    const u64 n_mod = static_cast<u64>(n) % mod;
    const u64 s0 = 3ULL % mod;
    const u64 s1 = c;
    const u64 s2 = static_cast<u64>((static_cast<u128>(c) * c) % mod);

    if (exponent == 0ULL) {
        return s0;
    }
    if (exponent == 1ULL) {
        return s1;
    }
    if (exponent == 2ULL) {
        return s2;
    }

    u64 a = s0;
    u64 b = s1;
    u64 cur = s2;
    for (u64 k = 3ULL; k <= exponent; ++k) {
        const u64 lhs = static_cast<u64>((static_cast<u128>(c) * cur) % mod);
        const u64 rhs = static_cast<u64>((static_cast<u128>(n_mod) * a) % mod);
        u64 next = lhs;
        if (next >= rhs) {
            next -= rhs;
        } else {
            next += mod - rhs;
        }
        a = b;
        b = cur;
        cur = next;
    }
    return cur;
}

u128 compute_newton_sum_exact_small(const u32 n, const u32 exponent) {
    u128 c = 1;
    for (u32 i = 0; i < n; ++i) {
        c *= 2U;
    }

    const u128 s0 = 3U;
    const u128 s1 = c;
    const u128 s2 = c * c;

    if (exponent == 0U) {
        return s0;
    }
    if (exponent == 1U) {
        return s1;
    }
    if (exponent == 2U) {
        return s2;
    }

    u128 a = s0;
    u128 b = s1;
    u128 cur = s2;
    for (u32 k = 3U; k <= exponent; ++k) {
        const u128 next = c * cur - static_cast<u128>(n) * a;
        a = b;
        b = cur;
        cur = next;
    }
    return cur;
}

long double evaluate_polynomial(const long double x, const long double c, const u32 n) {
    return x * x * x - c * x * x + static_cast<long double>(n);
}

long double largest_real_root(const u32 n) {
    const long double c = std::ldexp(1.0L, static_cast<int>(n));
    long double low = std::max(0.0L, c - 1.0L);
    long double high = c;

    while (evaluate_polynomial(low, c, n) > 0.0L && low > 0.0L) {
        low = std::max(0.0L, low - 1.0L);
    }

    for (int iter = 0; iter < 220; ++iter) {
        const long double mid = (low + high) * 0.5L;
        if (evaluate_polynomial(mid, c, n) <= 0.0L) {
            low = mid;
        } else {
            high = mid;
        }
    }
    return (low + high) * 0.5L;
}

u64 floor_power(const long double base, const u32 exponent) {
    long double value = 1.0L;
    for (u32 i = 0; i < exponent; ++i) {
        value *= base;
    }
    const long double floored = std::floor(value + 1e-18L);
    return (floored <= 0.0L) ? 0ULL : static_cast<u64>(floored);
}

u64 compute_term_mod(const u32 n, const u64 exponent) {
    const u64 sum_mod = compute_newton_sum_mod(n, exponent, kMod);
    return (sum_mod + kMod - 1ULL) % kMod;
}

u64 solve_range(const u32 n_min,
                const u32 n_max,
                const u64 exponent,
                const bool allow_multithreading,
                const unsigned requested_threads) {
    const std::size_t workload = static_cast<std::size_t>(n_max) - static_cast<std::size_t>(n_min) + 1ULL;
    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload);

    if (threads <= 1U) {
        u64 total = 0ULL;
        for (u32 n = n_min; n <= n_max; ++n) {
            total += compute_term_mod(n, exponent);
            if (total >= kMod) {
                total -= kMod;
            }
        }
        return total;
    }

    std::atomic<u32> next_n(n_min);
    std::vector<u64> partial(static_cast<std::size_t>(threads), 0ULL);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    auto worker = [&](const unsigned tid) {
        u64 local = 0ULL;
        while (true) {
            const u32 n = next_n.fetch_add(1U, std::memory_order_relaxed);
            if (n > n_max) {
                break;
            }
            local += compute_term_mod(n, exponent);
            if (local >= kMod) {
                local -= kMod;
            }
        }
        partial[static_cast<std::size_t>(tid)] = local;
    };

    for (unsigned t = 0U; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (std::thread& th : pool) {
        th.join();
    }

    u64 total = 0ULL;
    for (const u64 part : partial) {
        total += part;
        if (total >= kMod) {
            total -= kMod;
        }
    }
    return total;
}

u64 solve(const Options& options) {
    return solve_range(options.n_min,
                       options.n_max,
                       options.exponent,
                       options.allow_multithreading,
                       options.requested_threads);
}

bool run_validations(const Options& options) {
    // Checkpoint 1: statement example floor(a_2^2) = 14.
    {
        const long double a2 = largest_real_root(2U);
        const u64 floor_a2_sq = floor_power(a2, 2U);
        if (floor_a2_sq != 14ULL) {
            std::cerr << "Validation failed (1/4): floor(a_2^2) expected 14, got "
                      << floor_a2_sq << '\n';
            return false;
        }
        std::cout << "Validation 1/4 [PASS]\n";
    }

    // Checkpoint 2: matrix recurrence matches direct iterative recurrence.
    {
        for (u32 n = 1U; n <= 12U; ++n) {
            for (u64 k = 0ULL; k <= 220ULL; ++k) {
                const u64 by_matrix = compute_newton_sum_mod(n, k, kMod);
                const u64 by_iter = compute_newton_sum_mod_iterative(n, k, kMod);
                if (by_matrix != by_iter) {
                    std::cerr << "Validation failed (2/4): recurrence mismatch at n=" << n
                              << ", k=" << k << " (matrix=" << by_matrix
                              << ", iterative=" << by_iter << ")\n";
                    return false;
                }
            }
        }
        std::cout << "Validation 2/4 [PASS]\n";
    }

    // Checkpoint 3: floor(a_n^m) == S_m(n)-1 for small odd-exponent spot checks.
    {
        const u32 odd_exponents[] = {1U, 3U, 5U, 7U, 9U, 11U};
        for (u32 n = 1U; n <= 3U; ++n) {
            const long double alpha = largest_real_root(n);
            for (const u32 m : odd_exponents) {
                const u64 floor_numeric = floor_power(alpha, m);
                const u128 exact_sum = compute_newton_sum_exact_small(n, m);
                const u64 floor_via_identity = static_cast<u64>(exact_sum - 1U);
                if (floor_numeric != floor_via_identity) {
                    std::cerr << "Validation failed (3/4): floor identity mismatch at n=" << n
                              << ", m=" << m << " (numeric=" << floor_numeric
                              << ", identity=" << floor_via_identity << ")\n";
                    return false;
                }
            }
        }
        std::cout << "Validation 3/4 [PASS]\n";
    }

    // Checkpoint 4: thread consistency on a larger reduced workload.
    if (options.allow_multithreading) {
        unsigned hardware_threads = std::thread::hardware_concurrency();
        if (hardware_threads == 0U) {
            hardware_threads = 1U;
        }

        if (hardware_threads > 1U) {
            const u32 check_n_min = 1U;
            const u32 check_n_max = 400U;
            const u64 single = solve_range(check_n_min, check_n_max, options.exponent, false, 1U);
            const unsigned compare_threads =
                (options.requested_threads > 1U) ? options.requested_threads : 2U;
            const u64 multi =
                solve_range(check_n_min, check_n_max, options.exponent, true, compare_threads);
            if (single != multi) {
                std::cerr << "Validation failed (4/4): single-thread vs multi-thread mismatch"
                          << " (single=" << single << ", multi=" << multi << ")\n";
                return false;
            }
            std::cout << "Validation 4/4 [PASS]\n";
        } else {
            std::cout << "Validation 4/4 [SKIP: single-core environment]\n";
        }
    } else {
        std::cout << "Validation 4/4 [SKIP: --single-thread enabled]\n";
    }

    return true;
}

void print_usage() {
    std::cerr << "Usage: ./Euler356 [--skip-checks] [--single-thread] [--threads=<unsigned>] "
                 "[--n-min=<u32>] [--n-max=<u32>] [--exp=<odd-u64>]\n";
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        print_usage();
        return 1;
    }

    if (options.run_checks && !run_validations(options)) {
        return 1;
    }

    const u64 answer = solve(options);
    std::cout << std::setfill('0') << std::setw(8) << answer << '\n';
    return 0;
}

Python

def pow_mod(base, exp, mod):
    if mod == 1: return 0
    res = 1
    cur = base % mod
    while exp > 0:
        if exp & 1:
            res = (res * cur) % mod
        cur = (cur * cur) % mod
        exp >>= 1
    return res

def multiply_matrix(a, b, mod):
    out = [[0]*3 for _ in range(3)]
    for i in range(3):
        for j in range(3):
            s = 0
            for k in range(3):
                s += a[i][k] * b[k][j]
            out[i][j] = s % mod
    return out

def power_matrix(base, exp, mod):
    res = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    while exp > 0:
        if exp & 1:
            res = multiply_matrix(res, base, mod)
        base = multiply_matrix(base, base, mod)
        exp >>= 1
    return res

def compute_newton_sum_mod(n, exp, mod):
    c = pow_mod(2, n, mod)
    s0 = 3 % mod
    s1 = c
    s2 = (c * c) % mod
    
    if exp == 0: return s0
    if exp == 1: return s1
    if exp == 2: return s2
    
    transition = [
        [c, 0, (mod - (n % mod)) % mod],
        [1, 0, 0],
        [0, 1, 0]
    ]
    p = power_matrix(transition, exp - 2, mod)
    top = p[0][0] * s2 + p[0][1] * s1 + p[0][2] * s0
    return top % mod

def compute_term_mod(n, exp, mod):
    sum_mod = compute_newton_sum_mod(n, exp, mod)
    return (sum_mod + mod - 1) % mod

def solve():
    mod = 100000000
    total = 0
    for n in range(1, 31):
        total = (total + compute_term_mod(n, 987654321, mod)) % mod
    return str(total)

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

Java

public class Euler356 {
    static long powMod(long base, long exp, long mod) {
        if (mod == 1)
            return 0;
        long res = 1;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * cur) % mod;
            cur = (cur * cur) % mod;
            exp >>= 1;
        }
        return res;
    }

    static long[][] multiplyMatrix(long[][] a, long[][] b, long mod) {
        long[][] out = new long[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                long s = 0;
                for (int k = 0; k < 3; k++) {
                    s += a[i][k] * b[k][j];
                }
                out[i][j] = s % mod;
            }
        }
        return out;
    }

    static long[][] powerMatrix(long[][] base, long exp, long mod) {
        long[][] res = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = multiplyMatrix(res, base, mod);
            base = multiplyMatrix(base, base, mod);
            exp >>= 1;
        }
        return res;
    }

    static long computeNewtonSumMod(int n, long exp, long mod) {
        long c = powMod(2, n, mod);
        long s0 = 3 % mod;
        long s1 = c;
        long s2 = (c * c) % mod;

        if (exp == 0)
            return s0;
        if (exp == 1)
            return s1;
        if (exp == 2)
            return s2;

        long[][] transition = {
                { c, 0, (mod - (n % mod)) % mod },
                { 1, 0, 0 },
                { 0, 1, 0 }
        };
        long[][] p = powerMatrix(transition, exp - 2, mod);
        long top = p[0][0] * s2 + p[0][1] * s1 + p[0][2] * s0;
        return top % mod;
    }

    public static String solve() {
        long mod = 100000000;
        long total = 0;
        for (int n = 1; n <= 30; n++) {
            long term = (computeNewtonSumMod(n, 987654321, mod) + mod - 1) % mod;
            total = (total + term) % mod;
        }
        return String.format("%08d", total);
    }

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