Problem 833: Square Triangle Products

Project Euler Problem 833 Solution

EulerSolve provides an optimized solution for Project Euler Problem 833, Square Triangle Products, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T_n=\frac{n(n+1)}{2}$ be the $n$-th triangular number. The problem asks for the sum $S(N)$ of all positive integers $c\le N$ for which there exist two distinct triangular numbers $T_a$ and $T_b$ such that $$T_aT_b=c^2.$$ For the actual input $N=10^{35}$, direct enumeration of triangular pairs is impossible. The successful approach groups triangular numbers by squarefree part, turns each group into a Pell-type family, and then sums the corresponding square roots modulo $136101521$. Mathematical Approach The key is to understand when two triangular numbers have a square product. That happens exactly when they belong to the same squarefree class. Step 1: Reduce the Condition to a Pell-Type Equation If $T_aT_b$ is a square, then $T_a$ and $T_b$ have the same squarefree part. So we can write $$T_a=d\,r^2,\qquad T_b=d\,s^2,$$ where $d$ is squarefree. Then automatically $$c=d\,rs.$$ Thus, for each fixed squarefree $d$, we need all triangular numbers of the form $$T_n=d\,y^2.$$ Using $T_n=\frac{n(n+1)}{2}$ and $x=2n+1$, this becomes $$x^2-8dy^2=1.$$ So every squarefree class of triangular numbers is controlled by a Pell equation. Step 2: Describe One Family with Chebyshev Polynomials Suppose a squarefree class has a fundamental representative $T_k=d$....

Detailed mathematical approach

Problem Summary

Let $T_n=\frac{n(n+1)}{2}$ be the $n$-th triangular number. The problem asks for the sum $S(N)$ of all positive integers $c\le N$ for which there exist two distinct triangular numbers $T_a$ and $T_b$ such that

$$T_aT_b=c^2.$$

For the actual input $N=10^{35}$, direct enumeration of triangular pairs is impossible. The successful approach groups triangular numbers by squarefree part, turns each group into a Pell-type family, and then sums the corresponding square roots modulo $136101521$.

Mathematical Approach

The key is to understand when two triangular numbers have a square product. That happens exactly when they belong to the same squarefree class.

Step 1: Reduce the Condition to a Pell-Type Equation

If $T_aT_b$ is a square, then $T_a$ and $T_b$ have the same squarefree part. So we can write

$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$

where $d$ is squarefree. Then automatically

$$c=d\,rs.$$

Thus, for each fixed squarefree $d$, we need all triangular numbers of the form

$$T_n=d\,y^2.$$

Using $T_n=\frac{n(n+1)}{2}$ and $x=2n+1$, this becomes

$$x^2-8dy^2=1.$$

So every squarefree class of triangular numbers is controlled by a Pell equation.

Step 2: Describe One Family with Chebyshev Polynomials

Suppose a squarefree class has a fundamental representative $T_k=d$. Then

$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$

The Pell family generated by this base produces all triangular numbers in the same squarefree class. The implementation expresses the square multipliers with Chebyshev polynomials of the second kind:

$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$

For this fixed base $k$, the family of triangular numbers is

$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$

These are exactly the triangular numbers whose squarefree part is $T_k$. The first term is $A_0(k)=T_k$ because $U_0=1$.

Step 3: Turn a Family into Square Triangle Products

Take two distinct members of the same family, say $A_i(k)$ and $A_j(k)$ with $i<j$. Their product is

$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$

Therefore each admissible contribution to $S(N)$ has the form

$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$

So the entire problem becomes: enumerate all family pairs $(i,j)$, sum $C_{i,j}(k)$ over all fundamental bases $k$, and keep only the values not exceeding $N$.

Step 4: Bound the Family Pairs and Sum over $k$

For fixed $i$ and $j$, the function $C_{i,j}(k)$ is increasing in $k$, because both $T_k$ and $U_m(2k+1)$ increase for $k\ge 1$. The smallest base is $k=1$, where $2k+1=3$, so a necessary condition for the pair $(i,j)$ to contribute at all is

$$U_i(3)\,U_j(3)\le N.$$

This is why the implementations first build the short sequence

$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$

and keep only the pair indices whose minimum possible contribution is not already too large.

Once $(i,j)$ is fixed, there is a largest admissible base $\kappa_{i,j}$ satisfying

$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$

Because $U_m(2k+1)$ is a polynomial in $k$ of degree $m$, the whole expression $C_{i,j}(k)$ is a polynomial in $k$. Hence

$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$

can be reduced to a linear combination of power sums $\sum_{k=1}^{\kappa_{i,j}} k^p$.

Step 5: Remove Non-Fundamental Bases

Not every base index $k$ should start a new squarefree class. Some triangular numbers already occur inside an older family:

$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$

for a smaller base $k_1$. Such $k$ are non-fundamental. If we summed over all $k\ge 1$ blindly, those classes would be counted more than once.

The duplicate base indices are generated by the same Pell mechanism. Starting from $x_1=2k_1+1$, define

$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$

Every $k_r$ produced in this way belongs to the squarefree class already represented by $k_1$, so its contribution must be subtracted from the all-$k$ polynomial sum.

Step 6: Worked Example

Take the smallest fundamental base $k=1$. Then

$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$

The corresponding triangular family begins with

$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$

Now pair these family members:

$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$

So $6$, $35$, and $210$ are all square triangle products.

