Problem 721: High Powers of Irrational Numbers

Mathematical Foundation

Theorem 1 (Companion Integer Sequence). Let $a$ be a positive integer that is not a perfect square and set $m = \lceil\sqrt{a}\rceil$. Define $\alpha = m + \sqrt{a}$ and $\beta = m - \sqrt{a}$. Then the sequence $T_n = \alpha^n + \beta^n$ satisfies the linear recurrence

and $T_n \in \mathbb{Z}$ for all $n \ge 0$.

Proof. The values $\alpha$ and $\beta$ are the two roots of the quadratic

By Newton’s identity for power sums of roots, $T_n = (\alpha+\beta)T_{n-1} - \alpha\beta\,T_{n-2} = 2m\,T_{n-1} - (m^2-a)\,T_{n-2}$. The initial conditions are $T_0 = \alpha^0 + \beta^0 = 2$ and $T_1 = \alpha + \beta = 2m$. Since $2m$ and $m^2 - a$ are integers and $T_0, T_1 \in \mathbb{Z}$, induction on $n$ gives $T_n \in \mathbb{Z}$ for all $n \ge 0$. $\square$

Lemma 1 (Floor Extraction). For non-perfect-square $a$ and $n \ge 1$, we have $f(a,n) = T_n - 1$.

Proof. Since $a$ is not a perfect square, $m = \lceil\sqrt{a}\rceil$ satisfies $m > \sqrt{a}$, so $\beta = m - \sqrt{a} > 0$. Moreover, $m \le \sqrt{a} + 1$ (as $m$ is the least integer $\ge \sqrt{a}$), which gives $\beta = m - \sqrt{a} \le 1$. The inequality is strict because $a$ is not a perfect square, so $\sqrt{a} \notin \mathbb{Z}$ and thus $m - \sqrt{a} < 1$. Hence $0 < \beta < 1$, and therefore $0 < \beta^n < 1$ for all $n \ge 1$.

Now $\alpha^n = T_n - \beta^n$ where $T_n$ is a positive integer and $0 < \beta^n < 1$. Therefore $T_n - 1 < \alpha^n < T_n$, which gives $\lfloor \alpha^n \rfloor = T_n - 1$. $\square$

Lemma 2 (Matrix Recurrence). The recurrence $T_n = 2m\,T_{n-1} - c\,T_{n-2}$ (where $c = m^2 - a$) admits the matrix form

Proof. Direct verification: the matrix equation encodes $T_n = 2m\,T_{n-1} - c\,T_{n-2}$ in the first row, and $T_{n-1} = T_{n-1}$ in the second row. Induction on $n$ yields the stated power formula. $\square$

Editorial

Compute G(N) = sum_{a=1}^{N} f(a, a^2) mod 999999937 where f(a, n) = floor((ceil(sqrt(a)) + sqrt(a))^n) and a ranges over non-perfect-squares. Key insight: alpha = ceil(sqrt(a)) + sqrt(a), beta = ceil(sqrt(a)) - sqrt(a) satisfy a quadratic, so alpha^n + beta^n is an integer (via linear recurrence). Since 0 < beta < 1, f(a,n) = alpha^n + beta^n - 1.

Pseudocode

Compute T_{a^2} mod p via matrix exponentiation
if exp is odd

Complexity Analysis

• Time: For each non-perfect-square $a \le N$, matrix exponentiation computes $M^{a^2-1}$ using $O(\log(a^2)) = O(\log a)$ multiplications of $2 \times 2$ matrices over $\mathbb{Z}/p\mathbb{Z}$. Each multiplication is $O(1)$. There are $N - O(\sqrt{N})$ non-perfect-square values. Total: $O(N \log N)$.
• Space: $O(1)$ auxiliary space per value of $a$ (two $2\times 2$ matrices). Overall $O(\sqrt{N})$ for the perfect-square set, or $O(1)$ if squares are detected on the fly.

Answer

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

/*
 * Problem 721: High Powers of Irrational Numbers
 *
 * For non-perfect-square a, f(a,n) = floor((ceil(sqrt(a)) + sqrt(a))^n).
 * Using the companion algebraic integer beta = ceil(sqrt(a)) - sqrt(a),
 * we have alpha^n + beta^n = T_n (integer), and f(a,n) = T_n - 1.
 *
 * T_n satisfies: T_k = 2m * T_{k-1} - (m^2-a) * T_{k-2}
 * where m = ceil(sqrt(a)).
 *
 * We compute T_{a^2} mod p via 2x2 matrix exponentiation.
 */

const long long MOD = 999999937LL;

typedef vector<vector<long long>> Mat;

