Problem 611: Hallway of Square Steps

Project Euler Problem 611 Solution

EulerSolve provides an optimized solution for Project Euler Problem 611, Hallway of Square Steps, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every pair of distinct positive squares $a^2 \lt b^2$, Peter toggles door $n=a^2+b^2$, provided $n\le N$. After all actions, door $n$ remains open exactly when the number of strict representations $$u(n)=\#\left\{(a,b)\in \mathbb{Z}_{>0}^2 : a \lt b,\ a^2+b^2=n\right\}$$ is odd. Therefore $$F(N)=\#\left\{n\le N : u(n)\equiv 1 \pmod 2\right\}.$$ The real input is $N=10^{12}$, so enumerating all square pairs is far too slow. The implementation instead turns the problem into a classification of prime exponents. Mathematical Approach The central object is the classical sum-of-two-squares counting function over all signs and orders. Step 1: Relate Door Toggles to the Classical Representation Count Let $$r_2(n)=\#\left\{(x,y)\in \mathbb{Z}^2 : x^2+y^2=n\right\}.$$ Every strict positive representation $a \lt b$ generates eight ordered signed solutions: $(\pm a,\pm b)$ and $(\pm b,\pm a)$. Two degenerate situations contribute only four solutions: $$n=c^2 \quad \text{gives} \quad (\pm c,0),(0,\pm c),$$ $$n=2c^2 \quad \text{gives} \quad (\pm c,\pm c).$$ So if we define $$s(n)= \begin{cases} 1,&\text{if } n \text{ is a square or twice a square},\\ 0,&\text{otherwise}, \end{cases}$$ then $$r_2(n)=8u(n)+4s(n),$$ hence $$\frac{r_2(n)}{4}=2u(n)+s(n).$$ This identity is the bridge from door toggling to arithmetic parity....

Detailed mathematical approach

Problem Summary

For every pair of distinct positive squares $a^2 \lt b^2$, Peter toggles door $n=a^2+b^2$, provided $n\le N$. After all actions, door $n$ remains open exactly when the number of strict representations

$$u(n)=\#\left\{(a,b)\in \mathbb{Z}_{>0}^2 : a \lt b,\ a^2+b^2=n\right\}$$

is odd. Therefore

$$F(N)=\#\left\{n\le N : u(n)\equiv 1 \pmod 2\right\}.$$

The real input is $N=10^{12}$, so enumerating all square pairs is far too slow. The implementation instead turns the problem into a classification of prime exponents.

Mathematical Approach

The central object is the classical sum-of-two-squares counting function over all signs and orders.

Step 1: Relate Door Toggles to the Classical Representation Count

Let

$$r_2(n)=\#\left\{(x,y)\in \mathbb{Z}^2 : x^2+y^2=n\right\}.$$

Every strict positive representation $a \lt b$ generates eight ordered signed solutions: $(\pm a,\pm b)$ and $(\pm b,\pm a)$. Two degenerate situations contribute only four solutions:

$$n=c^2 \quad \text{gives} \quad (\pm c,0),(0,\pm c),$$

$$n=2c^2 \quad \text{gives} \quad (\pm c,\pm c).$$

So if we define

$$s(n)= \begin{cases} 1,&\text{if } n \text{ is a square or twice a square},\\ 0,&\text{otherwise}, \end{cases}$$

then

$$r_2(n)=8u(n)+4s(n),$$

hence

$$\frac{r_2(n)}{4}=2u(n)+s(n).$$

This identity is the bridge from door toggling to arithmetic parity.

Step 2: Use the Sum-of-Two-Squares Theorem

Write the factorization of $n$ as

$$n=2^a\prod_{p\equiv 1 \pmod 4} p^{\alpha_p}\prod_{q\equiv 3 \pmod 4} q^{\beta_q}.$$

If any prime $q\equiv 3 \pmod 4$ has odd exponent, then $n$ is not representable as a sum of two squares at all, so $r_2(n)=0$ and $u(n)$ is even. The only interesting numbers are therefore

$$n=2^a\prod_{p\equiv 1 \pmod 4} p^{\alpha_p}\prod_{q\equiv 3 \pmod 4} q^{2\gamma_q}.$$

For such $n$, the standard formula becomes

$$r_2(n)=4\prod_{p\equiv 1 \pmod 4}(\alpha_p+1).$$

Define

$$R(n)=\prod_{p\equiv 1 \pmod 4}(\alpha_p+1).$$

Then the parity condition is simply

