Problem 372: Pencils of Rays

Detailed mathematical approach

Problem Summary

For integers \(0\le M\lt N\), the case the problem always uses, let \(R(M,N)\) count the lattice points \((x,y)\) that satisfy \(M\lt x\le N\) and \(M\lt y\le N\) and for which the integer

$$\left\lfloor\frac{y^2}{x^2}\right\rfloor$$

is odd; here \(\lfloor t\rfloor\) is the greatest integer not exceeding \(t\). The statement supplies two checkpoints, \(R(0,100)=3019\) and \(R(100,10000)=29750422\), and asks for \(R(2\cdot10^6,10^9)\). ^[1] Throughout \(\ell=M+1\) and \(n=N\), so both coordinates run over \(\ell,\dots,n\); these are the l and n of the solution files.

The whole solution in one sentence: the condition says the slope \(y/x\) lies between \(\sqrt k\) and \(\sqrt{k+1}\) for an odd \(k\), so the points wanted fill a fan of angular sectors, and one sector costs only a few sums \(\sum\lceil x\sqrt d\rceil\), each of which a Euclidean-style recursion evaluates in a handful of steps.

Mathematical Approach

1. What the floor condition really says

Since \(x\) and \(y\) are positive, \(\lfloor y^2/x^2\rfloor=k\) holds exactly when

$$k\le\frac{y^2}{x^2}\lt k+1,$$

and taking square roots turns this into a statement about the slope of the ray from the origin through \((x,y)\):

$$\sqrt k\le\frac yx\lt\sqrt{k+1}.$$

So the rays \(y=x\sqrt k\), \(k=0,1,2,\dots\), cut the quadrant \(x\gt0\), \(y\ge0\) into angular sectors, and the sector with index \(k\) carries the floor value \(k\). We keep the odd sectors and discard the even ones — a pencil of rays through the origin, which is where the problem takes its name.

An odd \(k\) is at least \(1\), so every counted point has \(y\ge x\): the odd sectors lie on or above the diagonal, which is itself counted in full.

The floor values for \(1\le x,y\le5\) make the pattern visible: the rows \(x=1,\dots,5\) read \(1,4,9,16,25\); \(0,1,2,4,6\); \(0,0,1,1,2\); \(0,0,0,1,1\); \(0,0,0,0,1\). The odd entries sit at

$$(1,1),(1,3),(1,5),(2,2),(3,3),$$

$$(3,4),(4,4),(4,5),(5,5),$$

so \(R(0,5)=9\). The row \(x=1\) alone already contributes three odd values, \(1\), \(9\) and \(25\): small \(x\) sees many sectors, and that is the whole difficulty.

2. Grouping the points by sector

For odd \(k\ge1\) put

$$A_k=\{(x,y):\ell\le x,y\le n,\ \lfloor y^2/x^2\rfloor=k\}.$$

Each point has exactly one floor value, so the \(A_k\) are disjoint and

$$R(M,N)=\sum_{j\ge0}|A_{2j+1}|.$$

The sum is finite: \(x\ge\ell\) and \(y\le n\) give \(y^2/x^2\le n^2/\ell^2\), so every non-empty \(A_k\) has

$$k\le K=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor.$$

This is the step that makes the problem tractable, and the reason the statement starts the square at \(M=2\cdot10^6\) rather than \(0\), where the bound would be \(K=n^2=10^{18}\). With \(\ell=2000001\) and \(n=10^9\) we computed \(K=249999\), so there are \(125000\) odd sectors, while the square itself holds \(998000000^2\approx9.96\cdot10^{17}\) lattice points.

3. How many \(y\) belong to one sector and one \(x\)

One lemma is needed. For reals \(a\le b\) the integers in \([a,b)\) are \(\lceil a\rceil,\dots,\lceil b\rceil-1\), since an integer \(j\) has \(j\ge a\) exactly when \(j\ge\lceil a\rceil\), and \(j\lt b\) exactly when \(j\le\lceil b\rceil-1\). Hence

$$\#\bigl([a,b)\cap\mathbb Z\bigr)=\lceil b\rceil-\lceil a\rceil.$$

Now fix an odd \(k\) and an \(x\) with \(\ell\le x\le n\). By Section 1 the admissible \(y\) are the integers with