Mat mat_mul(const Mat& A, const Mat& B, long long mod) {
    int n = A.size();
    Mat C(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++) if (A[i][k])
            for (int j = 0; j < n; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
    return C;
}

Mat mat_pow(Mat M, long long p, long long mod) {
    int n = M.size();
    Mat result(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++) result[i][i] = 1;
    while (p > 0) {
        if (p & 1) result = mat_mul(result, M, mod);
        M = mat_mul(M, M, mod);
        p >>= 1;
    }
    return result;
}

long long compute_T(long long m, long long c, long long n, long long mod) {
    // T_k = 2m * T_{k-1} - c * T_{k-2}, T_0=2, T_1=2m
    if (n == 0) return 2 % mod;
    if (n == 1) return (2 * m) % mod;
    Mat M = {{(2*m) % mod, (mod - c % mod) % mod}, {1, 0}};
    Mat R = mat_pow(M, n - 1, mod);
    return (R[0][0] * (2*m % mod) % mod + R[0][1] * 2 % mod) % mod;
}

bool is_perfect_square(long long a) {
    long long s = (long long)sqrt((double)a);
    while (s * s > a) s--;
    while ((s+1)*(s+1) <= a) s++;
    return s * s == a;
}

int main() {
    long long N = 1000; // Change to 5000000 for full answer
    long long total = 0;

    for (long long a = 1; a <= N; a++) {
        if (is_perfect_square(a)) continue;
        long long s = (long long)sqrt((double)a);
        while (s * s > a) s--;
        long long m = s + 1;
        long long c = m * m - a;
        long long n = a * a;
        long long T = compute_T(m, c, n, MOD);
        total = (total + T - 1 + MOD) % MOD;
    }

    cout << total << endl;
    return 0;
}

"""
Problem 721: High Powers of Irrational Numbers

Compute G(N) = sum_{a=1}^{N} f(a, a^2) mod 999999937
where f(a, n) = floor((ceil(sqrt(a)) + sqrt(a))^n)
and a ranges over non-perfect-squares.

Key insight: alpha = ceil(sqrt(a)) + sqrt(a), beta = ceil(sqrt(a)) - sqrt(a)
satisfy a quadratic, so alpha^n + beta^n is an integer (via linear recurrence).
Since 0 < beta < 1, f(a,n) = alpha^n + beta^n - 1.
"""

from math import isqrt

MOD = 999_999_937

def mat_mul(A, B, mod):
    """Multiply two 2x2 matrices modulo mod."""
    return [
        [(A[0][0]*B[0][0] + A[0][1]*B[1][0]) % mod,
         (A[0][0]*B[0][1] + A[0][1]*B[1][1]) % mod],
        [(A[1][0]*B[0][0] + A[1][1]*B[1][0]) % mod,
         (A[1][0]*B[0][1] + A[1][1]*B[1][1]) % mod],
    ]

def mat_pow(M, n, mod):
    """Compute M^n mod p for a 2x2 matrix."""
    result = [[1, 0], [0, 1]]  # identity
    base = [row[:] for row in M]
    while n > 0:
        if n & 1:
            result = mat_mul(result, base, mod)
        base = mat_mul(base, base, mod)
        n >>= 1
    return result

def T_n(m, c, n, mod):
    """Compute T_n = alpha^n + beta^n mod p using matrix exponentiation.
    Recurrence: T_k = 2m * T_{k-1} - c * T_{k-2}, T_0=2, T_1=2m.
    Here c = m^2 - a.
    """
    if n == 0:
        return 2 % mod
    if n == 1:
        return (2 * m) % mod
    M = [[(2*m) % mod, (-c) % mod], [1, 0]]
    R = mat_pow(M, n - 1, mod)
    # T_n = R[0][0]*T_1 + R[0][1]*T_0
    return (R[0][0] * (2*m) + R[0][1] * 2) % mod

def is_perfect_square(a):
    s = isqrt(a)
    return s * s == a

# --- Method 1: Direct computation for small N ---
def solve_small(N):
    """Compute G(N) mod p for small N (brute force verification)."""
    total = 0
    for a in range(1, N + 1):
        if is_perfect_square(a):
            continue
        s = isqrt(a)
        m = s + 1  # ceil(sqrt(a))
        c = m * m - a
        n = a * a
        val = T_n(m, c, n, MOD)
        total = (total + val - 1) % MOD
    return total

# --- Method 2: Same algorithm, for full N ---
def solve(N):
    """Compute G(N) mod p."""
    return solve_small(N)

# --- Verification ---
# f(5, 2): m=3, c=4, T_2 = 28, f=27
m, c = 3, 4
assert T_n(m, c, 2, MOD) == 28
assert T_n(m, c, 2, MOD) - 1 == 27  # f(5,2) = 27

# f(5, 5): T_5 = 3936, f=3935
assert T_n(m, c, 5, MOD) == 3936
assert T_n(m, c, 5, MOD) - 1 == 3935  # f(5,5) = 3935

# G(1000) = 163861845
result_1000 = solve(1000)
assert result_1000 == 163861845, f"G(1000) = {result_1000}, expected 163861845"

print(f"G(1000) = {result_1000}")

# Full answer would be solve(5000000) -- takes a few minutes
# ans = solve(5000000)
# print(ans)