$$R(n)=2u(n)+s(n).$$

So the question becomes: for which exponent patterns is $(R(n)-s(n))/2$ odd?

Step 3: Family A, the Non-Degenerate Case

First assume $s(n)=0$, meaning $n$ is neither a square nor twice a square. Then $u(n)$ is odd exactly when

$$R(n)\equiv 2 \pmod 4.$$

Because each factor $\alpha_p+1$ is odd when $\alpha_p$ is even and even when $\alpha_p$ is odd, the product can be $2 \pmod 4$ only when exactly one factor contributes a single power of $2$ and every other factor remains odd. So there must be exactly one prime $p\equiv 1 \pmod 4$ with odd exponent, and that exponent must satisfy

$$\alpha_p+1\equiv 2 \pmod 4 \iff \alpha_p\equiv 1 \pmod 4.$$

Therefore the non-degenerate doors that remain open are exactly those of the form

$$n=p^{4t+1}m^2 \quad \text{or} \quad n=2p^{4t+1}m^2,$$

where $p\equiv 1 \pmod 4$ and $p\nmid m$. This is the first family counted by the implementation.

Step 4: Family B, the Square / Twice-Square Correction

Now assume $s(n)=1$. That means every exponent $\alpha_p$ is even, so $n$ is either a square or twice a square. We can write it uniquely as

$$n=2^a m^2,\qquad m \text{ odd}.$$

Since $u(n)$ is odd exactly when

$$R(n)\equiv 3 \pmod 4,$$

we inspect each factor with $\alpha_p=2\delta_p$. Then

$$\alpha_p+1=2\delta_p+1\equiv \begin{cases} 1 \pmod 4,&\delta_p \text{ even},\\ 3 \pmod 4,&\delta_p \text{ odd}. \end{cases}$$

So $R(n)\equiv 3 \pmod 4$ exactly when an odd number of primes $p\equiv 1 \pmod 4$ occur to odd exponent in the odd square root $m$. That produces the second family:

$$n=2^a m^2,\qquad m \text{ odd},$$

with the parity filter “the number of primes $p\equiv 1 \pmod 4$ occurring to odd exponent in $m$ is odd.”

Step 5: Count Family A Efficiently

Let $\pi_1(x)$ denote the number of primes $p\le x$ with $p\equiv 1 \pmod 4$. For the even-$2$-adic half of family A, we count numbers

$$n=p^{4t+1}m^2,\qquad p\equiv 1 \pmod 4,\ p\nmid m,\ n\le N.$$

For the first exponent layer $p^1$, the contribution of one prime is

$$\#\{m : p m^2\le N,\ p\nmid m\}=\left\lfloor\sqrt{\frac{N}{p}}\right\rfloor-\left\lfloor\sqrt{\frac{N}{p^3}}\right\rfloor.$$

The first term is grouped by equal square-root quotients:

$$\sum_{p\equiv 1 \pmod 4}\left\lfloor\sqrt{\frac{N}{p}}\right\rfloor =\sum_{t\ge 1} t\left(\pi_1\left(\left\lfloor\frac{N}{t^2}\right\rfloor\right)-\pi_1\left(\left\lfloor\frac{N}{(t+1)^2}\right\rfloor\right)\right).$$

The subtraction by $\left\lfloor\sqrt{N/p^3}\right\rfloor$ removes the forbidden choices where the square part already contains $p$. Higher valid exponents $p^5,p^9,\dots$ are sparse, so they are added directly through

$$\left\lfloor\sqrt{\frac{N}{p^{4t+1}}}\right\rfloor-\left\lfloor\sqrt{\frac{N}{p^{4t+3}}}\right\rfloor,\qquad t\ge 1.$$

If $v_2(n)$ is odd, then $n$ belongs to the second variant $2p^{4t+1}m^2$. Multiplying by $2$ gives a bijection with the even-$v_2$ case at bound $N/2$, so the same routine can be reused with $N$ replaced by $N/2$.

Step 6: Count Family B Efficiently

For family B, every valid odd $m\le \sqrt{N}$ contributes all numbers

$$m^2,\ 2m^2,\ 4m^2,\ \dots,\ 2^a m^2\le N.$$

The number of admissible powers of $2$ is

$$\left\lfloor \log_2\left(\frac{N}{m^2}\right)\right\rfloor+1.$$

So the remaining work is to know, for each odd $m$, whether the count of primes $p\equiv 1 \pmod 4$ appearing to odd exponent is odd or even.