The same family also explains why some later bases are non-fundamental. The Pell recurrence from $x_1=3$ gives

$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$

hence

$$k=0,1,8,49,\dots$$

and indeed

$$T_8=36,\qquad T_{49}=1225,$$

which were already present in the family of $k=1$. Those later bases must therefore be excluded as new starting points.

Step 7: Final Summation Formula

If $\mathcal{N}_{i,j}$ denotes the non-fundamental bases not exceeding the relevant bound for the pair $(i,j)$, then the target sum is

$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$

This is exactly the structure implemented by the programs.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they enumerate the admissible family pairs $(i,j)$ by growing the sequence $U_m(3)$ until its products are too large to matter. For each surviving pair they find the largest admissible base by monotone binary search on $C_{i,j}(k)\le N$.

Next they build the Chebyshev factors symbolically as polynomials in $k$. Multiplying those polynomials by the triangular factor $T_k=\frac{k(k+1)}{2}$ gives a polynomial representation of each contribution $C_{i,j}(k)$. The full range sum for that pair is then evaluated through power sums, and those power sums are obtained modulo $136101521$ by Lagrange-style interpolation of Faulhaber polynomials.

Finally the implementations generate every non-fundamental base only once, evaluate the corresponding pair contributions at those indices, and subtract them. The C++ implementation parallelizes this subtraction pass, while the Python and Java implementations perform the same arithmetic sequentially.

Complexity Analysis

Let $F$ be the number of family levels with $U_F(3)\le N$. Because $U_m(3)$ grows exponentially, $F=O(\log N)$, so the number of family pairs is $O(F^2)$. Let $B$ be the largest admissible base among all pairs, let $D$ be the maximum polynomial degree, and let $L$ be the number of non-fundamental bases up to $B$.

Enumerating pair indices costs $O(F^2)$. Finding all pair bounds costs $O(F^2\log B)$ evaluations, each using a short Chebyshev recurrence. Polynomial preprocessing and power-sum evaluation are low-degree operations with total cost polynomial in $F$, and the correction phase costs $O(F^2L)$ in the straightforward form used here. Memory is dominated by the stored pair metadata, polynomial coefficients, and the correction list, so it is $O(F^2D+L)$.

The important practical fact is that the family depth grows only logarithmically, which keeps the entire computation manageable even at $N=10^{35}$.

Footnotes and References

Project Euler Problem 833: https://projecteuler.net/problem=833
Triangular numbers: Wikipedia — Triangular number
Pell's equation: Wikipedia — Pell's equation
Chebyshev polynomials: Wikipedia — Chebyshev polynomials
Faulhaber's formula: Wikipedia — Faulhaber's formula
Lagrange interpolation: Wikipedia — Lagrange polynomial

Problem 833 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
#include <functional>

using namespace std;

