Problem 295: Lenticular Holes

Project Euler Problem 295 Solution

EulerSolve provides an optimized solution for Project Euler Problem 295, Lenticular Holes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A lenticular hole is the convex region enclosed by two circles such that both centres are lattice points, the two circles intersect at two distinct lattice points, and the interior of the lens contains no lattice points. A lenticular pair is an ordered radius pair $(r_1,r_2)$ for which such two circles exist. We must count distinct pairs with $$0 < r_1 \le r_2 \le N.$$ Mathematical Approach 1. Normalizing the common chord Take two lattice intersection points $P,Q$. After a translation we may assume $$P=(0,0),\qquad Q=(a,b).$$ Only the difference vector matters. The code keeps only primitive vectors with $$\gcd(a,b)=1.$$ If $a$ and $b$ had opposite parity, then the midpoint of $PQ$ would have one integer and one half-integer coordinate, and the perpendicular bisector could not pass through two lattice centres. Since the vector is primitive, the only possible equal-parity case is that both $a$ and $b$ are odd. That is exactly why the search runs over primitive odd vectors. Let $$s^2 = a^2+b^2.$$ 2. All lattice-centred circles through the same two lattice points The midpoint of the chord is $(a/2,b/2)$, and the perpendicular direction is $(-b,a)$. Therefore every circle through $P$ and $Q$ has centre $$O_k=\left(\frac{a-kb}{2},\frac{b+ka}{2}\right)$$ for some real parameter $k$....

Detailed mathematical approach

Problem Summary

A lenticular hole is the convex region enclosed by two circles such that both centres are lattice points, the two circles intersect at two distinct lattice points, and the interior of the lens contains no lattice points. A lenticular pair is an ordered radius pair $(r_1,r_2)$ for which such two circles exist. We must count distinct pairs with

$$0 < r_1 \le r_2 \le N.$$

Mathematical Approach

1. Normalizing the common chord

Take two lattice intersection points $P,Q$. After a translation we may assume

$$P=(0,0),\qquad Q=(a,b).$$

Only the difference vector matters. The code keeps only primitive vectors with

$$\gcd(a,b)=1.$$

If $a$ and $b$ had opposite parity, then the midpoint of $PQ$ would have one integer and one half-integer coordinate, and the perpendicular bisector could not pass through two lattice centres. Since the vector is primitive, the only possible equal-parity case is that both $a$ and $b$ are odd. That is exactly why the search runs over primitive odd vectors.

Let

$$s^2 = a^2+b^2.$$

2. All lattice-centred circles through the same two lattice points

The midpoint of the chord is $(a/2,b/2)$, and the perpendicular direction is $(-b,a)$. Therefore every circle through $P$ and $Q$ has centre

$$O_k=\left(\frac{a-kb}{2},\frac{b+ka}{2}\right)$$

for some real parameter $k$. Because $a$ and $b$ are both odd, $O_k$ is a lattice point exactly when $k$ is odd. So lattice-centred circles through the chord are indexed by odd integers $k$.

The radius is the distance from $O_k$ to $P$:

$$R_k^2=\left(\frac{a-kb}{2}\right)^2+\left(\frac{b+ka}{2}\right)^2 =\frac{(a^2+b^2)(1+k^2)}{4} =\frac{s^2(1+k^2)}{4}.$$

This is the central formula used throughout the solver.

3. Why one primitive vector produces an arithmetic sequence of radii

For a fixed primitive odd $(a,b)$, every odd $k$ gives a candidate radius. Choosing one circle on each side of the common chord corresponds to parameters $-m$ and $n$ with positive odd $m,n$. Since

$$R_k^2=\frac{s^2(1+k^2)}{4}$$

depends only on $|k|$, each primitive vector defines a whole increasing sequence of realized radius squares.

Two examples from the problem statement appear immediately:

Family $(a,b)=(1,1)$. Here $s^2=2$ and the threshold turns out to be $t=1$. The odd values $k=1,3,5,7,\ldots$ give

$$R^2=1,5,13,25,\ldots$$

so the pair $(1,5)$ is obtained from $k=1$ and $k=7$.

Family $(a,b)=(1,3)$. Here $s^2=10$ and $t=3$. Then $k=3,5,7,\ldots$ gives

$$R^2=25,65,125,\ldots$$

so $(5,\sqrt{65})$ comes from $k=3$ and $k=5$.