Worked Example: $N=100$

The checkpoint $F(100)=27$ follows cleanly from the classification.

Family A with even $v_2$ gives

$$5,13,17,20,29,37,41,45,52,53,61,68,73,80,89,97,$$

so there are $16$ such doors.

Family A with odd $v_2$ is obtained by doubling the even-$v_2$ doors up to $50$, namely

$$10,26,34,40,58,74,82,90,$$

which contributes $8$ more doors.

For family B, the odd candidates $m\le 10$ are $1,3,5,7,9$. Only $m=5$ has an odd number of $1 \pmod 4$ primes occurring to odd exponent, so it contributes

$$25,\ 50,\ 100,$$

hence $3$ more doors. Therefore

$$F(100)=16+8+3=27.$$

How the Code Works

The C++, Python, and Java implementations all follow the same arithmetic decomposition. The Python entry point is only a thin wrapper around that same numerical strategy, so the mathematics is identical across languages.

The implementation first builds a prime-counting structure over the quotient values $\lfloor N/i\rfloor$. It stores both the ordinary prime count $\pi(x)$ and a weighted prime sum using the Dirichlet character

$$\chi_4(n)= \begin{cases} 0,&n \text{ even},\\ 1,&n\equiv 1 \pmod 4,\\ -1,&n\equiv 3 \pmod 4. \end{cases}$$

From these two tables it recovers

$$\pi_1(x)=\frac{(\pi(x)-1)+\sum_{p\le x}\chi_4(p)}{2},$$

which is exactly the count of primes $p\le x$ with $p\equiv 1 \pmod 4$.

With that tool in place, the implementation counts family A in three pieces: the grouped $p^1$ layer, the correction that subtracts roots already divisible by the chosen prime, and the direct enumeration of the sparse higher layers $p^5,p^9,\dots$. It then reuses the same routine at $N/2$ to count the odd-$v_2$ branch.

For family B, the implementation precomputes smallest prime factors up to $\lfloor\sqrt{N}\rfloor$. Factoring each odd $m$ reveals whether the number of primes $p\equiv 1 \pmod 4$ occurring to odd exponent is odd. Whenever the parity condition holds, it adds

$$\left\lfloor \log_2\left(\frac{N}{m^2}\right)\right\rfloor+1$$

to the answer.

Finally, the implementation adds the even-$v_2$ branch of family A, the corresponding odd-$v_2$ branch obtained from the $N/2$ bijection, and the family-B contribution.

Complexity Analysis

Let $M=\lfloor\sqrt{N}\rfloor$. The smallest-prime-factor sieve and the parity table for odd $m$ cost $O(M\log\log M)$ time and $O(M)$ memory. The family-B scan is linear in $M$. The higher-power part of family A is very small because $p^5$ already grows quickly.

The dominant work is the quotient-class prime-counting structure used to answer many values of $\pi_1(x)$ efficiently. It stores $O(M)$ quotient values and keeps the overall method practical for $N=10^{12}$. So the implementation uses $O(M)$ memory, and the runtime is dominated by the prime-counting preprocessing plus the linear-in-$M$ auxiliary passes.

Footnotes and References

Problem page: https://projecteuler.net/problem=611
Fermat's theorem on sums of two squares: Wikipedia — Fermat's theorem on sums of two squares
Sum of two squares function: Wikipedia — Sum of two squares function
Dirichlet character: Wikipedia — Dirichlet character
Prime-counting function: Wikipedia — Prime-counting function

Problem 611 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

// Project Euler 611: count n <= N where #{0<a<b : a^2 + b^2 = n} is odd.
// Reduce to two arithmetic families via r2(n) and prime-exponent parity; count with a prime-counting sieve for p ≡ 1 (mod 4).

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

static u64 isqrt_u64(u64 x) {
    u64 r = (u64)std::sqrt((long double)x);
    while ((u128)(r + 1) * (r + 1) <= x) ++r;
    while ((u128)r * r > x) --r;
    return r;
}

static int ilog2_u64(u64 x) {
    int r = 0;
    while ((1ULL << (r + 1)) <= x) ++r;
    return r;
}