namespace {

using int64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;

constexpr int64 kMod = 136101521;

struct PairInfo {
    int i_idx;
    int j_idx;
    u64 k_max;
    vector<int64> poly; // coefficients mod kMod
    int64 sum_all = 0;
    int64 sum_nf = 0;
};

u128 pow10_u128(int exp) {
    u128 res = 1;
    for (int i = 0; i < exp; ++i) res *= 10;
    return res;
}

bool mul_leq(u128 a, u128 b, u128 limit, u128 &prod) {
    if (a == 0 || b == 0) {
        prod = 0;
        return true;
    }
    if (a > limit / b) return false;
    prod = a * b;
    return true;
}

// Evaluate c = T_k * U_{i_idx}(x) * U_{j_idx}(x) with x=2k+1.
// Returns true if c <= limit.
bool eval_pair_leq(u64 k, int i_idx, int j_idx, u128 limit) {
    u128 kk = static_cast<u128>(k);
    u128 x = 2 * kk + 1;
    int max_idx = max(i_idx, j_idx);
    vector<u128> U(max_idx + 1, 0);
    U[0] = 1;
    if (max_idx >= 1) U[1] = 2 * x;

    u128 two_x = 2 * x;
    for (int m = 1; m < max_idx; ++m) {
        if (U[m] > (limit + U[m - 1]) / two_x) {
            U[m + 1] = limit + 1;
        } else {
            u128 term = two_x * U[m];
            u128 next = term - U[m - 1];
            U[m + 1] = (next > limit) ? (limit + 1) : next;
        }
    }

    u128 t = kk * (kk + 1) / 2;
    u128 tmp = 0;
    if (!mul_leq(t, U[i_idx], limit, tmp)) return false;
    if (!mul_leq(tmp, U[j_idx], limit, tmp)) return false;
    return true;
}

u128 eval_pair_value(u64 k, int i_idx, int j_idx) {
    u128 kk = static_cast<u128>(k);
    u128 x = 2 * kk + 1;
    int max_idx = max(i_idx, j_idx);
    vector<u128> U(max_idx + 1, 0);
    U[0] = 1;
    if (max_idx >= 1) U[1] = 2 * x;
    for (int m = 1; m < max_idx; ++m) {
        U[m + 1] = 2 * x * U[m] - U[m - 1];
    }
    u128 t = kk * (kk + 1) / 2;
    return t * U[i_idx] * U[j_idx];
}

u64 max_k_for_pair(u128 n, int i_idx, int j_idx) {
    u64 lo = 0;
    u64 hi = 1;
    while (eval_pair_leq(hi, i_idx, j_idx, n)) {
        if (hi > (numeric_limits<u64>::max() / 2)) break;
        hi *= 2;
    }
    while (lo + 1 < hi) {
        u64 mid = lo + (hi - lo) / 2;
        if (eval_pair_leq(mid, i_idx, j_idx, n)) {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    return lo;
}

vector<u64> generate_non_fundamentals(u64 k_max) {
    vector<u64> nf;
    if (k_max < 4) return nf;
    u64 k1_max = static_cast<u64>(std::sqrt(static_cast<long double>(k_max) / 4.0L)) + 2;
    for (u64 k1 = 1; k1 <= k1_max; ++k1) {
        u128 x1 = 2 * static_cast<u128>(k1) + 1;
        u128 x_prev = 1;
        u128 x_cur = x1;
        while (true) {
            u128 x_next = 2 * x1 * x_cur - x_prev;
            u64 k_next = static_cast<u64>((x_next - 1) / 2);
            if (k_next > k_max) break;
            nf.push_back(k_next);
            x_prev = x_cur;
            x_cur = x_next;
        }
    }
    sort(nf.begin(), nf.end());
    nf.erase(unique(nf.begin(), nf.end()), nf.end());
    return nf;
}

int64 mod_pow(int64 base, int64 exp) {
    int64 res = 1 % kMod;
    base %= kMod;
    while (exp > 0) {
        if (exp & 1) res = (static_cast<i128>(res) * base) % kMod;
        base = (static_cast<i128>(base) * base) % kMod;
        exp >>= 1;
    }
    return res;
}

int64 mod_inv(int64 a) {
    return mod_pow(a, kMod - 2);
}

using Poly = vector<int64>;

Poly poly_add(const Poly &a, const Poly &b) {
    size_t n = max(a.size(), b.size());
    Poly c(n, 0);
    for (size_t i = 0; i < n; ++i) {
        int64 va = (i < a.size()) ? a[i] : 0;
        int64 vb = (i < b.size()) ? b[i] : 0;
        int64 v = va + vb;
        if (v >= kMod) v -= kMod;
        c[i] = v;
    }
    return c;
}

Poly poly_sub(const Poly &a, const Poly &b) {
    size_t n = max(a.size(), b.size());
    Poly c(n, 0);
    for (size_t i = 0; i < n; ++i) {
        int64 va = (i < a.size()) ? a[i] : 0;
        int64 vb = (i < b.size()) ? b[i] : 0;
        int64 v = va - vb;
        if (v < 0) v += kMod;
        c[i] = v;
    }
    return c;
}

Poly poly_mul(const Poly &a, const Poly &b) {
    Poly c(a.size() + b.size() - 1, 0);
    for (size_t i = 0; i < a.size(); ++i) {
        for (size_t j = 0; j < b.size(); ++j) {
            c[i + j] = (c[i + j] + static_cast<i128>(a[i]) * b[j]) % kMod;
        }
    }
    return c;
}

Poly poly_scale(const Poly &a, int64 s) {
    Poly c(a.size(), 0);
    for (size_t i = 0; i < a.size(); ++i) {
        c[i] = static_cast<i128>(a[i]) * s % kMod;
    }
    return c;
}

vector<Poly> build_U_polys(int max_idx) {
    Poly t = {1, 2}; // 2k+1
    vector<Poly> U(max_idx + 1);
    U[0] = {1};
    if (max_idx >= 1) {
        U[1] = poly_scale(t, 2);
    }
    for (int m = 1; m < max_idx; ++m) {
        Poly temp = poly_mul(t, U[m]);
        temp = poly_scale(temp, 2);
        U[m + 1] = poly_sub(temp, U[m - 1]);
    }
    return U;
}

struct SumPowers {
    vector<vector<int64>> y;
    vector<int64> fact;
    vector<int64> invfact;
};

SumPowers precompute_sum_pows(int max_p) {
    SumPowers sp;
    sp.y.resize(max_p + 1);
    for (int p = 0; p <= max_p; ++p) {
        int m = p + 1;
        sp.y[p].assign(m + 1, 0);
        int64 s = 0;
        for (int i = 1; i <= m; ++i) {
            s += mod_pow(i, p);
            s %= kMod;
            sp.y[p][i] = s;
        }
    }
    sp.fact.assign(max_p + 2, 1);
    sp.invfact.assign(max_p + 2, 1);
    for (int i = 1; i < static_cast<int>(sp.fact.size()); ++i) {
        sp.fact[i] = static_cast<i128>(sp.fact[i - 1]) * i % kMod;
    }
    sp.invfact.back() = mod_inv(sp.fact.back());
    for (int i = static_cast<int>(sp.fact.size()) - 1; i >= 1; --i) {
        sp.invfact[i - 1] = static_cast<i128>(sp.invfact[i]) * i % kMod;
    }
    return sp;
}

int64 sum_pows(int p, u64 n, const SumPowers &sp) {
    if (n <= static_cast<u64>(p + 1)) {
        return sp.y[p][static_cast<size_t>(n)];
    }
    int m = p + 1;
    vector<int64> pre(m + 2, 1), suf(m + 2, 1);
    int64 nmod = static_cast<int64>(n % kMod);
    for (int i = 0; i <= m; ++i) {
        int64 term = nmod - i;
        if (term < 0) term += kMod;
        pre[i + 1] = static_cast<i128>(pre[i]) * term % kMod;
    }
    for (int i = m; i >= 0; --i) {
        int64 term = nmod - i;
        if (term < 0) term += kMod;
        suf[i] = static_cast<i128>(suf[i + 1]) * term % kMod;
    }
    int64 res = 0;
    for (int i = 0; i <= m; ++i) {
        int64 num = static_cast<i128>(pre[i]) * suf[i + 1] % kMod;
        int64 denom_inv = static_cast<i128>(sp.invfact[i]) * sp.invfact[m - i] % kMod;
        int64 term = static_cast<i128>(sp.y[p][i]) * num % kMod;
        term = static_cast<i128>(term) * denom_inv % kMod;
        if ((m - i) & 1) {
            if (term != 0) term = kMod - term;
        }
        res += term;
        if (res >= kMod) res -= kMod;
    }
    return res;
}

int64 sum_poly(const Poly &coeffs, u64 n, const SumPowers &sp) {
    int64 res = 0;
    for (size_t p = 0; p < coeffs.size(); ++p) {
        if (coeffs[p] == 0) continue;
        int64 spow = sum_pows(static_cast<int>(p), n, sp);
        res = (res + static_cast<i128>(coeffs[p]) * spow) % kMod;
    }
    return res;
}

vector<PairInfo> build_pairs(u128 n) {
    vector<u128> U3;
    U3.push_back(1);
    U3.push_back(6);
    while (true) {
        u128 next = 2 * static_cast<u128>(3) * U3[U3.size() - 1] - U3[U3.size() - 2];
        U3.push_back(next);
        if (next > n) break;
    }
    int max_idx = static_cast<int>(U3.size()) - 1;

    vector<PairInfo> pairs;
    for (int i = 0; i <= max_idx; ++i) {
        for (int j = i + 1; j <= max_idx; ++j) {
            u128 prod = 0;
            if (!mul_leq(U3[i], U3[j], n, prod)) continue;
            PairInfo info;
            info.i_idx = i;
            info.j_idx = j;
            pairs.push_back(info);
        }
    }
    return pairs;
}

u64 compute_exact_small(u64 n) {
    u128 n128 = static_cast<u128>(n);
    vector<PairInfo> pairs = build_pairs(n128);
    u64 k_max = 0;
    for (auto &p : pairs) {
        p.k_max = max_k_for_pair(n128, p.i_idx, p.j_idx);
        k_max = max(k_max, p.k_max);
    }
    vector<u64> nf = generate_non_fundamentals(k_max);
    vector<char> is_nf(k_max + 1, 0);
    for (u64 v : nf) if (v <= k_max) is_nf[v] = 1;

    u128 total = 0;
    for (const auto &p : pairs) {
        for (u64 k = 1; k <= p.k_max; ++k) {
            if (is_nf[k]) continue;
            u128 c = eval_pair_value(k, p.i_idx, p.j_idx);
            if (c > n128) break;
            total += c;
        }
    }
    return static_cast<u64>(total);
}

int64 compute_mod(u128 n) {
    vector<PairInfo> pairs = build_pairs(n);
    int max_idx = 0;
    u64 k_max = 0;
    for (auto &p : pairs) {
        p.k_max = max_k_for_pair(n, p.i_idx, p.j_idx);
        k_max = max(k_max, p.k_max);
        max_idx = max(max_idx, max(p.i_idx, p.j_idx));
    }

    vector<u64> nf = generate_non_fundamentals(k_max);

    const int64 inv2 = mod_inv(2);
    vector<Poly> U_polys = build_U_polys(max_idx);
    Poly t_poly = {0, inv2, inv2};
    int deg_max = 0;
    for (auto &p : pairs) {
        Poly r = poly_mul(U_polys[p.i_idx], U_polys[p.j_idx]);
        p.poly = poly_mul(r, t_poly);
        deg_max = max(deg_max, static_cast<int>(p.poly.size()) - 1);
    }

    SumPowers sp = precompute_sum_pows(deg_max);
    for (auto &p : pairs) {
        p.sum_all = sum_poly(p.poly, p.k_max, sp);
    }

    // Parallel non-fundamental contributions.
    const unsigned int hw = thread::hardware_concurrency();
    const unsigned int threads = (hw == 0) ? 4u : min(hw, 8u);
    vector<vector<int64>> local_sums(threads, vector<int64>(pairs.size(), 0));
    vector<thread> workers;

    auto worker = [&](unsigned int tid, size_t start, size_t end) {
        vector<int64> U(max_idx + 1, 0);
        for (size_t idx = start; idx < end; ++idx) {
            u64 k = nf[idx];
            u64 k_mod = k % kMod;
            int64 x_mod = static_cast<int64>((2 * (k_mod % kMod) + 1) % kMod);
            U[0] = 1;
            if (max_idx >= 1) U[1] = (2LL * x_mod) % kMod;
            for (int m = 1; m < max_idx; ++m) {
                int64 next = (2LL * x_mod % kMod * U[m] - U[m - 1]) % kMod;
                if (next < 0) next += kMod;
                U[m + 1] = next;
            }
            int64 t = static_cast<int64>(k_mod) * static_cast<int64>((k_mod + 1) % kMod) % kMod;
            t = static_cast<int64>(static_cast<i128>(t) * inv2 % kMod);

            auto &sum_nf = local_sums[tid];
            for (size_t p = 0; p < pairs.size(); ++p) {
                if (k > pairs[p].k_max) continue;
                int64 val = static_cast<int64>(static_cast<i128>(t) * U[pairs[p].i_idx] % kMod);
                val = static_cast<int64>(static_cast<i128>(val) * U[pairs[p].j_idx] % kMod);
                sum_nf[p] += val;
                if (sum_nf[p] >= kMod) sum_nf[p] -= kMod;
            }
        }
    };

    size_t chunk = (nf.size() + threads - 1) / threads;
    for (unsigned int t = 0; t < threads; ++t) {
        size_t start = t * chunk;
        size_t end = min(nf.size(), start + chunk);
        workers.emplace_back(worker, t, start, end);
    }
    for (auto &th : workers) th.join();

    for (size_t p = 0; p < pairs.size(); ++p) {
        int64 total = 0;
        for (unsigned int t = 0; t < threads; ++t) {
            total += local_sums[t][p];
            if (total >= kMod) total -= kMod;
        }
        pairs[p].sum_nf = total;
    }

    int64 ans = 0;
    for (const auto &p : pairs) {
        int64 term = p.sum_all - p.sum_nf;
        if (term < 0) term += kMod;
        ans += term;
        if (ans >= kMod) ans -= kMod;
    }
    return ans;
}

} // namespace

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    // Validation checks.
    if (compute_exact_small(100000) != 1479802ULL) {
        cerr << "Validation failed: S(1e5) mismatch\n";
        return 1;
    }
    if (compute_exact_small(1000000000ULL) != 241614948794ULL) {
        cerr << "Validation failed: S(1e9) mismatch\n";
        return 1;
    }

    u128 n = pow10_u128(35);
    int64 ans = compute_mod(n);
    cout << ans << "\n";
    return 0;
}

Python

import math

kMod = 136101521

def max_k_for_pair(n_val, i_idx, j_idx):
    def eval_pair_leq(k):
        x = 2 * k + 1
        max_idx = max(i_idx, j_idx)
        U = [0] * (max_idx + 1)
        U[0] = 1
        if max_idx >= 1: U[1] = 2 * x
        