$$x\sqrt k\le y\lt x\sqrt{k+1},$$

subject to \(\ell\le y\le n\). The lower bound costs nothing: \(k\ge1\) gives \(y\ge x\sqrt k\ge x\ge\ell\), so \(M\lt y\) never has to be enforced; the upper one is \(y\lt n+1\). Applying the lemma to \([x\sqrt k,\ \min(x\sqrt{k+1},n+1))\), with \(\lceil\min(t,m)\rceil=\min(\lceil t\rceil,m)\) for integer \(m\), the count is

$$C_k(x)=\max\bigl(0,\,U_k(x)-\lceil x\sqrt k\rceil\bigr),$$

$$U_k(x)=\min\bigl(\lceil x\sqrt{k+1}\rceil,\,n+1\bigr).$$

The outer maximum matters: the lemma needs ordered endpoints, and for large \(x\) the lower endpoint passes the upper one. With \(n=3\), \(k=3\), \(x=3\) the subtraction alone gives \(4-6=-2\), where the true count is \(0\).

4. Splitting the \(x\)-range in two

The count \(C_k(x)\) has two shapes, depending on which term the minimum picks. Define

$$u_k=\left\lfloor\frac n{\sqrt k}\right\rfloor,\qquad v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor,$$

so \(u_k\) is the largest integer \(x\) with \(x\sqrt k\le n\), and \(v_k\) the largest with \(x\sqrt{k+1}\le n+1\). Three facts follow.

• If \(x\gt u_k\) then \(x\sqrt k\gt n\), so no admissible \(y\) exists and \(C_k(x)=0\): only \(x\le u_k\) contributes.
• If \(x\le v_k\) then \(U_k(x)=\lceil x\sqrt{k+1}\rceil\), so \(C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt k\rceil\), already non-negative as a count of integers in an interval.
• If \(v_k\lt x\le u_k\) then \(U_k(x)=n+1\) and \(C_k(x)=(n+1)-\lceil x\sqrt k\rceil\), at least \(1\) because \(\lceil x\sqrt k\rceil\le n\).

So, with \(x_1=\min(u_k,v_k)\) and \(x_2=\max(\ell,v_k+1)\), the interval \([\ell,u_k]\) splits into \([\ell,x_1]\), where the second bullet applies, and \([x_2,u_k]\), where the third does. The pieces are disjoint and cover \([\ell,u_k]\): if \(v_k\ge u_k\) the second is empty, if \(v_k\lt\ell\) the first is, and otherwise the first ends at \(v_k\) and the second starts at \(v_k+1\). On both the clamp of Section 3 is inactive, so no maximum survives below. Finally \(u_k\ge\ell\) for every \(k\le K\), since \(k\le n^2/\ell^2\) gives \(\ell\sqrt k\le n\), so the guard that would skip a whole sector never fires — checked for all \(125000\) sectors of the target.

5. Everything reduces to sums of \(\lceil x\sqrt d\rceil\)

Write

$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt d\rceil,$$

with \(T_d(a,b)=0\) when \(a\gt b\); the closed forms below are for \(a\le b\). Let \(h_k=\max(0,u_k-x_2+1)\) be the length of the second piece. Summing \(C_k(x)\) over both pieces,

$$\begin{aligned}|A_k|={}&T_{k+1}(\ell,x_1)-T_k(\ell,x_1)\\&+h_k(n+1)-T_k(x_2,u_k).\end{aligned}$$

The maximum is needed: with \(\ell=1\), \(n=3\), \(k=3\) we get \(u_k=1\), \(x_2=3\), and \(u_k-x_2+1=-1\) would charge the empty piece \(-4\) instead of \(0\). At most three interval sums per sector, and nothing else: counting over \(10^{18}\) lattice points becomes \(374998\) evaluations of \(T_d\), two short of \(3\cdot125000\) because the first piece is empty for \(k=249999\).

6. Doing it in exact integers

Every quantity above involves \(\sqrt k\), yet all are integers, and rounding a floating-point square root the wrong way at a boundary would corrupt a whole sector. Two identities remove floating point. First, for positive integers \(t,d\),

$$\left\lfloor\frac t{\sqrt d}\right\rfloor=\left\lfloor\sqrt{\left\lfloor\frac{t^2}d\right\rfloor}\right\rfloor,$$

