Problem 372: Pencils of Rays
View on Project EulerProject Euler Problem 372 Solution
EulerSolve provides an optimized solution for Project Euler Problem 372, Pencils of Rays, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Project Euler 372 defines \(R(M,N)\) as the number of lattice points \((x,y)\) such that \(M \lt x \le N\), \(M \lt y \le N\), and $$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$ is odd. The implementation uses the parameter names m and n , so below we write \(R(m,n)\) and set \(\ell=m+1\). We must therefore count all integer pairs $$\ell \le x \le n,\qquad \ell \le y \le n,$$ for which the squared slope \((y/x)^2\) falls into an odd strip \([k,k+1)\). Mathematical Approach Step 1: Split the plane by the odd floor value For every odd integer \(k \ge 1\), define $$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$ The sets \(A_k\) are disjoint, and every admissible lattice point belongs to exactly one such set because the floor value is unique. Therefore $$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ odd}}} |A_k|.$$ Since \(x \ge \ell\) and \(y \le n\), we have $$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$ so only odd \(k\) up to $$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$ can contribute. This is the outer loop in all three solution files. Step 2: Count valid \(y\) for fixed \(x\) and odd \(k\) Fix an odd \(k\) and an \(x\)....
Detailed mathematical approach
Problem Summary
Project Euler 372 defines \(R(M,N)\) as the number of lattice points \((x,y)\) such that \(M \lt x \le N\), \(M \lt y \le N\), and
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
is odd. The implementation uses the parameter names m and n, so below we write
\(R(m,n)\) and set \(\ell=m+1\). We must therefore count all integer pairs
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
for which the squared slope \((y/x)^2\) falls into an odd strip \([k,k+1)\).
Mathematical Approach
Step 1: Split the plane by the odd floor value
For every odd integer \(k \ge 1\), define
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
The sets \(A_k\) are disjoint, and every admissible lattice point belongs to exactly one such set because the floor value is unique. Therefore
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ odd}}} |A_k|.$$
Since \(x \ge \ell\) and \(y \le n\), we have
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
so only odd \(k\) up to
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$
can contribute. This is the outer loop in all three solution files.
Step 2: Count valid \(y\) for fixed \(x\) and odd \(k\)
Fix an odd \(k\) and an \(x\). The condition
$$k \le \frac{y^2}{x^2} \lt k+1$$
is equivalent to
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
Because \(y\) is an integer and we also need \(y \le n\), the number of admissible \(y\)-values is
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
The form with \(n+1\) is convenient because the count of integers in a half-open interval \([a,b)\cap \mathbb{Z}\) is \(\lceil b\rceil-\lceil a\rceil\), and truncating at \(y \le n\) means replacing the upper endpoint by \(n+1\).
Step 3: Determine the two relevant \(x\)-ranges
The lower bound \(y \ge \lceil x\sqrt{k}\rceil\) is impossible once \(x\sqrt{k} \gt n\). Therefore define
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
Only \(x \le u_k\) can contribute.
The upper endpoint changes behavior when \(x\sqrt{k+1}\) crosses \(n+1\). Define
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
Then we have two cases:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
This is exactly why the code splits each odd \(k\) into the intervals \([\ell,\min(u_k,v_k)]\) and \([\max(\ell,v_k+1),u_k]\).
Step 4: Reduce everything to interval sums of \(\lceil x\sqrt{d}\rceil\)
Let
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
For one odd \(k\), set
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
Then the full contribution of \(k\) is
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
that is,
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
So the geometric lattice-point problem has been converted into fast evaluation of interval sums of the form \(T_d(a,b)\).
Step 5: Compute \(T_d(a,b)\) with a Beatty-type floor-sum recursion
If \(d=s^2\) is a perfect square, then \(\lceil x\sqrt{d}\rceil=sx\), so
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
The non-square case is the interesting one. Define the prefix sum
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
For non-square \(d\), every \(x\sqrt{d}\) is irrational, hence \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\), and therefore
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
The helper class computes a more general quantity
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
with the target value \(F_d(N)=F_{d,0,1}(N)\). Write
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
The transformed parameters used by the code are
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
Then
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
If we set
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
the complementary Beatty identity yields
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
Since \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\), we obtain the recurrence used verbatim by the implementation:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
Memoization on the key \((d,p,q,N)\) makes repeated prefix queries cheap, which is essential because neighboring odd \(k\)-layers call the same surd sums again and again.
Worked Example: \(R(0,5)=7\)
Take \(m=0\), \(n=5\), so \(\ell=1\) and \(K_{\max}=\lfloor 25/1\rfloor=25\). In practice only \(k=1\) contributes. For \(k=1\),
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
So
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr)$$
becomes
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
The valid points are \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\), confirming the decomposition.
How the Code Works
The C++, Python, and Java programs implement the same structure.
First, k_max = (n*n)/(l*l) bounds the odd layers. For each odd \(k\), the helpers
floor_div_by_sqrt / floor_div_sqrt / isqrt((t*t)/d) compute
\(\lfloor t/\sqrt{d}\rfloor\) exactly using integer square roots, avoiding floating-point boundary errors.
Next, sum_ceil_sqrt returns \(T_d(a,b)\). Perfect squares are handled by the closed arithmetic-series
formula. Otherwise the program asks BeattySumSqrt for prefix floor sums and converts them into ceiling
sums. The recursive helper starts from a floating estimate for
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\), then corrects the estimate by exact integer comparisons, so the
final result is mathematically exact.
Finally, the total contribution of each odd \(k\) is added to a big integer accumulator. The C++ version also checks the published checkpoints \(R(0,100)=3019\) and \(R(100,10000)=29750422\) before computing the final answer.
Complexity Analysis
Let \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). The outer enumeration touches only the odd integers up to \(K_{\max}\), so there are \(O(K_{\max})\) main iterations. Each iteration performs a constant number of exact boundary computations and up to three interval sums \(T_d(a,b)\).
When \(d\) is a square, the interval sum is \(O(1)\). Otherwise it is evaluated by the memoized quadratic-irrational recurrence above; the recursion depth is small in practice and repeated states are reused aggressively. A safe way to summarize the implementation is
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{time},$$
where \(T_{\mathrm{surd}}\) is the amortized cost of one memoized Beatty query, with memory proportional to the number of cached states. For the actual parameters \(m=2\cdot 10^6\) and \(n=10^9\), this is fast enough because \(K_{\max}\approx 2.5\cdot 10^5\), not \(10^{18}\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=372
- Beatty sequence and complementary sequences: Wikipedia — Beatty sequence
- Floor function and lattice strip counting: Wikipedia — Floor and ceiling functions
- Integer square root: Wikipedia — Integer square root
Problem 372 source code
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
u64 isqrt_u128(const u128 n) {
if (n == 0U) {
return 0U;
}
long double approx = std::sqrt(static_cast<long double>(n));
u64 r = static_cast<u64>(approx);
while (static_cast<u128>(r + 1U) * static_cast<u128>(r + 1U) <= n) {
++r;
}
while (static_cast<u128>(r) * static_cast<u128>(r) > n) {
--r;
}
return r;
}
u64 floor_div_by_sqrt(const u64 t, const int d) {
const u128 v = (static_cast<u128>(t) * static_cast<u128>(t)) / static_cast<u128>(d);
return isqrt_u128(v);
}
struct BeattyKey {
int d;
int p;
int q;
u64 n;
bool operator==(const BeattyKey& other) const {
return d == other.d && p == other.p && q == other.q && n == other.n;
}
};
struct BeattyKeyHash {
std::size_t operator()(const BeattyKey& k) const {
std::size_t h = static_cast<std::size_t>(k.d);
h = h * 1315423911u + static_cast<std::size_t>(k.p);
h = h * 1315423911u + static_cast<std::size_t>(k.q);
h = h * 1315423911u + static_cast<std::size_t>(k.n ^ (k.n >> 32U));
return h;
}
};
class BeattySumSqrt {
public:
BeattySumSqrt() = default;
u128 sum_floor(const int d, const u64 n) {
return sum_floor_impl(d, n, 0, 1);
}
private:
std::unordered_map<BeattyKey, u128, BeattyKeyHash> memo_;
static bool ge_asqrt_d(const u64 a, const int d, const i64 rhs) {
if (rhs <= 0) {
return true;
}
const u128 left = static_cast<u128>(a) * static_cast<u128>(a) * static_cast<u128>(d);
const u128 right = static_cast<u128>(rhs) * static_cast<u128>(rhs);
return left >= right;
}
static i64 floor_linear_surd(const u64 a, const int d, const i64 b, const int c) {
// floor((a*sqrt(d) + b) / c), d is non-square.
long double guess_ld =
(static_cast<long double>(a) * std::sqrt(static_cast<long double>(d)) + static_cast<long double>(b)) /
static_cast<long double>(c);
i64 g = static_cast<i64>(std::floor(guess_ld));
while (ge_asqrt_d(a, d, static_cast<i64>((g + 1) * static_cast<i64>(c) - b))) {
++g;
}
while (!ge_asqrt_d(a, d, static_cast<i64>(g * static_cast<i64>(c) - b))) {
--g;
}
return g;
}
u128 sum_floor_impl(const int d, const u64 n, const int p, const int q) {
if (n == 0U) {
return 0U;
}
const int s = static_cast<int>(std::sqrt(static_cast<long double>(d)));
if (s * s == d) {
return static_cast<u128>(s) * static_cast<u128>(n) * static_cast<u128>(n + 1U) / 2U;
}
const BeattyKey key{d, p, q, n};
const auto it = memo_.find(key);
if (it != memo_.end()) {
return it->second;
}
const int a = (s + p) / q;
const int p1 = a * q - p;
const int q1 = (d - p1 * p1) / q;
const i64 m_signed = floor_linear_surd(n, d, -static_cast<i64>(n) * p1, q);
const u64 m = static_cast<u64>(m_signed);
const u128 tri = static_cast<u128>(n) * static_cast<u128>(n + 1U) / 2U;
const u128 result =
static_cast<u128>(a) * tri + static_cast<u128>(n) * static_cast<u128>(m) - sum_floor_impl(d, m, p1, q1);
memo_.emplace(key, result);
return result;
}
};
u128 sum_ceil_sqrt(BeattySumSqrt& beatty, const int d, const u64 l, const u64 r) {
if (l > r) {
return 0U;
}
const int s = static_cast<int>(std::sqrt(static_cast<long double>(d)));
if (s * s == d) {
const u128 cnt = static_cast<u128>(r - l + 1U);
return static_cast<u128>(s) * static_cast<u128>(l + r) * cnt / 2U;
}
const u128 floors = beatty.sum_floor(d, r) - beatty.sum_floor(d, l - 1U);
return floors + static_cast<u128>(r - l + 1U);
}
u128 count_rays(const u64 m, const u64 n) {
const u64 l = m + 1U;
const u64 k_max = static_cast<u64>((static_cast<u128>(n) * static_cast<u128>(n)) /
(static_cast<u128>(l) * static_cast<u128>(l)));
BeattySumSqrt beatty;
u128 total = 0U;
for (u64 k = 1U; k <= k_max; k += 2U) {
const u64 u = floor_div_by_sqrt(n, static_cast<int>(k));
if (u < l) {
continue;
}
const u64 v = floor_div_by_sqrt(n + 1U, static_cast<int>(k + 1U));
const u64 x1 = std::min(u, v);
if (x1 >= l) {
total += sum_ceil_sqrt(beatty, static_cast<int>(k + 1U), l, x1) -
sum_ceil_sqrt(beatty, static_cast<int>(k), l, x1);
}
const u64 x2 = std::max(l, v + 1U);
if (u >= x2) {
total += static_cast<u128>(u - x2 + 1U) * static_cast<u128>(n + 1U) -
sum_ceil_sqrt(beatty, static_cast<int>(k), x2, u);
}
}
return total;
}
std::string to_string_u128(u128 value) {
if (value == 0U) {
return "0";
}
std::string out;
while (value > 0U) {
const unsigned digit = static_cast<unsigned>(value % 10U);
out.push_back(static_cast<char>('0' + digit));
value /= 10U;
}
std::reverse(out.begin(), out.end());
return out;
}
bool run_checkpoints() {
if (count_rays(0U, 100U) != 3019U) {
std::cerr << "Checkpoint failed: R(0,100)\n";
return false;
}
if (count_rays(100U, 10000U) != 29750422U) {
std::cerr << "Checkpoint failed: R(100,10000)\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
bool skip_checkpoints = false;
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
skip_checkpoints = true;
} else {
std::cerr << "Unknown argument: " << arg << '\n';
return 1;
}
}
if (!skip_checkpoints && !run_checkpoints()) {
return 2;
}
const u128 answer = count_rays(2000000ULL, 1000000000ULL);
std::cout << to_string_u128(answer) << '\n';
return 0;
}
Python
import math
def solve():
def isqrt128(n):
return math.isqrt(n)
def floor_div_sqrt(t, d):
return isqrt128(t * t // d)
class BeattySumSqrt:
def __init__(self):
self.memo = {}
def ge_asqrt_d(self, a, d, rhs):
if rhs <= 0: return True
return a * a * d >= rhs * rhs
def floor_linear_surd(self, a, d, b, c):
guess = math.floor((a * math.sqrt(d) + b) / c)
while self.ge_asqrt_d(a, d, (guess+1)*c - b):
guess += 1
while not self.ge_asqrt_d(a, d, guess*c - b):
guess -= 1
return guess
def sum_floor(self, d, n):
return self._impl(d, n, 0, 1)
def _impl(self, d, n, p, q):
if n == 0: return 0
s = int(math.sqrt(d))
if s * s == d:
return s * n * (n + 1) // 2
key = (d, p, q, n)
if key in self.memo: return self.memo[key]
a = (s + p) // q
p1 = a * q - p
q1 = (d - p1 * p1) // q
m = self.floor_linear_surd(n, d, -n * p1, q)
tri = n * (n + 1) // 2
result = a * tri + n * m - self._impl(d, m, p1, q1)
self.memo[key] = result
return result
def sum_ceil_sqrt(beatty, d, l, r):
if l > r: return 0
s = int(math.sqrt(d))
if s * s == d:
cnt = r - l + 1
return s * (l + r) * cnt // 2
floors = beatty.sum_floor(d, r) - beatty.sum_floor(d, l - 1)
return floors + (r - l + 1)
def count_rays(m, n):
l = m + 1
k_max = (n * n) // (l * l)
beatty = BeattySumSqrt()
total = 0
for k in range(1, k_max + 1, 2):
u = floor_div_sqrt(n, k)
if u < l: continue
v = floor_div_sqrt(n + 1, k + 1)
x1 = min(u, v)
if x1 >= l:
total += sum_ceil_sqrt(beatty, k+1, l, x1) - sum_ceil_sqrt(beatty, k, l, x1)
x2 = max(l, v + 1)
if u >= x2:
total += (u - x2 + 1) * (n + 1) - sum_ceil_sqrt(beatty, k, x2, u)
return total
return str(count_rays(2000000, 1000000000))
if __name__ == '__main__':
print(solve())
Java
import java.math.BigInteger;
import java.util.HashMap;
import java.util.Map;
import java.util.Objects;
public class Euler372 {
static long isqrt(long n) {
if (n == 0)
return 0;
long approx = (long) Math.sqrt(n);
while ((approx + 1) * (approx + 1) <= n)
approx++;
while (approx * approx > n)
approx--;
return approx;
}
static class BeattyKey {
int d, p, q;
long n;
BeattyKey(int d, int p, int q, long n) {
this.d = d;
this.p = p;
this.q = q;
this.n = n;
}
@Override
public boolean equals(Object o) {
if (this == o)
return true;
if (o == null || getClass() != o.getClass())
return false;
BeattyKey that = (BeattyKey) o;
return d == that.d && p == that.p && q == that.q && n == that.n;
}
@Override
public int hashCode() {
int h = d;
h = h * 1315423911 + p;
h = h * 1315423911 + q;
h = h * 1315423911 + (int) (n ^ (n >>> 32));
return h;
}
}
static class BeattySumSqrt {
Map<BeattyKey, BigInteger> memo = new HashMap<>();
BigInteger sumFloor(int d, long n) {
return sumFloorImpl(d, n, 0, 1);
}
boolean geAsqrtD(long a, int d, long rhs) {
if (rhs <= 0)
return true;
BigInteger left = BigInteger.valueOf(a).multiply(BigInteger.valueOf(a)).multiply(BigInteger.valueOf(d));
BigInteger right = BigInteger.valueOf(rhs).multiply(BigInteger.valueOf(rhs));
return left.compareTo(right) >= 0;
}
long floorLinearSurd(long a, int d, long b, int c) {
double guessLd = (a * Math.sqrt(d) + b) / c;
long g = (long) Math.floor(guessLd);
while (geAsqrtD(a, d, (g + 1) * c - b))
g++;
while (!geAsqrtD(a, d, g * c - b))
g--;
return g;
}
BigInteger sumFloorImpl(int d, long n, int p, int q) {
if (n == 0)
return BigInteger.ZERO;
int s = (int) Math.sqrt(d);
if (s * s == d) {
return BigInteger.valueOf(s).multiply(BigInteger.valueOf(n)).multiply(BigInteger.valueOf(n + 1))
.divide(BigInteger.valueOf(2));
}
BeattyKey key = new BeattyKey(d, p, q, n);
if (memo.containsKey(key))
return memo.get(key);
int a = (s + p) / q;
int p1 = a * q - p;
int q1 = (d - p1 * p1) / q;
long mSigned = floorLinearSurd(n, d, -n * p1, q);
long m = mSigned;
BigInteger tri = BigInteger.valueOf(n).multiply(BigInteger.valueOf(n + 1)).divide(BigInteger.valueOf(2));
BigInteger result = BigInteger.valueOf(a).multiply(tri)
.add(BigInteger.valueOf(n).multiply(BigInteger.valueOf(m)))
.subtract(sumFloorImpl(d, m, p1, q1));
memo.put(key, result);
return result;
}
}
static BigInteger sumCeilSqrt(BeattySumSqrt beatty, int d, long l, long r) {
if (l > r)
return BigInteger.ZERO;
int s = (int) Math.sqrt(d);
if (s * s == d) {
BigInteger cnt = BigInteger.valueOf(r - l + 1);
return BigInteger.valueOf(s).multiply(BigInteger.valueOf(l + r)).multiply(cnt)
.divide(BigInteger.valueOf(2));
}
BigInteger floors = beatty.sumFloor(d, r).subtract(beatty.sumFloor(d, l - 1));
return floors.add(BigInteger.valueOf(r - l + 1));
}
static BigInteger countRays(long m, long n) {
long l = m + 1;
long kMax = (n * n) / (l * l);
BeattySumSqrt beatty = new BeattySumSqrt();
BigInteger total = BigInteger.ZERO;
for (long k = 1; k <= kMax; k += 2) {
long u = isqrt((n * n) / k);
if (u < l)
continue;
long v = isqrt(((n + 1) * (n + 1)) / (k + 1));
long x1 = Math.min(u, v);
if (x1 >= l) {
total = total.add(sumCeilSqrt(beatty, (int) (k + 1), l, x1))
.subtract(sumCeilSqrt(beatty, (int) k, l, x1));
}
long x2 = Math.max(l, v + 1);
if (u >= x2) {
total = total.add(BigInteger.valueOf(u - x2 + 1).multiply(BigInteger.valueOf(n + 1)))
.subtract(sumCeilSqrt(beatty, (int) k, x2, u));
}
}
return total;
}
static String solve() {
return countRays(2000000L, 1000000000L).toString();
}
public static void main(String[] args) {
System.out.println(solve());
}
}