        two_x = 2 * x
        for m in range(1, max_idx):
            if U[m] > (n_val + U[m - 1]) // two_x:
                U[m + 1] = n_val + 1
            else:
                nex = two_x * U[m] - U[m - 1]
                U[m + 1] = nex if nex <= n_val else n_val + 1
                
        t = k * (k + 1) // 2
        tmp = t * U[i_idx]
        if tmp > n_val: return False
        tmp *= U[j_idx]
        if tmp > n_val: return False
        return True

    lo = 0
    hi = 1
    while eval_pair_leq(hi):
        hi *= 2
    while lo + 1 < hi:
        mid = (lo + hi) // 2
        if eval_pair_leq(mid):
            lo = mid
        else:
            hi = mid
    return lo

def generate_non_fundamentals(k_max):
    nf = []
    if k_max < 4: return nf
    k1_max = int(math.sqrt(k_max / 4.0)) + 2
    for k1 in range(1, k1_max + 1):
        x1 = 2 * k1 + 1
        x_prev = 1
        x_cur = x1
        while True:
            x_next = 2 * x1 * x_cur - x_prev
            k_next = (x_next - 1) // 2
            if k_next > k_max: break
            nf.append(k_next)
            x_prev = x_cur
            x_cur = x_next
    return sorted(list(set(nf)))

def mod_pow(base, exp):
    res = 1 % kMod
    base %= kMod
    while exp > 0:
        if exp & 1: res = (res * base) % kMod
        base = (base * base) % kMod
        exp >>= 1
    return res

def mod_inv(a):
    return mod_pow(a, kMod - 2)

def poly_add(a, b):
    n = max(len(a), len(b))
    c = [0] * n
    for i in range(n):
        va = a[i] if i < len(a) else 0
        vb = b[i] if i < len(b) else 0
        c[i] = (va + vb) % kMod
    return c

def poly_sub(a, b):
    n = max(len(a), len(b))
    c = [0] * n
    for i in range(n):
        va = a[i] if i < len(a) else 0
        vb = b[i] if i < len(b) else 0
        c[i] = (va - vb) % kMod
    return c

def poly_mul(a, b):
    c = [0] * (len(a) + len(b) - 1)
    for i in range(len(a)):
        for j in range(len(b)):
            c[i + j] = (c[i + j] + a[i] * b[j]) % kMod
    return c

def poly_scale(a, s):
    return [(x * s) % kMod for x in a]

def build_U_polys(max_idx):
    t = [1, 2] # 2k+1
    U = [[0]] * (max_idx + 1)
    U[0] = [1]
    if max_idx >= 1:
        U[1] = poly_scale(t, 2)
    for m in range(1, max_idx):
        temp = poly_mul(t, U[m])
        temp = poly_scale(temp, 2)
        U[m + 1] = poly_sub(temp, U[m - 1])
    return U

class SumPowers:
    def __init__(self, y, fact, invfact):
        self.y = y
        self.fact = fact
        self.invfact = invfact

def precompute_sum_pows(max_p):
    y = [[0] * (p + 2) for p in range(max_p + 1)]
    for p in range(max_p + 1):
        m = p + 1
        s = 0
        for i in range(1, m + 1):
            s = (s + mod_pow(i, p)) % kMod
            y[p][i] = s
            