because for a non-negative integer \(x\) — squaring needs the sign —

$$\begin{aligned}x\sqrt d\le t&\iff x^2d\le t^2\\&\iff x^2\le\lfloor t^2/d\rfloor,\end{aligned}$$

the last step being legitimate because \(x^2\) is an integer. Both \(u_k\) and \(v_k\) come this way from an integer square root. ^[2] Second,

$$\lceil x\sqrt d\rceil=\bigl\lfloor\sqrt{x^2d}\bigr\rfloor+\varepsilon,$$

with \(\varepsilon=0\) if \(x^2d\) is a perfect square and \(1\) if not — the fact that \(\lceil t\rceil-\lfloor t\rfloor\) is \(0\) for integral \(t\) and \(1\) otherwise. ^[3]

7. Square and non-square \(d\)

If \(d=s^2\) is a square then \(\lceil x\sqrt d\rceil=sx\) and the interval sum collapses to an arithmetic series,

$$T_d(a,b)=s\cdot\frac{(a+b)(b-a+1)}2.$$

Among the \(125000\) odd \(k\le K\) of the target, \(250\) are squares and \(250\) have \(k+1\) square. A square \(k\) takes the shortcut in both pieces, \(500\) calls; a square \(k+1\) in one, adding \(249\) rather than \(250\), since the first piece missing at \(k=249999\) is just the one whose successor \(250000\) is square. The shortcut fires \(749\) times.

Otherwise \(\sqrt d\) is irrational, and so is \(x\sqrt d\) for every \(x\ge1\). Hence \(\lceil x\sqrt d\rceil=\lfloor x\sqrt d\rfloor+1\) exactly, and with the prefix sum

$$F_d(m)=\sum_{x=1}^{m}\lfloor x\sqrt d\rfloor$$

we get

$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$

The sequence \(\lfloor x\sqrt d\rfloor\) is the Beatty sequence of \(\sqrt d\). ^[4] Everything now rests on computing \(F_d(m)\) fast: a direct loop costs \(m\) steps, and \(m\) reaches \(10^9\).

8. A generalised floor sum and the surd step

Fix a non-square \(d\gt0\). The recursion needs a larger family: for integers \(p\), \(q\gt0\) put

$$F_{d,p,q}(m)=\sum_{i=1}^{m}\left\lfloor i\,\frac{\sqrt d+p}q\right\rfloor,$$

so \(F_d(m)=F_{d,0,1}(m)\). Put \(\alpha=(\sqrt d+p)/q\) and \(s=\lfloor\sqrt d\rfloor\), and split off the integer part: because \(\lfloor(\lfloor t\rfloor+p)/q\rfloor=\lfloor(t+p)/q\rfloor\), ^[3]

$$a=\lfloor\alpha\rfloor=\left\lfloor\frac{s+p}q\right\rfloor,$$

by integer division, once \(s\) is known. Put \(p_1=aq-p\); the fractional part is then

$$\beta=\alpha-a=\frac{\sqrt d-p_1}q\in(0,1),$$

and rationalising the denominator of \(1/\beta\) gives

$$\begin{aligned}\frac1\beta&=\frac q{\sqrt d-p_1}=\frac{\sqrt d+p_1}{q_1},\\q_1&=\frac{d-p_1^2}q.\end{aligned}$$

So the reciprocal of the fractional part has the same shape as \(\alpha\), with new parameters \((p_1,q_1)\). That \(q_1\) is an integer follows from the invariant \(q\mid d-p^2\), which holds at the start, where \(p=0\), \(q=1\), and is preserved because

$$d-p_1^2=(d-p^2)+2apq-a^2q^2$$

is divisible by \(q\) whenever \(d-p^2\) is. This map is exactly the classical step of the continued-fraction expansion of \(\sqrt d\). ^[5]