static std::vector<int> sieve_primes_int(int n) {
    std::vector<bool> is_prime(n + 1, true);
    if (n >= 0) is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;
    for (int p = 2; (u64)p * p <= (u64)n; ++p) {
        if (!is_prime[p]) continue;
        for (u64 j = (u64)p * p; j <= (u64)n; j += (u64)p) is_prime[(size_t)j] = false;
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

struct PrimeCountMod4 {
    u64 n = 0;
    u64 sq = 0;
    std::vector<u64> vals;
    std::vector<u64> g_pi;
    std::vector<i64> g_chi;
    std::vector<int> id1;
    std::vector<int> id2;
    std::vector<int> primes;
    std::vector<i64> pref_chi_prime;

    static i64 chi_int(u64 x) {
        if ((x & 1ULL) == 0) return 0;
        return (x & 3ULL) == 1 ? 1 : -1;
    }

    static i64 sum_chi_1_to(u64 m) {
        // chi(k)=0 for even; for odd k, +1 if 1 mod4 else -1.
        const u64 c1 = (m + 3) / 4;
        const u64 c3 = (m + 1) / 4;
        return (i64)c1 - (i64)c3;
    }

    explicit PrimeCountMod4(u64 n_) : n(n_) {
        sq = (u64)std::sqrt((long double)n);
        primes = sieve_primes_int((int)sq);

        for (u64 l = 1; l <= n;) {
            const u64 w = n / l;
            vals.push_back(w);
            l = n / w + 1;
        }
        const int m = (int)vals.size();

        g_pi.resize(m);
        g_chi.resize(m);
        id1.assign((size_t)sq + 1, -1);
        id2.assign((size_t)sq + 1, -1);

        for (int i = 0; i < m; ++i) {
            const u64 w = vals[i];
            g_pi[i] = (w >= 2) ? (w - 1) : 0; // count of integers in [2..w]
            i64 s = (w >= 1) ? (sum_chi_1_to(w) - 1) : 0; // sum_{k=2..w} chi(k)
            g_chi[i] = s;
            if (w <= sq) id1[w] = i;
            else id2[n / w] = i;
        }

        pref_chi_prime.assign(primes.size() + 1, 0);
        for (std::size_t i = 0; i < primes.size(); ++i) {
            const int p = primes[i];
            pref_chi_prime[i + 1] = pref_chi_prime[i] + chi_int((u64)p);
        }

        for (std::size_t i = 0; i < primes.size(); ++i) {
            const u64 p = (u64)primes[i];
            const u64 p2 = p * p;
            if (p2 > n) break;
            const i64 chi_p = chi_int(p);

            for (int j = 0; j < m && vals[j] >= p2; ++j) {
                const u64 w = vals[j];
                const int idx = id(w / p);
                g_pi[j] -= g_pi[idx] - (u64)i;
                if (chi_p != 0) g_chi[j] -= chi_p * (g_chi[idx] - pref_chi_prime[i]);
            }
        }
    }

    inline int id(u64 x) const { return (x <= sq) ? id1[x] : id2[n / x]; }

    inline u64 pi(u64 x) const { return g_pi[id(x)]; }

    inline i64 sum_chi_primes(u64 x) const { return g_chi[id(x)]; }

    inline u64 pi1(u64 x) const {
        if (x < 5) return 0;
        const u64 pix = pi(x);
        const u64 odd_primes = pix - 1; // exclude prime 2
        const i64 s = sum_chi_primes(x);
        return (u64)((odd_primes + s) / 2);
    }
};

static std::vector<int> sieve_spf(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);
    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const long long v = 1LL * p * i;
            if (v > n) break;
            spf[(int)v] = p;
            if (p == spf[i]) break;
        }
    }
    return spf;
}

static u64 count_B(u64 N, const std::vector<uint8_t>& parity_odd1mod4) {
    const u64 lim = isqrt_u64(N);
    u64 ans = 0;
    // Unique parameterization for squares / twice-squares: n = 2^a * m^2 with m odd.
    for (u64 m = 1; m <= lim; m += 2) {
        if (!parity_odd1mod4[(size_t)m]) continue;
        const u64 base = m * m;
        const u64 t = N / base;
        ans += (u64)(ilog2_u64(t) + 1);
    }
    return ans;
}