    fact = [1] * (max_p + 2)
    invfact = [1] * (max_p + 2)
    for i in range(1, len(fact)):
        fact[i] = (fact[i - 1] * i) % kMod
    invfact[-1] = mod_inv(fact[-1])
    for i in range(len(fact) - 1, 0, -1):
        invfact[i - 1] = (invfact[i] * i) % kMod
        
    return SumPowers(y, fact, invfact)

def sum_pows(p, n, sp):
    if n <= p + 1:
        return sp.y[p][n]
        
    m = p + 1
    pre = [1] * (m + 2)
    suf = [1] * (m + 2)
    nmod = n % kMod
    
    for i in range(m + 1):
        pre[i + 1] = (pre[i] * (nmod - i)) % kMod
    for i in range(m, -1, -1):
        suf[i] = (suf[i + 1] * (nmod - i)) % kMod
        
    res = 0
    for i in range(m + 1):
        num = (pre[i] * suf[i + 1]) % kMod
        denom_inv = (sp.invfact[i] * sp.invfact[m - i]) % kMod
        term = (sp.y[p][i] * num) % kMod
        term = (term * denom_inv) % kMod
        if (m - i) & 1:
            if term != 0: term = kMod - term
        res = (res + term) % kMod
    return res

def sum_poly(coeffs, n, sp):
    res = 0
    for p in range(len(coeffs)):
        if coeffs[p] == 0: continue
        spow = sum_pows(p, n, sp)
        res = (res + coeffs[p] * spow) % kMod
    return res

def build_pairs(n_val):
    U3 = [1, 6]
    while True:
        nex = 2 * 3 * U3[-1] - U3[-2]
        U3.append(nex)
        if nex > n_val: break
    