Positivity of \(q_1\) needs an argument, and brings bounds with it. Write \(\alpha'\) for the conjugate of a surd, got by replacing \(\sqrt d\) with \(-\sqrt d\). The first step turns \((0,1)\) into \(1/\beta=(\sqrt d+s)/(d-s^2)\gt1\), whose conjugate \(-1/(\sqrt d+s)\) lies in \((-1,0)\); such a surd is called reduced. Reducedness propagates: if \(\alpha\gt1\) and \(\alpha'\in(-1,0)\) then \(a=\lfloor\alpha\rfloor\ge1\), the next value \(1/(\alpha-a)\) exceeds \(1\), and its conjugate \(1/(\alpha'-a)\) lies in \((-1,0)\) because \(\alpha'-a\lt-1\). A reduced value minus its conjugate is \(2\sqrt d/q_1\), and positive, so \(q_1\gt0\); the two inequalities read

$$\sqrt d-p\lt q\lt\sqrt d+p,\qquad p\lt\sqrt d,$$

giving \(0\lt p\lt\sqrt d\) and \(0\lt q\lt2\sqrt d\) at every later state. With \(d\le250000\) that means \(p\le499\) and \(q\le999\), and a scan of all non-square \(d\le250000\) attains \(p=499\) and \(q=998\): well inside 32-bit integers.

9. The reflection identity

Here is the step that makes \(F_{d,p,q}\) cheap. First the idea on four columns, with \(\beta=\sqrt2-1\) and \(m=4\): the column counts \(\lfloor i\beta\rfloor\) are \(0,0,1,1\), two points in all, and the only occupied row is \(j=1\), holding the columns \(i=3,4\) because \(\lfloor1/\beta\rfloor=2\) — counting by rows gives \(4-2=2\) again. In general, let \(\beta\in(0,1)\) be irrational, \(m\ge1\) and \(r=\lfloor m\beta\rfloor\). Then

$$\sum_{i=1}^{m}\lfloor i\beta\rfloor=mr-\sum_{j=1}^{r}\left\lfloor\frac j\beta\right\rfloor.$$

Proof. Both sides count the lattice points \((i,j)\) with \(1\le i\le m\) and \(1\le j\le i\beta\); by columns one gets the left side. By rows: \(j\) occurs only for \(j\le m\beta\), that is \(j\le r\), and then \(j\le i\beta\) reads \(i\ge j/\beta\). As \(\beta\) is irrational, \(j/\beta\) is not an integer, so the admissible \(i\) are \(\lfloor j/\beta\rfloor+1,\dots,m\), which is \(m-\lfloor j/\beta\rfloor\) of them. Summing over \(j\) gives the right side. \(\square\)

Because \(\lfloor i\alpha\rfloor=ai+\lfloor i\beta\rfloor\) and \(1/\beta=(\sqrt d+p_1)/q_1\), the identity becomes a recursion whose limit shrinks:

$$F_{d,p,q}(m)=a\frac{m(m+1)}2+mr-F_{d,p_1,q_1}(r),$$

$$r=\left\lfloor\frac{m\sqrt d-mp_1}q\right\rfloor.$$

Since \(\beta\lt1\) we have \(r\le m-1\), so the recursion terminates; the complexity section proves that the limit in fact halves every two levels. Each level needs one more exact floor \(\lfloor(m\sqrt d+b)/q\rfloor\), which the code obtains from a floating-point guess corrected by exact integer comparisons.

10. Worked example: \(R(2,8)=16\)

Take \(M=2\), \(N=8\): then \(\ell=3\), \(n=8\), \(K=\lfloor64/9\rfloor=7\), and the odd sectors are \(k=1,3,5,7\).

For \(k=1\): \(u_1=8\), \(v_1=\lfloor9/\sqrt2\rfloor=6\), so \(x_1=6\), \(x_2=7\). On \(x=3,\dots,6\) the ceilings \(\lceil x\sqrt2\rceil\) are \(5,6,8,9\) against \(\lceil x\rceil=x\), giving counts \(2,2,3,3\) with \(y\) sets \(\{3,4\},\{4,5\},\{5,6,7\},\{6,7,8\}\) — \(10\) points. On \(x=7,8\) the count is \(9-x\), so \(2\) and \(1\), from \(y\in\{7,8\}\) and \(y=8\). Sector total \(13\).

For \(k=3\): \(u_3=\lfloor8/\sqrt3\rfloor=4\) and \(v_3=\lfloor9/2\rfloor=4\), so \(x_1=4\) and the second piece is empty. Here \(k+1=4\) is a square, so \(\lceil2x\rceil=2x\) and the counts are \(6-\lceil3\sqrt3\rceil=0\) and \(8-\lceil4\sqrt3\rceil=1\): the point \((4,7)\), and indeed \(49/16=3.0625\).

For \(k=5\): \(u_5=v_5=3\), only \(x=3\) survives, and \(\lceil3\sqrt6\rceil-\lceil3\sqrt5\rceil=8-7=1\), the point \((3,7)\) with \(49/9\approx5.44\). For \(k=7\): again \(u_7=v_7=3\), and \(\lceil3\sqrt8\rceil-\lceil3\sqrt7\rceil=9-8=1\), the point \((3,8)\) with \(64/9\approx7.11\).

So \(13+1+1+1=16\), and a direct scan of the \(36\) pairs confirms those points. The later sectors contribute only three, all with small \(x\) and large \(y\) — the thin high-slope slivers.

11. Worked example: the recursion for \(F_2(10)\)

Take \(d=2\), \(p=0\), \(q=1\), \(m=10\), so \(s=1\), \(a=1\), \(p_1=1\), \(q_1=(2-1)/1=1\) and \(\beta=\sqrt2-1\approx0.41421\). Then \(r=\lfloor10\beta\rfloor=4\) and

$$F_{2,0,1}(10)=1\cdot55+10\cdot4-F_{2,1,1}(4).$$

At the next level \(a=2\), the parameters repeat as \(p_1=q_1=1\), and \(r=\lfloor4\beta\rfloor=1\), so

$$F_{2,1,1}(4)=2\cdot10+4\cdot1-F_{2,1,1}(1)=22,$$

with \(F_{2,1,1}(1)=2\cdot1+1\cdot0-0=2\). So \(F_2(10)=55+40-22=73\), and summing \(\lfloor x\sqrt2\rfloor\) for \(x=1,\dots,10\) — that is \(1,2,4,5,7,8,9,11,12,14\) — indeed gives \(73\). ^[6]

The example also shows why the implementations cache on \((d,p,q,m)\): the pair \((p,q)=(1,1)\) recurs, because the continued fraction of \(\sqrt2\) is periodic, and by Lagrange's theorem so is that of every \(\sqrt d\). ^[5] Only the limit really changes down a chain, which is what makes the cache so thin a win below.

12. A continuous cross-check

The odd sectors form a cone, so their share of a large square follows from area. At height \(y\in(0,1]\) the sector \(k\) spans \(x\) from \(y/\sqrt{k+1}\) to \(y/\sqrt k\), a width proportional to \(y\); as \(\int_0^1y\,dy=\tfrac12\), its area in the unit square is \(\tfrac12(k^{-1/2}-(k+1)^{-1/2})\). Adding over odd \(k\) regroups into the alternating series

$$\frac12\sum_{j=1}^{\infty}\frac{(-1)^{j-1}}{\sqrt j}=\frac{\eta(1/2)}2,$$

with \(\eta\) the Dirichlet eta function, \(\eta(s)=(1-2^{1-s})\zeta(s)\), and \(\zeta\) the Riemann zeta function. ^[7] The fraction is \(0.3024493\dots\), close to \(3019/100^2=0.3019\). A shifted square needs the clipped area, not the cone density: restricting each interval to \((100,10000]\) before integrating gives \(29746832.70\), about \(0.0121\) per cent from the count \(29750422\). That is a plausibility check, not a proof.

How the Code Works

The C++ file is the reference implementation. ^[8] isqrt_u128 walks a long double estimate until \(r^2\le n\lt(r+1)^2\) holds exactly, so the estimate need only be close, and floor_div_by_sqrt(t, d) applies Section 6's identity as isqrt_u128(t*t/d) in 128 bits, supplying \(u_k\) and \(v_k\).

count_rays(m, n) is the outer loop of Section 2. It sets l = m+1 and k_max = (n*n)/(l*l) in 128-bit arithmetic, then walks \(k=1,3,5,\dots\), computing u and v, forming x1 = min(u,v) and x2 = max(l, v+1), and adding each piece of Section 5 only when its range is non-empty, which is how it avoids the negative terms warned about in Sections 3 and 5. The accumulator is an unsigned __int128, a GCC extension, ^[9] printed by to_string_u128 through repeated division by ten. sum_ceil_sqrt(beatty, d, l, r) is \(T_d(l,r)\): it compares s*s with d and either returns the arithmetic series of Section 7 or asks for two prefix floor sums and adds the term count.

The class BeattySumSqrt holds the recursion of Sections 8 and 9. Its public sum_floor(d, n) calls the private sum_floor_impl(d, n, p, q) with \(p=0\), \(q=1\); that function returns \(0\) for \(n=0\), takes the square shortcut when s*s == d, and otherwise looks up a BeattyKey — the quadruple \((d,p,q,n)\), with its hash functor BeattyKeyHash — in the table memo_. On a miss it computes a, p1 and q1 as in Section 8, obtains \(r\) from floor_linear_surd, recurses, stores and returns. floor_linear_surd is the exact floor of \((a\sqrt d+b)/c\): a long double guess, then two correction loops driven by ge_asqrt_d, which returns true at once when its rhs argument is not positive and otherwise compares the squares in 128 bits — the guess-then-verify shape of the integer-square-root algorithms. ^[2] Lastly run_checkpoints checks count_rays(0,100) and count_rays(100,10000) against \(3019\) and \(29750422\), failing with a nonzero exit status unless main is given --skip-checkpoints.

The Python file mirrors it function for function inside one solve. It needs no wide-integer machinery, Python integers having unlimited precision, ^[10] and reaches math.isqrt through the one-line wrapper isqrt128. It still calls math.sqrt twice over, to classify a radicand as square and to seed floor_linear_surd; only the second is corrected by exact comparisons, and the first is safe because \(d\le250000\).

The Java file keeps the structure with BigInteger for every accumulated value ^[11] and a BeattyKey class with hand-written equals and hashCode so it can serve as a HashMap key. Two places differ: its isqrt works in signed 64-bit arithmetic, safe because its arguments stay below \((10^9+1)^2\), and countRays forms (n*n)/k and ((n+1)*(n+1))/(k+1) in long where the C++ widens. Both narrow the sector index to int, safe because Section 2 bounds \(k+1\) by \(250000\). Only the C++ runs the published checkpoints; we ran all three and they print the same value.

Complexity and Verification

The outer loop runs over the odd \(k\le K=\lfloor n^2/\ell^2\rfloor\), that is \(K/2+O(1)\) sectors, each costing at most three interval sums and so six prefix sums. A square radicand is \(O(1)\); a non-square one recurses, with logarithmic depth. Let \(\beta\) be the current fractional part and \(m\) the limit. If \(\beta\le1/2\) the next limit \(\lfloor m\beta\rfloor\) is at most \(m/2\); if \(\beta\gt1/2\) then \(1\lt1/\beta\lt2\), the next fractional part is \(1/\beta-1\), and two steps bring the limit below \(m\beta(1/\beta-1)=m(1-\beta)\lt m/2\). Either way it halves every two levels, so the depth is \(O(\log n)\). In arithmetic operations, with expected constant-time hashing, the total is

$$O\!\left(\frac{N^2}{(M+1)^2}\log N\right)$$

time, with \(O(K\log n)\) stored states and \(O(\log n)\) stack — counts of stored integers, to which Python and Java add the bit cost of their big integers. The target counts confirm the scale: \(125000\) sectors, \(748498\) top-level prefix requests, \(7558320\) recursive calls, about \(10.1\) per request, maximum depth \(32\).

Memory is the weak point, and the memo is the reason. At the target it retains \(6809822\) states and supplies \(316456\) hits — \(4.4\) per cent of the \(7126278\) lookups, or \(4.2\) per cent of all calls, since \(432042\) calls return at limit zero without consulting the table. The reuse is systematic but narrow: whenever both pieces of a sector are non-empty, \(x_2-1=x_1\), so the two interval sums with radicand \(k\) ask for the same \(F_k(x_1)\), one guaranteed hit in each of the \(124749\) sectors with non-square \(k\), and deeper paths occasionally meet too. Between sectors there is nothing to gain, since odd \(k\) uses radicands \(k\) and \(k+1\) and different odd \(k\) give disjoint pairs. On the machine used here the C++ run took \(2.4\) seconds and about \(743\) MB, while deleting the lookup and the insertion from a copy — for measurement only — left the printed value unchanged at \(0.31\) seconds and \(10\) MB, for \(10425737\) calls instead of \(7558320\). There the table is a net loss, though that balance is machine-dependent. Python took \(5.0\) seconds and \(1.75\) GB; Java, \(3.0\) seconds of wall time and \(4\) GB under an \(8\) GB heap.

The following were verified for this article by running code. A direct double loop over all \(1830\) pairs \((M,N)\) with \(0\le M\lt N\le60\) agreed in every case, and both published checkpoints reproduce. The recursion of Section 9 agrees with direct integer-square-root summation of \(\lfloor x\sqrt d\rfloor\) for all \(4180\) combinations of the \(380\) non-square \(d\) with \(2\le d\lt400\) and the limits \(0,1,2,3,5,10,30,100,500,1000,5000\), with \(q\mid d-p^2\) and \(q\gt0\) asserted at every recursive entry, and \(F_2(10)=73\) is the tenth published partial sum. ^[6] The identity \(\lceil x\sqrt d\rceil=\lfloor x\sqrt d\rfloor+1\) holds for every non-square \(d\lt200\) and \(x\lt200\); Section 10's table matches a direct scan point for point; Section 12's clipped area differs from the second checkpoint by about \(0.0121\) per cent. One structural detail: by Section 4 the test \(u_k\lt\ell\) never fires, so it is dead code rather than a safeguard.

Footnotes and References

These supply the background facts used above. The sector decomposition, the split of the \(x\)-range, the reflection identity, the surd bounds and all measured counts are derived here.

1. Project Euler 372 — Pencils of Rays. The official statement defines R(M, N) as the lattice points with M < x ≤ N, M < y ≤ N and ⌊y²/x²⌋ odd, gives R(0, 100) = 3019 and R(100, 10000) = 29750422, and asks for R(2·10⁶, 10⁹).
2. Wikipedia — Integer square root. Defines it as the greatest integer not exceeding √n, characterised by x² ≤ n < (x+1)², and corrects a candidate with while loops — the test and shape the solution files use.
3. Wikipedia — Floor and ceiling functions. States that ⌈t⌉ − ⌊t⌋ is 0 for integer t and 1 otherwise, that ⌊t + k⌋ = ⌊t⌋ + k, and that ⌊t/q⌋ = ⌊⌊t⌋/q⌋: the three rules of Sections 6 and 8.
4. Wikipedia — Beatty sequence. The Beatty sequence of an irrational r > 1 is the sequence of ⌊ir⌋, what Section 7 sums. Its own result is Rayleigh's theorem on r and r/(r−1), not Section 9's reflection identity, which is proved here.
5. Wikipedia — Periodic continued fraction. Its canonical algorithm for √S iterates m′ = d·a − m, d′ = (S − m′²)/d, a′ = ⌊(a₀ + m′)/d′⌋ — Section 8's map with (m, d) as (p, q). It also gives Lagrange's theorem on the periodicity of quadratic irrationals and Galois's reduced surd (P + √D)/Q: value above 1, conjugate in (−1, 0).
6. OEIS A194102 — partial sums of the Beatty sequence of √2. Defined as the sum of ⌊j√2⌋ for j = 1 to n, it begins 1, 3, 7, 12, 19, 27, 36, 47, 59, 73, whose tenth term is Section 11's value; the summands are A001951.
7. Wikipedia — Dirichlet eta function. Defines η(s) as the alternating series with terms (−1)ⁿ⁻¹/nˢ, convergent for positive real part, and gives η(s) = (1 − 2¹⁻ˢ)ζ(s) — how Section 12 evaluates its series.
8. C++ — Euler372.cpp. The linked immutable revision holds isqrt_u128, floor_div_by_sqrt, BeattySumSqrt, sum_ceil_sqrt, count_rays and run_checkpoints, and takes only --skip-checkpoints. It is the implementation record; the cache-free comparison and the scans above are extra validation, not tests in that file.
9. GCC manual — 128-bit integers. States that __int128 is an extension on targets with a 128-bit integer mode, with unsigned __int128 its unsigned form: the accumulator type of the C++ file.
10. Python — numeric types. States that integers have unlimited precision and that // rounds towards minus infinity, so the Python version needs no overflow analysis.
11. Java — BigInteger. Immutable arbitrary-precision integers with analogues of the primitive integer operators, used by the Java version for every accumulated sum.