static u64 count_A_even_v2(u64 N, const PrimeCountMod4& pc, const std::vector<int>& primes_small) {
    // a=0 term: sum_{p≡1(4)} (T - T/p), T=floor(sqrt(N/p)).
    u128 S1 = 0;
    const u64 tmax = isqrt_u64(N);
    for (u64 t = 1; t <= tmax; ++t) {
        const u64 R = N / (t * t);
        if (R < 5) break;
        const u64 L = N / ((t + 1) * (t + 1));
        const u64 cnt = pc.pi1(R) - pc.pi1(L);
        S1 += (u128)t * cnt;
    }

    u128 corr = 0;
    const u64 cbrtN = 10000; // for N up to 1e12 (exact)
    for (int p : primes_small) {
        if ((u64)p > cbrtN) break;
        if (p % 4 != 1) continue;
        const u64 T = isqrt_u64(N / (u64)p);
        corr += (u64)(T / (u64)p);
    }

    u128 total = S1 - corr;

    // higher exponents: p^{4a+1} with a>=1 (start at p^5)
    for (int p : primes_small) {
        if (p % 4 != 1) continue;
        u128 pow = 1;
        for (int i = 0; i < 5; ++i) pow *= (u64)p;
        if (pow > N) break;
        const u128 p4 = (u128)p * p * p * p;
        while (pow <= N) {
            const u64 T = isqrt_u64((u64)(N / (u64)pow));
            total += (u128)(T - T / (u64)p);
            if (pow > (u128)N / p4) break;
            pow *= p4;
        }
    }

    return (u64)total;
}

static u64 F(u64 N, const PrimeCountMod4& pc, const std::vector<int>& primes_small, const std::vector<int>& spf,
             const std::vector<uint8_t>& parity_odd1mod4) {
    (void)spf;
    const u64 a_even = count_A_even_v2(N, pc, primes_small);
    const u64 a_odd = count_A_even_v2(N / 2, pc, primes_small); // multiply-by-2 bijection
    return a_even + a_odd + count_B(N, parity_odd1mod4);
}