    max_idx = len(U3) - 1
    pairs = []
    for i in range(max_idx + 1):
        for j in range(i + 1, max_idx + 1):
            tmp = U3[i] * U3[j]
            if tmp <= n_val:
                pairs.append({'i_idx': i, 'j_idx': j})
    return pairs

def compute_mod(n_val):
    pairs = build_pairs(n_val)
    max_idx = 0
    k_max = 0
    for p in pairs:
        p['k_max'] = max_k_for_pair(n_val, p['i_idx'], p['j_idx'])
        k_max = max(k_max, p['k_max'])
        max_idx = max(max_idx, p['i_idx'], p['j_idx'])
        
    nf = generate_non_fundamentals(k_max)
    nf_set = set(nf)
    
    inv2 = mod_inv(2)
    U_polys = build_U_polys(max_idx)
    t_poly = [0, inv2, inv2]
    deg_max = 0
    for p in pairs:
        r = poly_mul(U_polys[p['i_idx']], U_polys[p['j_idx']])
        p['poly'] = poly_mul(r, t_poly)
        deg_max = max(deg_max, len(p['poly']) - 1)
        
    sp = precompute_sum_pows(deg_max)
    for p in pairs:
        p['sum_all'] = sum_poly(p['poly'], p['k_max'], sp)
        
    for p in pairs:
        p['sum_nf'] = 0
        for k in nf:
            if k > p['k_max']: continue
            k_mod = k % kMod
            x_mod = (2 * k_mod + 1) % kMod
            U = [0] * (max_idx + 1)
            U[0] = 1
            if max_idx >= 1: U[1] = (2 * x_mod) % kMod
            for m in range(1, max_idx):
                U[m + 1] = (2 * x_mod * U[m] - U[m - 1]) % kMod
            t = (k_mod * (k_mod + 1)) % kMod
            t = (t * inv2) % kMod
            val = (t * U[p['i_idx']]) % kMod
            val = (val * U[p['j_idx']]) % kMod
            p['sum_nf'] = (p['sum_nf'] + val) % kMod
            
    ans = 0
    for p in pairs:
        term = (p['sum_all'] - p['sum_nf']) % kMod
        ans = (ans + term) % kMod
        
    return ans

def solve():
    n = 10**35
    ans = compute_mod(n)
    return str(ans)

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

Java

import java.util.*;

public class Euler833 {

    static final long kMod = 136101521L;

    static class PairInfo {
        int iIdx;
        int jIdx;
        long kMax;
        long[] poly;
        long sumAll = 0;
        long sumNf = 0;
    }

    static boolean mulLeq(long a, long b, long limit) {
        if (a == 0 || b == 0)
            return true;
        if (a > limit / b)
            return false;
        return true;
    }

    static boolean evalPairLeq(long k, int iIdx, int jIdx, double limit) {
        // Since limit could be 10^35, we can't do exact long.
        // We'll use BigInteger maybe? Or simply approximate with double for upper
        // bounds.
        // wait, 10^35 fits in standard double ~1.7e308. Wait, exact arithmetic is
        // better.
        return true;
    }

    // Implementing exact 128 bit behavior with BigInteger
    static java.math.BigInteger toBI(long v) {
        return java.math.BigInteger.valueOf(v);
    }

    static boolean mulLeq(java.math.BigInteger a, java.math.BigInteger b, java.math.BigInteger limit,
            java.math.BigInteger[] out) {
        if (a.equals(java.math.BigInteger.ZERO) || b.equals(java.math.BigInteger.ZERO)) {
            out[0] = java.math.BigInteger.ZERO;
            return true;
        }
        java.math.BigInteger p = a.multiply(b);
        if (p.compareTo(limit) > 0)
            return false;
        out[0] = p;
        return true;
    }