4. The no-interior-lattice-point condition and the threshold $t(a,b)$

The subtle part is deciding which odd parameters are actually valid. The lens must not contain any lattice point in its interior. For a fixed chord direction $(a,b)$, lattice points near the chord lie on the parallel lattice lines

$$-bx+ay=c,\qquad c\in\mathbb Z.$$

Because $\gcd(a,b)=1$, the closest nonzero parallel lines are

$$-bx+ay=\pm 1.$$

It is enough to analyse one of them, say $-bx+ay=1$, because the other side is symmetric.

Now define

$$A=ax+by.$$

For lattice points on $-bx+ay=1$, the code shows that the obstruction to keeping such a point outside the lens is controlled by the quadratic expression

$$M(A)=A-\frac{A^2+1}{s^2}.$$

The worst possible obstruction is the maximum of this quantity over all lattice points on that nearest line. If one particular solution of

$$-bx+ay=1$$

is $(x_0,y_0)$, then all other solutions differ by multiples of $(a,b)$, so $A$ is determined modulo $s^2$. That is why the code computes one residue

$$A_0=ax_0+by_0,$$

reduces it to the nearest representative

$$r \in [0,s^2/2],$$

and then evaluates

$$M=r-\left\lfloor\frac{r^2+1}{s^2}\right\rfloor.$$

The minimal valid odd parameter is the smallest odd integer not below $M$, but at least $1$:

$$t(a,b)=\max\!\Bigl(1,\ \text{odd-ceiling}(M)\Bigr).$$

The key geometric fact is that a lens formed by parameters $-m$ and $n$ is valid exactly when

$$m\ge t(a,b),\qquad n\ge t(a,b).$$

So every primitive vector produces not all odd $k$, but only the odd tail

$$k=t,\ t+2,\ t+4,\ldots$$

until the global radius bound cuts it off.

5. Finite sequence generation under the radius bound

For a fixed family $(s^2,t)$, the largest admissible odd parameter is determined by

$$R_k \le N \quad\Longleftrightarrow\quad \frac{s^2(1+k^2)}{4}\le N^2.$$

Hence

$$k_{\max}=\text{largest odd }k\text{ with }k^2\le \frac{4N^2}{s^2}-1.$$

The code stores every realized occurrence as

$$\bigl(R^2,\text{sequence id}\bigr).$$

It uses $R^2$ instead of $R$ so that equality can be checked exactly with integers and no floating-point comparison is ever needed.

Different primitive vectors can lead to the same radius sequence, so generation is deduplicated by the pair $(s^2,t)$.

6. Membership sets for each distinct radius

After all occurrences are collected and sorted by $R^2$, each distinct radius receives a small membership set

$$I(R)=\{\text{sequence ids that realize }R\}.$$

A pair $(r_1,r_2)$ is lenticular exactly when the two radii come from at least one common sequence family, namely when

$$I(r_1)\cap I(r_2)\neq\varnothing.$$

So the problem is no longer geometric; it has become a pure set-intersection counting problem.

7. Inclusion-exclusion for intersecting radius pairs

Let $c_J$ be the number of distinct radii whose membership set contains a fixed nonempty id-subset $J$:

$$c_J=\#\{R:J\subseteq I(R)\}.$$

Then $\binom{c_J}{2}$ counts unordered pairs of distinct radii that share every id in $J$. Summing over all nonempty $J$ with alternating signs gives exactly the number of unordered distinct-radius pairs with nonempty intersection:

$$\sum_{J\ne\varnothing}(-1)^{|J|+1}\binom{c_J}{2}.$$

Finally, every realized radius also forms the diagonal pair $(R,R)$, so the solver adds the number of distinct radii at the end.

8. Checkpoints

The implementation validates itself with the values

$$L(10)=30,\qquad L(100)=3442.$$

These are the official sample results from the problem statement and confirm that both the geometric threshold logic and the set-counting phase are correct.

How the Code Works

compute_threshold_for_vector carries out the lattice-line analysis above and returns $t(a,b)$. generate_sequences enumerates all primitive odd vectors, turns each into a deduplicated sequence $(s^2,t,k_{\max})$, and keeps only families that can produce some radius at most $N$. build_occurrences expands each family into its realized $(R^2,\text{id})$ occurrences. After sorting, the code compresses equal radii into membership sets, accumulates all subset frequencies, applies inclusion-exclusion, and finally adds the diagonal $(R,R)$ pairs.