int main() {
    const u64 Nmax = 1000000000000ULL;
    PrimeCountMod4 pc(Nmax);

    const int LIM = 1000000;
    const std::vector<int> spf = sieve_spf(LIM);
    const std::vector<int> primes_small = sieve_primes_int(LIM);

    // parity[m] = 1 iff m has an odd number of primes ≡1 (mod 4) to odd exponent.
    std::vector<uint8_t> parity((size_t)LIM + 1, 0);
    parity[1] = 0;
    for (int x = 2; x <= LIM; ++x) {
        int n = x;
        uint8_t par = 0;
        while (n > 1) {
            const int p = spf[n];
            int e = 0;
            while (n % p == 0) {
                n /= p;
                e ^= 1;
            }
            if (e && (p % 4 == 1)) par ^= 1;
        }
        parity[x] = par;
    }

    // Statement validations.
    assert(F(5, pc, primes_small, spf, parity) == 1ULL);
    assert(F(100, pc, primes_small, spf, parity) == 27ULL);
    assert(F(1000, pc, primes_small, spf, parity) == 233ULL);
    assert(F(1000000ULL, pc, primes_small, spf, parity) == 112168ULL);

    std::cout << F(Nmax, pc, primes_small, spf, parity) << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler611 {
    static long isqrt(long x) {
        long r = (long) Math.sqrt(x);
        while ((r + 1) * (r + 1) <= x)
            r++;
        while (r * r > x)
            r--;
        return r;
    }

    static int ilog2(long x) {
        if (x <= 0)
            return 0;
        return 63 - Long.numberOfLeadingZeros(x);
    }

    static List<Integer> sievePrimesInt(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int p = 2; (long) p * p <= n; p++) {
            if (isPrime[p]) {
                for (int j = p * p; j <= n; j += p) {
                    isPrime[j] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static long chiInt(long x) {
        if ((x & 1) == 0)
            return 0;
        return (x & 3) == 1 ? 1 : -1;
    }

    static long sumChi1To(long m) {
        long c1 = (m + 3) / 4;
        long c3 = (m + 1) / 4;
        return c1 - c3;
    }

    static class PrimeCountMod4 {
        long n;
        long sq;
        List<Long> vals;
        long[] gPi;
        long[] gChi;
        int[] id1;
        int[] id2;
        List<Integer> primes;
        long[] prefChiPrime;

        PrimeCountMod4(long n) {
            this.n = n;
            this.sq = isqrt(n);
            this.primes = sievePrimesInt((int) sq);

            this.vals = new ArrayList<>();
            for (long l = 1; l <= n;) {
                long w = n / l;
                vals.add(w);
                l = n / w + 1;
            }

            int m = vals.size();
            gPi = new long[m];
            gChi = new long[m];
            id1 = new int[(int) sq + 1];
            id2 = new int[(int) sq + 1];

            for (int i = 0; i < m; i++) {
                long w = vals.get(i);
                gPi[i] = (w >= 2) ? (w - 1) : 0;
                gChi[i] = (w >= 1) ? (sumChi1To(w) - 1) : 0;
                if (w <= sq)
                    id1[(int) w] = i;
                else
                    id2[(int) (n / w)] = i;
            }

            prefChiPrime = new long[primes.size() + 1];
            for (int i = 0; i < primes.size(); i++) {
                prefChiPrime[i + 1] = prefChiPrime[i] + chiInt(primes.get(i));
            }

            for (int i = 0; i < primes.size(); i++) {
                long p = primes.get(i);
                long p2 = p * p;
                if (p2 > n)
                    break;
                long chiP = chiInt(p);

                for (int j = 0; j < m; j++) {
                    long w = vals.get(j);
                    if (w < p2)
                        break;
                    long q = w / p;
                    int idx = (q <= sq) ? id1[(int) q] : id2[(int) (n / q)];
                    gPi[j] -= gPi[idx] - i;
                    if (chiP != 0) {
                        gChi[j] -= chiP * (gChi[idx] - prefChiPrime[i]);
                    }
                }
            }
        }

        int id(long x) {
            return (x <= sq) ? id1[(int) x] : id2[(int) (n / x)];
        }

        long pi(long x) {
            return gPi[id(x)];
        }

        long sumChiPrimes(long x) {
            return gChi[id(x)];
        }

        long pi1(long x) {
            if (x < 5)
                return 0;
            long pix = pi(x);
            long oddPrimes = pix - 1;
            long s = sumChiPrimes(x);
            return (oddPrimes + s) / 2;
        }
    }

    static int[] sieveSpf(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > n)
                    break;
                spf[(int) v] = p;
                if (p == spf[i])
                    break;
            }
        }
        return spf;
    }

    static long countB(long N, byte[] parityOdd1Mod4) {
        long lim = isqrt(N);
        long ans = 0;
        for (long m = 1; m <= lim; m += 2) {
            if (parityOdd1Mod4[(int) m] == 0)
                continue;
            long base = m * m;
            long t = N / base;
            ans += ilog2(t) + 1;
        }
        return ans;
    }

    static long countAEvenV2(long N, PrimeCountMod4 pc, List<Integer> primesSmall) {
        long S1 = 0;
        long tmax = isqrt(N);
        for (long t = 1; t <= tmax; t++) {
            long R = N / (t * t);
            if (R < 5)
                break;
            long L = N / ((t + 1) * (t + 1));
            long cnt = pc.pi1(R) - pc.pi1(L);
            S1 += t * cnt;
        }

        long corr = 0;
        long cbrtN = 10000;
        for (int p : primesSmall) {
            if (p > cbrtN)
                break;
            if (p % 4 != 1)
                continue;
            long T = isqrt(N / p);
            corr += T / p;
        }

        long total = S1 - corr;

        for (int p : primesSmall) {
            if (p % 4 != 1)
                continue;
            long powVal = (long) p * p * p * p * p;
            if (powVal > N || powVal < 0)
                break;
            long p4 = (long) p * p * p * p;
            while (powVal <= N && powVal > 0) {
                long T = isqrt(N / powVal);
                total += T - T / p;
                if (powVal > N / p4)
                    break;
                powVal *= p4;
            }
        }
        return total;
    }

    static long F(long N, PrimeCountMod4 pc, List<Integer> primesSmall, int[] spf, byte[] parity) {
        long aEven = countAEvenV2(N, pc, primesSmall);
        long aOdd = countAEvenV2(N / 2, pc, primesSmall);
        return aEven + aOdd + countB(N, parity);
    }

    public static String solve() {
        long Nmax = 1000000000000L;
        PrimeCountMod4 pc = new PrimeCountMod4(Nmax);

        int LIM = 1000000;
        int[] spf = sieveSpf(LIM);
        List<Integer> primesSmall = sievePrimesInt(LIM);

        byte[] parity = new byte[LIM + 1];
        for (int x = 2; x <= LIM; x++) {
            int n = x;
            byte par = 0;
            while (n > 1) {
                int p = spf[n];
                int e = 0;
                while (n % p == 0) {
                    n /= p;
                    e ^= 1;
                }
                if (e != 0 && (p % 4 == 1)) {
                    par ^= 1;
                }
            }
            parity[x] = par;
        }

        long ans = F(Nmax, pc, primesSmall, spf, parity);
        return Long.toString(ans);
    }

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