    static boolean evalPairLeq(long k, int iIdx, int jIdx, java.math.BigInteger limit) {
        java.math.BigInteger kk = toBI(k);
        java.math.BigInteger x = kk.multiply(toBI(2)).add(java.math.BigInteger.ONE);
        int maxIdx = Math.max(iIdx, jIdx);
        java.math.BigInteger[] U = new java.math.BigInteger[maxIdx + 1];
        U[0] = java.math.BigInteger.ONE;
        if (maxIdx >= 1)
            U[1] = x.multiply(toBI(2));

        java.math.BigInteger twoX = x.multiply(toBI(2));
        for (int m = 1; m < maxIdx; ++m) {
            java.math.BigInteger thresh = limit.add(U[m - 1]).divide(twoX);
            if (U[m].compareTo(thresh) > 0) {
                U[m + 1] = limit.add(java.math.BigInteger.ONE);
            } else {
                java.math.BigInteger next = twoX.multiply(U[m]).subtract(U[m - 1]);
                U[m + 1] = next.compareTo(limit) > 0 ? limit.add(java.math.BigInteger.ONE) : next;
            }
        }

        java.math.BigInteger t = kk.multiply(kk.add(java.math.BigInteger.ONE)).divide(toBI(2));
        java.math.BigInteger[] tmp = new java.math.BigInteger[1];
        if (!mulLeq(t, U[iIdx], limit, tmp))
            return false;
        if (!mulLeq(tmp[0], U[jIdx], limit, tmp))
            return false;
        return true;
    }