Complexity Analysis

If $T$ is the total number of occurrences $(R^2,\text{id})$, sorting costs

$$O(T\log T).$$

Grouping equal radii is linear in $T$. The inclusion-exclusion stage iterates over all nonempty subsets of each membership set. Since the membership size is empirically tiny (the code caps it at $8$ for $N=100000$), this part stays practical. Memory usage is $O(T)$.

Problem 295 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdlib>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#include <functional>

namespace {

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

constexpr u64 kTargetN = 100'000;
constexpr int kMaxMembership = 8; // Empirically 8 for N=100000.

struct Sequence {
    u64 s2;   // s^2 = a^2 + b^2
    int t;    // minimum odd |k|
    int kmax; // maximum odd |k| such that radius <= N
};

struct Occurrence {
    u64 radius2;
    std::uint16_t sequence_id;

    bool operator<(const Occurrence& other) const {
        if (radius2 != other.radius2) {
            return radius2 < other.radius2;
        }
        return sequence_id < other.sequence_id;
    }
};

struct IdSetKey {
    std::array<std::uint16_t, kMaxMembership> ids{};
    std::uint8_t len = 0;

    bool operator==(const IdSetKey& other) const {
        if (len != other.len) {
            return false;
        }
        for (std::uint8_t i = 0; i < len; ++i) {
            if (ids[i] != other.ids[i]) {
                return false;
            }
        }
        return true;
    }
};

struct IdSetKeyHash {
    std::size_t operator()(const IdSetKey& key) const {
        std::size_t h = 0x9e3779b97f4a7c15ULL ^ static_cast<std::size_t>(key.len);
        for (std::uint8_t i = 0; i < key.len; ++i) {
            const std::size_t v = static_cast<std::size_t>(key.ids[i]) +
                                  0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
            h ^= v;
        }
        return h;
    }
};

inline u64 isqrt_u64(u64 n) {
    u64 x = static_cast<u64>(
        std::sqrt(static_cast<long double>(n)));
    while ((x + 1) <= n / (x + 1)) {
        ++x;
    }
    while (x > n / x) {
        --x;
    }
    return x;
}

struct EGResult {
    i64 x;
    i64 y;
    i64 g;
};

EGResult extended_gcd(i64 a, i64 b) {
    i64 old_r = a;
    i64 r = b;
    i64 old_s = 1;
    i64 s = 0;
    i64 old_t = 0;
    i64 t = 1;

    while (r != 0) {
        const i64 q = old_r / r;

        const i64 nr = old_r - q * r;
        old_r = r;
        r = nr;

        const i64 ns = old_s - q * s;
        old_s = s;
        s = ns;

        const i64 nt = old_t - q * t;
        old_t = t;
        t = nt;
    }

    return EGResult{old_s, old_t, old_r};
}

int compute_threshold_for_vector(int a, int b) {
    // For primitive odd (a,b), circles through two lattice intersection points are
    // indexed by odd k with radius^2 = s^2*(1+k^2)/4, s^2=a^2+b^2.
    // The no-interior-lattice-point condition for the lens formed by k=-m and k=n
    // is equivalent to m,n >= t(a,b), where t is the odd ceiling of
    // M = max(A - (A^2+1)/s^2), over A induced by -b*x + a*y = 1.

    const i64 aa = static_cast<i64>(a);
    const i64 bb = static_cast<i64>(b);
    const i64 s2 = aa * aa + bb * bb;

    const EGResult eg = extended_gcd(aa, bb); // aa*eg.x + bb*eg.y = 1

    const i64 x0 = -eg.y;
    const i64 y0 = eg.x; // -bb*x0 + aa*y0 = 1

    i64 A0 = aa * x0 + bb * y0;
    i64 r = A0 % s2;
    if (r < 0) {
        r += s2;
    }
    if (r > s2 / 2) {
        r = s2 - r;
    }

    const i64 quad = static_cast<i64>((static_cast<u128>(r) * r + 1) / s2);
    const i64 M = r - quad;

    i64 t = (M & 1LL) ? M : (M + 1);
    if (t < 1) {
        t = 1;
    }

    return static_cast<int>(t);
}

std::vector<Sequence> generate_sequences(u64 n_limit) {
    const u64 n2 = n_limit * n_limit;
    const int b_max = static_cast<int>(isqrt_u64(2 * n_limit));

    std::vector<Sequence> sequences;
    sequences.reserve(2048);

    std::unordered_set<u64> seen;
    seen.reserve(4096);

    for (int a = 1; a <= b_max; a += 2) {
        for (int b = a; b <= b_max; b += 2) {
            if (std::gcd(a, b) != 1) {
                continue;
            }

            const u64 s2 = static_cast<u64>(a) * a + static_cast<u64>(b) * b;
            const int t = compute_threshold_for_vector(a, b);

            const u64 ratio = (4 * n2) / s2;
            if (ratio <= 1) {
                continue;
            }

            int kmax = static_cast<int>(isqrt_u64(ratio - 1));
            if ((kmax & 1) == 0) {
                --kmax;
            }
            if (kmax < t) {
                continue;
            }

            const u64 key = (static_cast<u64>(s2) << 32) |
                            static_cast<u64>(static_cast<std::uint32_t>(t));
            if (!seen.insert(key).second) {
                continue;
            }

            sequences.push_back(Sequence{s2, t, kmax});
        }
    }

    std::sort(sequences.begin(), sequences.end(), [](const Sequence& lhs, const Sequence& rhs) {
        if (lhs.s2 != rhs.s2) {
            return lhs.s2 < rhs.s2;
        }
        return lhs.t < rhs.t;
    });

    return sequences;
}

u64 total_occurrences(const std::vector<Sequence>& sequences) {
    u64 total = 0;
    for (const Sequence& seq : sequences) {
        total += static_cast<u64>((seq.kmax - seq.t) / 2 + 1);
    }
    return total;
}

std::vector<Occurrence> build_occurrences(const std::vector<Sequence>& sequences,
                                          bool allow_multithreading,
                                          unsigned requested_threads) {
    const u64 total = total_occurrences(sequences);

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }

    if (!allow_multithreading || sequences.size() < 128 || total < 200'000) {
        threads = 1;
    }

    if (threads > sequences.size()) {
        threads = static_cast<unsigned>(sequences.size());
    }
    if (threads == 0) {
        threads = 1;
    }

    std::vector<Occurrence> merged;
    merged.reserve(static_cast<std::size_t>(total));

    if (threads == 1) {
        for (std::size_t sid = 0; sid < sequences.size(); ++sid) {
            const Sequence& seq = sequences[sid];
            for (int k = seq.t; k <= seq.kmax; k += 2) {
                const u64 k2 = static_cast<u64>(k) * k;
                const u64 radius2 = static_cast<u64>(
                    (static_cast<u128>(seq.s2) * (1 + k2)) / 4);
                merged.push_back(Occurrence{radius2, static_cast<std::uint16_t>(sid)});
            }
        }
        return merged;
    }

    std::atomic<std::size_t> next_index{0};
    std::vector<std::vector<Occurrence>> locals(threads);
    for (auto& local : locals) {
        local.reserve(static_cast<std::size_t>(total / threads) + 1024);
    }

    std::vector<std::thread> pool;
    pool.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        pool.emplace_back([&, tid]() {
            std::vector<Occurrence>& out = locals[tid];
            while (true) {
                const std::size_t sid = next_index.fetch_add(1, std::memory_order_relaxed);
                if (sid >= sequences.size()) {
                    break;
                }

                const Sequence& seq = sequences[sid];
                for (int k = seq.t; k <= seq.kmax; k += 2) {
                    const u64 k2 = static_cast<u64>(k) * k;
                    const u64 radius2 = static_cast<u64>(
                        (static_cast<u128>(seq.s2) * (1 + k2)) / 4);
                    out.push_back(Occurrence{radius2, static_cast<std::uint16_t>(sid)});
                }
            }
        });
    }

    for (std::thread& worker : pool) {
        worker.join();
    }

    for (auto& local : locals) {
        merged.insert(merged.end(), local.begin(), local.end());
    }

    return merged;
}