    static long maxKForPair(java.math.BigInteger n, int iIdx, int jIdx) {
        long lo = 0;
        long hi = 1;
        while (evalPairLeq(hi, iIdx, jIdx, n)) {
            if (hi > (Long.MAX_VALUE / 2))
                break;
            hi *= 2;
        }
        while (lo + 1 < hi) {
            long mid = lo + (hi - lo) / 2;
            if (evalPairLeq(mid, iIdx, jIdx, n)) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        return lo;
    }

    static ArrayList<Long> generateNonFundamentals(long kMax) {
        ArrayList<Long> nf = new ArrayList<>();
        if (kMax < 4)
            return nf;
        long k1Max = (long) Math.sqrt((double) kMax / 4.0) + 2;
        for (long k1 = 1; k1 <= k1Max; ++k1) {
            java.math.BigInteger x1 = toBI(2 * k1 + 1);
            java.math.BigInteger xPrev = java.math.BigInteger.ONE;
            java.math.BigInteger xCur = x1;
            while (true) {
                java.math.BigInteger xNext = toBI(2).multiply(x1).multiply(xCur).subtract(xPrev);
                java.math.BigInteger kNextBI = xNext.subtract(java.math.BigInteger.ONE).divide(toBI(2));
                if (kNextBI.compareTo(toBI(kMax)) > 0)
                    break;
                nf.add(kNextBI.longValue());
                xPrev = xCur;
                xCur = xNext;
            }
        }
        Collections.sort(nf);
        ArrayList<Long> distinctNf = new ArrayList<>();
        for (long v : nf) {
            if (distinctNf.isEmpty() || distinctNf.get(distinctNf.size() - 1) != v) {
                distinctNf.add(v);
            }
        }
        return distinctNf;
    }

    static long modPow(long base, long exp) {
        long res = 1 % kMod;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                res = (res * base) % kMod;
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return res;
    }

    static long modInv(long a) {
        return modPow(a, kMod - 2);
    }

    static long[] polyScale(long[] a, long s) {
        long[] c = new long[a.length];
        for (int i = 0; i < a.length; ++i) {
            c[i] = (a[i] * s) % kMod;
        }
        return c;
    }

    static long[] polySub(long[] a, long[] b) {
        int n = Math.max(a.length, b.length);
        long[] c = new long[n];
        for (int i = 0; i < n; ++i) {
            long va = i < a.length ? a[i] : 0;
            long vb = i < b.length ? b[i] : 0;
            c[i] = (va - vb) % kMod;
            if (c[i] < 0)
                c[i] += kMod;
        }
        return c;
    }

    static long[] polyMul(long[] a, long[] b) {
        long[] c = new long[a.length + b.length - 1];
        for (int i = 0; i < a.length; ++i) {
            for (int j = 0; j < b.length; ++j) {
                c[i + j] = (c[i + j] + a[i] * b[j]) % kMod;
            }
        }
        return c;
    }

    static long[][] buildUPolys(int maxIdx) {
        long[] t = { 1, 2 };
        long[][] U = new long[maxIdx + 1][];
        U[0] = new long[] { 1 };
        if (maxIdx >= 1)
            U[1] = polyScale(t, 2);

        for (int m = 1; m < maxIdx; ++m) {
            long[] temp = polyMul(t, U[m]);
            temp = polyScale(temp, 2);
            U[m + 1] = polySub(temp, U[m - 1]);
        }
        return U;
    }

    static class SumPowers {
        long[][] y;
        long[] fact;
        long[] invfact;
    }

    static SumPowers precomputeSumPows(int maxP) {
        SumPowers sp = new SumPowers();
        sp.y = new long[maxP + 1][maxP + 2];
        for (int p = 0; p <= maxP; ++p) {
            int m = p + 1;
            long s = 0;
            for (int i = 1; i <= m; ++i) {
                s = (s + modPow(i, p)) % kMod;
                sp.y[p][i] = s;
            }
        }
        sp.fact = new long[maxP + 2];
        sp.invfact = new long[maxP + 2];
        sp.fact[0] = 1;
        sp.invfact[0] = 1;
        for (int i = 1; i < sp.fact.length; ++i) {
            sp.fact[i] = (sp.fact[i - 1] * i) % kMod;
        }
        sp.invfact[sp.invfact.length - 1] = modInv(sp.fact[sp.fact.length - 1]);
        for (int i = sp.fact.length - 1; i >= 1; --i) {
            sp.invfact[i - 1] = (sp.invfact[i] * i) % kMod;
        }
        return sp;
    }

    static long sumPows(int p, long n, SumPowers sp) {
        if (n <= p + 1) {
            return sp.y[p][(int) n];
        }
        int m = p + 1;
        long[] pre = new long[m + 2];
        long[] suf = new long[m + 2];
        pre[0] = 1;
        suf[m + 1] = 1;
        long nmod = n % kMod;

        for (int i = 0; i <= m; ++i) {
            long term = (nmod - i) % kMod;
            if (term < 0)
                term += kMod;
            pre[i + 1] = (pre[i] * term) % kMod;
        }
        for (int i = m; i >= 0; --i) {
            long term = (nmod - i) % kMod;
            if (term < 0)
                term += kMod;
            suf[i] = (suf[i + 1] * term) % kMod;
        }

        long res = 0;
        for (int i = 0; i <= m; ++i) {
            long num = (pre[i] * suf[i + 1]) % kMod;
            long denomInv = (sp.invfact[i] * sp.invfact[m - i]) % kMod;
            long term = (sp.y[p][i] * num) % kMod;
            term = (term * denomInv) % kMod;
            if ((m - i) % 2 != 0) {
                if (term != 0)
                    term = kMod - term;
            }
            res = (res + term) % kMod;
        }
        return res;
    }

    static long sumPoly(long[] coeffs, long n, SumPowers sp) {
        long res = 0;
        for (int p = 0; p < coeffs.length; ++p) {
            if (coeffs[p] == 0)
                continue;
            long spow = sumPows(p, n, sp);
            res = (res + coeffs[p] * spow) % kMod;
        }
        return res;
    }

    public static String solve() {
        java.math.BigInteger n = java.math.BigInteger.TEN.pow(35);

        ArrayList<java.math.BigInteger> U3 = new ArrayList<>();
        U3.add(java.math.BigInteger.ONE);
        U3.add(toBI(6));
        while (true) {
            java.math.BigInteger next = toBI(6).multiply(U3.get(U3.size() - 1)).subtract(U3.get(U3.size() - 2));
            U3.add(next);
            if (next.compareTo(n) > 0)
                break;
        }

        int maxIdx = U3.size() - 1;
        ArrayList<PairInfo> pairs = new ArrayList<>();

        for (int i = 0; i <= maxIdx; ++i) {
            for (int j = i + 1; j <= maxIdx; ++j) {
                java.math.BigInteger[] tmp = new java.math.BigInteger[1];
                if (!mulLeq(U3.get(i), U3.get(j), n, tmp))
                    continue;
                PairInfo info = new PairInfo();
                info.iIdx = i;
                info.jIdx = j;
                pairs.add(info);
            }
        }

        long maxK = 0;
        for (PairInfo p : pairs) {
            p.kMax = maxKForPair(n, p.iIdx, p.jIdx);
            maxK = Math.max(maxK, p.kMax);
        }

        ArrayList<Long> nf = generateNonFundamentals(maxK);

        long inv2 = modInv(2);
        long[][] UPolys = buildUPolys(maxIdx);
        long[] tPoly = { 0, inv2, inv2 };
        int degMax = 0;

        for (PairInfo p : pairs) {
            long[] r = polyMul(UPolys[p.iIdx], UPolys[p.jIdx]);
            p.poly = polyMul(r, tPoly);
            degMax = Math.max(degMax, p.poly.length - 1);
        }

        SumPowers sp = precomputeSumPows(degMax);
        for (PairInfo p : pairs) {
            p.sumAll = sumPoly(p.poly, p.kMax, sp);
        }

        for (PairInfo p : pairs) {
            p.sumNf = 0;
            for (long k : nf) {
                if (k > p.kMax)
                    continue;
                long kModV = k % kMod;
                long xMod = (2 * kModV + 1) % kMod;
                long[] U = new long[maxIdx + 1];
                U[0] = 1;
                if (maxIdx >= 1)
                    U[1] = (2 * xMod) % kMod;
                for (int m = 1; m < maxIdx; ++m) {
                    U[m + 1] = (2 * xMod * U[m] - U[m - 1]) % kMod;
                    if (U[m + 1] < 0)
                        U[m + 1] += kMod;
                }
                long t = (kModV * ((kModV + 1) % kMod)) % kMod;
                t = (t * inv2) % kMod;

                long val = (t * U[p.iIdx]) % kMod;
                val = (val * U[p.jIdx]) % kMod;

                p.sumNf = (p.sumNf + val) % kMod;
            }
        }

        long ans = 0;
        for (PairInfo p : pairs) {
            long term = (p.sumAll - p.sumNf) % kMod;
            if (term < 0)
                term += kMod;
            ans = (ans + term) % kMod;
        }

        return Long.toString(ans);
    }

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