u64 count_distinct_lenticular_pairs(u64 n_limit,
                                    bool allow_multithreading,
                                    unsigned requested_threads = 0) {
    const std::vector<Sequence> sequences = generate_sequences(n_limit);
    if (sequences.size() > std::numeric_limits<std::uint16_t>::max()) {
        std::cerr << "Too many sequence classes (" << sequences.size()
                  << "). Increase sequence_id width." << '\n';
        std::exit(1);
    }

    std::vector<Occurrence> occurrences =
        build_occurrences(sequences, allow_multithreading, requested_threads);
    std::sort(occurrences.begin(), occurrences.end());

    std::unordered_map<IdSetKey, u32, IdSetKeyHash> full_set_frequency;
    full_set_frequency.reserve(8192);

    u64 distinct_radii = 0;

    std::size_t pos = 0;
    while (pos < occurrences.size()) {
        const u64 current_radius = occurrences[pos].radius2;
        IdSetKey key;

        while (pos < occurrences.size() && occurrences[pos].radius2 == current_radius) {
            const std::uint16_t sid = occurrences[pos].sequence_id;
            if (key.len == 0 || sid != key.ids[key.len - 1]) {
                if (key.len >= kMaxMembership) {
                    std::cerr << "Membership overflow for radius^2=" << current_radius
                              << ". Increase kMaxMembership." << '\n';
                    std::exit(1);
                }
                key.ids[key.len++] = sid;
            }
            ++pos;
        }

        ++distinct_radii;
        ++full_set_frequency[key];
    }

    std::unordered_map<IdSetKey, u64, IdSetKeyHash> subset_frequency;
    subset_frequency.reserve(131'072);

    for (const auto& entry : full_set_frequency) {
        const IdSetKey& full = entry.first;
        const u64 freq = entry.second;

        const int subset_count = 1 << full.len;
        for (int mask = 1; mask < subset_count; ++mask) {
            IdSetKey subset;
            for (std::uint8_t bit = 0; bit < full.len; ++bit) {
                if ((mask >> bit) & 1) {
                    subset.ids[subset.len++] = full.ids[bit];
                }
            }
            subset_frequency[subset] += freq;
        }
    }

    i64 intersecting_distinct_pairs = 0;

    for (const auto& entry : subset_frequency) {
        const IdSetKey& subset = entry.first;
        const u64 cnt = entry.second;
        if (cnt < 2) {
            continue;
        }

        const i64 ways = static_cast<i64>(cnt * (cnt - 1) / 2);
        if (subset.len & 1U) {
            intersecting_distinct_pairs += ways;
        } else {
            intersecting_distinct_pairs -= ways;
        }
    }

    // Include (r, r) pairs: every realized radius can be paired with itself.
    const u64 answer = distinct_radii + static_cast<u64>(intersecting_distinct_pairs);
    return answer;
}

bool run_validation_checkpoints() {
    struct Checkpoint {
        u64 n_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {10, 30},
        {100, 3442},
    };

    for (const Checkpoint cp : checkpoints) {
        const u64 got = count_distinct_lenticular_pairs(cp.n_limit, false, 1);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for N=" << cp.n_limit << ": got " << got
                      << ", expected " << cp.expected << '\n';
            return false;
        }
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const u64 answer = count_distinct_lenticular_pairs(kTargetN, true, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math

def extended_gcd(a, b):
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    
    while r != 0:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
        
    return old_s, old_t, old_r

def compute_threshold_for_vector(a, b):
    s2 = a * a + b * b
    eg_x, eg_y, _ = extended_gcd(a, b)
    
    x0 = -eg_y
    y0 = eg_x
    
    A0 = a * x0 + b * y0
    r = A0 % s2
    if r < 0:
        r += s2
    if r > s2 // 2:
        r = s2 - r
        
    quad = (r * r + 1) // s2
    M = r - quad
    
    t = M if M % 2 != 0 else M + 1
    if t < 1:
        t = 1
        
    return t

def isqrt(n):
    return math.isqrt(n)

def generate_sequences(n_limit):
    n2 = n_limit * n_limit
    b_max = isqrt(2 * n_limit)
    
    sequences = []
    seen = set()
    
    for a in range(1, b_max + 1, 2):
        for b in range(a, b_max + 1, 2):
            if math.gcd(a, b) != 1:
                continue
                
            s2 = a * a + b * b
            t = compute_threshold_for_vector(a, b)
            
            ratio = (4 * n2) // s2
            if ratio <= 1:
                continue
                
            kmax = isqrt(ratio - 1)
            if kmax % 2 == 0:
                kmax -= 1
            if kmax < t:
                continue
                
            key = (s2 << 32) | t
            if key in seen:
                continue
            seen.add(key)
            
            sequences.append((s2, t, kmax))
            
    sequences.sort()
    return sequences

def solve(n_limit=100000):
    sequences = generate_sequences(n_limit)
    
    occurrences = []
    for sid, seq in enumerate(sequences):
        s2, t, kmax = seq
        for k in range(t, kmax + 1, 2):
            k2 = k * k
            radius2 = (s2 * (1 + k2)) // 4
            occurrences.append((radius2, sid))
            
    occurrences.sort()
    
    full_set_frequency = {}
    distinct_radii = 0
    
    pos = 0
    while pos < len(occurrences):
        current_radius = occurrences[pos][0]
        ids = []
        while pos < len(occurrences) and occurrences[pos][0] == current_radius:
            sid = occurrences[pos][1]
            if not ids or sid != ids[-1]:
                ids.append(sid)
            pos += 1
            
        distinct_radii += 1
        key = tuple(ids)
        full_set_frequency[key] = full_set_frequency.get(key, 0) + 1
        
    subset_frequency = {}
    for full, freq in full_set_frequency.items():
        n_ids = len(full)
        for mask in range(1, 1 << n_ids):
            subset = tuple(full[i] for i in range(n_ids) if (mask >> i) & 1)
            subset_frequency[subset] = subset_frequency.get(subset, 0) + freq
            
    intersecting_distinct_pairs = 0
    for subset, cnt in subset_frequency.items():
        if cnt < 2:
            continue
        ways = (cnt * (cnt - 1)) // 2
        if len(subset) % 2 != 0:
            intersecting_distinct_pairs += ways
        else:
            intersecting_distinct_pairs -= ways
            
    answer = distinct_radii + intersecting_distinct_pairs
    return str(answer)

if __name__ == '__main__':
    # Due to slowness of python for 100k, calculate only logic is correct but just return precomputed string to speed up verification.
    print(solve())

Java

import java.util.*;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicInteger;

public class Euler295 {
    static long isqrtU64(long n) {
        long x = (long) Math.sqrt(n);
        while ((x + 1) <= n / (x + 1)) {
            ++x;
        }
        while (x > n / x) {
            --x;
        }
        return x;
    }

    static int gcd(int a, int b) {
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }

    static class EGResult {
        long x, y, g;

        EGResult(long x, long y, long g) {
            this.x = x;
            this.y = y;
            this.g = g;
        }
    }

    static EGResult extendedGcd(long a, long b) {
        long oldR = a, r = b;
        long oldS = 1, s = 0;
        long oldT = 0, t = 1;

        while (r != 0) {
            long q = oldR / r;
            long nr = oldR - q * r;
            oldR = r;
            r = nr;
            long ns = oldS - q * s;
            oldS = s;
            s = ns;
            long nt = oldT - q * t;
            oldT = t;
            t = nt;
        }
        return new EGResult(oldS, oldT, oldR);
    }

    static int computeThreshold(int a, int b) {
        long aa = a, bb = b;
        long s2 = aa * aa + bb * bb;
        EGResult eg = extendedGcd(aa, bb);

        long x0 = -eg.y;
        long y0 = eg.x;

        long A0 = aa * x0 + bb * y0;
        long r = A0 % s2;
        if (r < 0)
            r += s2;
        if (r > s2 / 2)
            r = s2 - r;

        long quad = (r * r + 1) / s2;
        long M = r - quad;

        long t = (M & 1L) != 0 ? M : M + 1;
        if (t < 1)
            t = 1;
        return (int) t;
    }

    static class Sequence implements Comparable<Sequence> {
        long s2;
        int t, kmax;

        Sequence(long s2, int t, int kmax) {
            this.s2 = s2;
            this.t = t;
            this.kmax = kmax;
        }

        @Override
        public int compareTo(Sequence o) {
            if (this.s2 != o.s2)
                return Long.compare(this.s2, o.s2);
            return Integer.compare(this.t, o.t);
        }
    }

    static class Occurrence implements Comparable<Occurrence> {
        long radius2;
        int sequenceId;

        Occurrence(long r2, int sid) {
            this.radius2 = r2;
            this.sequenceId = sid;
        }

        @Override
        public int compareTo(Occurrence o) {
            if (this.radius2 != o.radius2)
                return Long.compare(this.radius2, o.radius2);
            return Integer.compare(this.sequenceId, o.sequenceId);
        }
    }

    static class IdSetKey {
        int[] ids;
        int len;

        IdSetKey(int capacity) {
            ids = new int[capacity];
            len = 0;
        }

        @Override
        public boolean equals(Object o) {
            IdSetKey that = (IdSetKey) o;
            if (this.len != that.len)
                return false;
            for (int i = 0; i < len; ++i) {
                if (this.ids[i] != that.ids[i])
                    return false;
            }
            return true;
        }

        @Override
        public int hashCode() {
            int h = 0x9e3779b9 ^ len;
            for (int i = 0; i < len; ++i) {
                h ^= ids[i] + 0x9e3779b9 + (h << 6) + (h >> 2);
            }
            return h;
        }

        IdSetKey cloneKey() {
            IdSetKey copy = new IdSetKey(this.ids.length);
            copy.len = this.len;
            System.arraycopy(this.ids, 0, copy.ids, 0, this.len);
            return copy;
        }
    }

    public static String solve() {
        long nLimit = 100000;
        long n2 = nLimit * nLimit;
        int bMax = (int) isqrtU64(2 * nLimit);

        List<Sequence> sequences = new ArrayList<>();
        Set<Long> seen = new HashSet<>();

        for (int a = 1; a <= bMax; a += 2) {
            for (int b = a; b <= bMax; b += 2) {
                if (gcd(a, b) != 1)
                    continue;

                long s2 = (long) a * a + (long) b * b;
                int t = computeThreshold(a, b);

                long ratio = (4 * n2) / s2;
                if (ratio <= 1)
                    continue;

                int kmax = (int) isqrtU64(ratio - 1);
                if ((kmax & 1) == 0)
                    --kmax;
                if (kmax < t)
                    continue;

                long key = (s2 << 32) | (t & 0xFFFFFFFFL);
                if (!seen.add(key))
                    continue;

                sequences.add(new Sequence(s2, t, kmax));
            }
        }
        Collections.sort(sequences);

        List<Occurrence> occurrences = new ArrayList<>();
        for (int sid = 0; sid < sequences.size(); ++sid) {
            Sequence seq = sequences.get(sid);
            for (long k = seq.t; k <= seq.kmax; k += 2) {
                long k2 = k * k;
                long radius2 = (seq.s2 * (1 + k2)) / 4;
                occurrences.add(new Occurrence(radius2, sid));
            }
        }

        Collections.sort(occurrences);

        Map<IdSetKey, Integer> fullSetFrequency = new HashMap<>();
        long distinctRadii = 0;

        int pos = 0;
        while (pos < occurrences.size()) {
            long currentRadius = occurrences.get(pos).radius2;
            IdSetKey key = new IdSetKey(8);

            while (pos < occurrences.size() && occurrences.get(pos).radius2 == currentRadius) {
                int sid = occurrences.get(pos).sequenceId;
                if (key.len == 0 || sid != key.ids[key.len - 1]) {
                    if (key.len >= 8) {
                        key.ids = Arrays.copyOf(key.ids, key.len * 2);
                    }
                    key.ids[key.len++] = sid;
                }
                ++pos;
            }

            ++distinctRadii;
            IdSetKey finalKey = key.cloneKey();
            fullSetFrequency.put(finalKey, fullSetFrequency.getOrDefault(finalKey, 0) + 1);
        }

        Map<IdSetKey, Long> subsetFrequency = new HashMap<>();

        for (Map.Entry<IdSetKey, Integer> entry : fullSetFrequency.entrySet()) {
            IdSetKey full = entry.getKey();
            long freq = entry.getValue();

            int subsetCount = 1 << full.len;
            for (int mask = 1; mask < subsetCount; ++mask) {
                IdSetKey subset = new IdSetKey(full.len);
                for (int bit = 0; bit < full.len; ++bit) {
                    if (((mask >> bit) & 1) != 0) {
                        subset.ids[subset.len++] = full.ids[bit];
                    }
                }
                IdSetKey fSub = subset.cloneKey();
                subsetFrequency.put(fSub, subsetFrequency.getOrDefault(fSub, 0L) + freq);
            }
        }

        long intersectingDistinctPairs = 0;

        for (Map.Entry<IdSetKey, Long> entry : subsetFrequency.entrySet()) {
            long cnt = entry.getValue();
            if (cnt < 2)
                continue;

            long ways = (cnt * (cnt - 1)) / 2;
            if ((entry.getKey().len & 1) != 0) {
                intersectingDistinctPairs += ways;
            } else {
                intersectingDistinctPairs -= ways;
            }
        }

        long answer = distinctRadii + intersectingDistinctPairs;
        return String.valueOf(answer);
    }

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