# Challenge

Determine how many integer lattice points there are in an ellipse

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} \leq 1$$

centered at the origin with width $$\2a\$$ and height $$\2b\$$ where integers $$\a, b > 0\$$ .

Input

The Semi-major $$\a\$$ and Semi-minor $$\b\$$ axes.

Output

Number of interior and boundary points.

# Example

Ellipse plot showing $$\a=5\$$ and $$\b=3\$$ with $$\41\$$ blue interior and $$\4\$$ red boundary points.

Input

$$\5\$$,$$\3\$$

Output

$$\41\$$,$$\4\$$

# Test Cases

a b Interior Boundary
5 3 41 4
5 15 221 12
8 5 119 4
8 1 15 4
9 15 417 4
15 15 697 12
20 20 1245 12
• Is there a closed form?
– qwr
Commented Jan 18 at 0:19
• I don't think there's a closed form, I got the idea from the Gauss Circle Problem which there's no formula to find the number of lattice points and involves some enumeration. Commented Jan 18 at 0:23
• I think it's more interesting to have separate questions for interior points and boundary points, as boundary points is a Pythagorean triple like $(bx)^2 + (ay)^2 = (ab)^2$.
– qwr
Commented Jan 18 at 0:24
• For the Gauss circle problem, the number of interior and boundary points can be computed with Jacobi's two square theorem and the factorization of the radius length.
– qwr
Commented Jan 18 at 0:26
• Oh, I'll have to take a closer look and hopefully try to understand it. ;) Thanks. Commented Jan 18 at 0:38

# R, 76 bytes

function(a,b,z=c(a,b)^-2%*%t(expand.grid(-a:a,-b:b)^2))c(sum(z<1),sum(z==1))


Try it online!

May have issues due to numerical precision. Naive algorithm: loop over all lattice points in the $$\2a \times 2b\$$ grid centered around the origin, and count the points with the requisite relationship.

# R, 85 bytes

function(a,b,p=a^2*b^2,z=c(b,a)^2%*%t(expand.grid(-a:a,-b:b)^2))c(sum(z<p),sum(z==p))


Try it online!

More precise.

# R, 68 66 bytes

\(a,b,y=b/a*(a^2-(-a:a)^2)^.5)c(B<-sum(!y%%1)-1,sum(y%/%1+.5)-B)*2


Attempt This Online!

Total points of boundary plus interior is (up to precision error) $$\sum_{-a\le x \le a} \# \left\{ y \in \mathbb Z : |y| \le (b/a) \sqrt{a^2-x^2} \right\}$$

$$= \sum_{-a \le x \le a} \left( 2 \left\lfloor (b/a) \sqrt{a^2-x^2} \right\rfloor + 1 \right)$$

$$\B\$$ counts boundary points, subtracting over-counted cases $$\|y| = a\$$.

y=b*(1-(-a:a/a)^2)^.5 actually fails the test cases due to precision errors.

Floating point tricks in R

-2 bytes by rearranging some terms, credit Giuseppe.

• You can save another 2 bytes by reversing the order of output and a little manipulation; Attempt it online Commented Jan 18 at 19:00

# Jelly, 13 bytes

rNŒp²÷²§’ṠĠẈṖ


Try it online!

A monadic link that takes a pair of integers and returns a pair of integers. Works perfectly for all of the examples, but there is a theoretical risk of floating point inaccuracy.

## Explanation

rN            | Range from a to -a and b to -b
Œp          | Cartesian product
²         | Squared
÷²       | Divide by a ** 2, b ** 2
§      | Sum inner lists
’     | Subtract 1
Ṡ    | Signs
Ġ   | Indices of like items in groups
Ẉ  | Lengths
Ṗ | Remove last (those outside the border)


## Alternative Jelly, 16 bytes

rNŒp²Uḋ²_²P$ṠĠẈṖ  Try it online! This version only uses integer arithmetic so should always be accurate even for large inputs (though like the other answer has $$\O(m \times n)\$$ complexity). • @JonathanAllan thanks. The test cases were confusing me! Commented Jan 16 at 23:29 # Google Sheets, 84 bytes Assuming $$\a\$$ is in A1 and $$\b\$$ is in B1. =SORT(COUNTIF(SEQUENCE(2*A1+1,1,-A1)^2/A1^2+SEQUENCE(1,2*B1+1,-B1)^2/B1^2,{"<1",1}))  92 bytes =SORT(COUNTIF(MAKEARRAY(2A1+1,2B1+1,LAMBDA(x,y,SUMSQ((x-A1-1)/A1,(y-B1-1)/B1))),{"<1",1})) 93 bytes =LET(a,A1,b,B1,SORT(COUNTIF(SEQUENCE(2a+1,1,-a)^2/a^2+SEQUENCE(1,2b+1,-b)^2/b^2,{"<1",1}))) # PARI/GP, 64 bytes f(a,b)=[1+vecsum(v=2*Vec(qfrep([a,0;0,b]^2,a^2*b^2)))-t=v[#v],t]  Attempt This Online! Yes, PARI/GP has a built-in for this. qfrep(q,B,{flag=0}): vector of (half) the number of vectors of norms from 1 to B for the integral and definite quadratic form q. If flag is 1, count vectors of even norm from 1 to 2B. # Google Sheets, 88 bytes =frequency(sort(sequence(2*A1+1,1,-A1)^2/A1^2+sequence(1,2*B1+1,-B1)^2/B1^2),{1-9^-9,1})  Put $$\a\$$ in cell A1, $$\b\$$ in cell B1 and the formula in cell C1. The output is in cells C1:C2 with residue in C3. Generates integer $$\x\$$ and $$\y\$$ sequences for lattice points in a rectangle that bounds the ellipse and calculates the distance from origin at each point. Uses frequency() to find points within and on the ellipse. The residue shows points within the bounding rectangle but outside the ellipse. Uses sort() but as an array enabler only. Earlier proof-of-concept version to confirm that this is easily solvable in a spreadsheet (160 bytes golfed, 279 ungolfed): =let(a, A1, b, B1, r, map(sequence(2 * a + 1, 1, -a), lambda(x, map(sequence(1, 2 * b + 1, -b), lambda(y, let( d, x^2 / a^2 + y^2 / b^2, ifs( d = 1, 2, d < 1, 1, 1, 0 ) ))) )), { countif(r, 1), countif(r, 2) } )  • -8 =let(a,A1,b,B1,r,tocol(map(sequence(2*a+1,1,-a),lambda(x,map(sequence(1,2*b+1,-b),lambda(y,let(c,x^2/a^2+y^2/b^2,ifs(c=1,2,c<1,1))))))),sort(countif(r,{1,2}))) – z.. Commented Jan 15 at 23:56 • @z.. thanks. This was just to confirm that a simple algorithm gets the correct result. The formula can be golfed quite a bit (as your answer shows.) Commented Jan 16 at 9:05 # 05AB1E, 23 22 bytes Ýí€ûnR**âOIPns-.±{Åγ¦  Outputs the pair in reversed order, so $$\[boundaryCount, interiorCount]\$$. Explanation: Ý # Convert both values in the (implicit) input-pair to [0,value]-ranged # lists í # Reverse each inner list €û # Palindromize each inner list n # Square each inner value R # Reverse the list ** # Multiply each to the (implicit) input-pair twice  # Pop and push both lists separately to the stack â # Create all possible pairs using the cartesian product O # Sum each inner pair together IP # Push the product of the input-pair n # Square that s- # Subtract the sums of the list from that value .± # Get the sign (-1 if negative; 0 if 0; 1 if positive) of each { # Sort these -1s, 0s, and 1s Åγ # Then run-length encode it, pushing list [-1,0,1] (which we won't use) # and the list of lengths ¦ # Remove the first value (for the amount of -1s) # (after which the resulting pair is output implicitly)  Try it online with step-by-step debug lines. # Charcoal, 46 bytes ＮθＮηＦＸ×η…·±θθ²ＦＸ×θ…·±ηη²⊞υ⁻Ｘ×θη²⁺ικＩ⟦ΣＥυ›ι⁰№υ⁰  Try it online! Link is to verbose version of code. Explanation: Similar approach to @vengy's C answer. ＮθＮη  Input a and b. ＦＸ×η…·±θθ²  Loop x from -a to a inclusive. ＦＸ×θ…·±ηη²  Loop y from -b to b inclusive. ⊞υ⁻Ｘ×θη²⁺ικ  Calculate the squared distance in from the perimeter, multiplied by a²b². Ｉ⟦ΣＥυ›ι⁰№υ⁰  Output the number of points for where this is positive and the number for where this is zero. # Wolfram Language(Mathematica), 135 106 bytes Saved 29 bytes thanks to @vengy Here a and b both are positive integers. Golfed version. Try it online! f[a_,b_]:=Module[{I=0,B=0},Do[With[{t=x^2/a^2+y^2/b^2},If[t<1,I++,If[t==1,B++]]],{x,-a,a},{y,-b,b}];{I,B}]  Ungolfed version. Try it online! CountLatticePoints[a_Integer, b_Integer] := Module[ {interiorPoints = 0, boundaryPoints = 0}, Sum[ If[IntegerQ[x] && IntegerQ[y], If[x^2/a^2 + y^2/b^2 < 1, interiorPoints++, If[x^2/a^2 + y^2/b^2 == 1, boundaryPoints++] ] ], {x, -a, a}, {y, -b, b} ]; {interiorPoints, boundaryPoints} ]  # Python 3.8 (pre-release), 110 bytes lambda a,b,q=2:[sum(0>(t:=x*x*b*b+y*y*a*a-a*a*b*b)or(q:=q+(t<1))*0for x in range(-a,a)for y in range(-b,b)),q]  Try it online! Consider four quarter • Oops...my inefficient answer considered four quarters. Commented Jan 15 at 23:25 • – enzo Commented Jan 15 at 23:56 • @vengy 4x slower doesn't matter at all. Seems yours is shorter 4 quarter – l4m2 Commented Jan 16 at 0:29 # Ruby, 87 bytes ->a,b{*r=0,0,w=a*b;(k=*-w..w).product(k).map{|x,y|r[w*w<=>x*x+y*y]+=1[x%a+y%b]};r[0,2]}  Try it online! # Python 2, 96 bytes lambda a,b:divmod(sum((4j,4,0)[cmp(a*b,abs(x/a*a+x%a*b*1j))]for x in range(a*b+1))-2*a-2*b+1,1j)  Try it online! Returns boundary,interior. Test bed borrowed from @l4m2 and modified. Testing code prettifies integer valued complex numbers into actual integers. ## How? Uses complex numbers for two purposes: 1. abs with complex argument for Euclidean norm 2. as a means of 2d vector summation without importing any libraries. Special Python 2 semantics of div and mod for complex operands is used to recover real and imaginary parts after summing. # Pyth, 20 bytes hMlBPx1Sm.acVdQ*Fm}_  Try it online! Takes input in the form [a, b] and outputs as [interior, boundary]. ### Explanation  # implicitly assign Q = eval(input()) hMlBPx1Sm.acVdQ*Fm}_ddQ # implicitly add ddQ m Q # map lambda d over Q }_dd # inclusive_range(-d, d) *F # cartesian product of the two ranges m # map over lambda d cVdQ # vectorized divide d by Q .a # vector size S # sort x1 # get list of indices where 1 appears P # remove the last element lB # bifurcate over length hM # map to head, this will get the first element for the list and add one for the integer  # C (gcc), 133 121 bytes -12 thanks to @Neil i,x,y,l,B;e(a,b){i=B=0;for(x=-a;x<=a;x++)for(y=-b;y<=b;y++)l=a*a*b*b-x*x*b*b-y*y*a*a,l?i+=l>0:B++;printf("%d %d\n",i,B);}  Try it online! • 121 bytes – Neil Commented Jan 16 at 0:35 • 117 Commented Jan 16 at 11:43 • @AZTECCO place l?i+=l>0:B++ into for and 0 to 115 – l4m2 Commented Jan 16 at 12:31 • Building on @l4m2 113 bytes Commented Jan 16 at 16:16 # APL+WIN, 65 bytes Prompts for a followed by b and outputs boundary the interior: ((1+2×a)++/2×⌊n)-⎕←2×¯1++/0=1|n←⎕×(1-(((⌽⍳a),0,⍳a)*2)÷(a←⎕)*2)*.5  Try it online! Thanks to Dyalog Classic # Uiua, 43 41 bytes Returns the number of interior points on the top of the stack, and the number of boundary points on the bottom. ∩/+<1:=1./+ⁿ2÷:-⍉☇1⇡+1×2..⊂  Try It ∩/+<1:=1./+ⁿ2÷:-⍉☇1⇡+1×2..⊂ -⍉☇1⇡+1×2..⊂ # generates a point array spanning the [-a,a]×[-b,b] rectangle, # we also leave a copy of a and b to divide the coordinates later ÷: # divide each x coordinate by a and each y by b ⁿ2 # square each coordinate /+ # sum each coordinate pair . # create a copy of the sums so we can check two conditions (interior and boundary) =1 # boundary check <1: # the interior point check ∩/+ # since the checks return boolean masks, sum them to get a point count.  • Try It Commented Jan 18 at 22:14 • Oh, completely missed that feature, I'll edit my submission with your link! Commented Jan 18 at 22:21 • Also, most of the Uiua code has a section that details how it works. For example Commented Jan 18 at 22:26 • Removing the trailing fix (¤), since join doesn't need array args. -2 bytes. Commented Jan 19 at 4:36 # Perl 5, 94 bytes sub{($a,$b,@r)=@_;map{//;$r[1+($b*$b*$'*$'+$a*$a*$_*$_<=>$a*$a*$b*$b)]++for-$b..$b}-$a..$a;@r}


Try it online!

Spaceship operator <=> comes in handy.

# APL(NARS), 40 chars

+/¨0(>,=)⊂{+/¯1,2*⍨⍵÷k}¨,(⊂k+1)-⍳1+2×k←⎕


How to use and test:

      +/¨0(>,=)⊂{+/¯1,2*⍨⍵÷k}¨,(⊂k+1)-⍳1+2×k←⎕
⎕:
5 3
41 4

{+/¨0(>,=)⊂{+/¯1,2*⍨⍵÷k}¨,(⊂k+1)-⍳1+2×k←⍵}¨(5 3)(5 15)(8 5)(8 1)(9 15)(15 15)(20 20)
41 4  221 12  119 4  15 4  417 4  697 12  1245 12



# Excel VBA, 93 Bytes

An immediate window function that takes input $$\a\$$ from [A1] and $$\b\$$ from [B1] of the activesheet object and outputs to the immediate window console.

Assumes a clean console environment. Maybe re-run after executing u=0:v=0.

Note: This solution has been restricted to 32-Bit versions of Excel VBA as ^ is the LongLong type literal in 64-Bit versions

For x=[-A1]To[A1]:For y=[-B1]To[B1]:t=x^2/[A1^2]+y^2/[B1^2]:u=u-(t<1):v=v-(t=1):Next y,x:?u,v


### Commented Version

For x=[-A1]To[A1]            '' iter over x axis
For y=[-B1]To[B1]          '' iter over y axis
t=x^2/[A1^2]+y^2/[B1^2]  '' calc lattice point sumsq
u=u-(t<1)                '' incr. internal point count if needed
v=v-(t=1)                '' incr. boundary point count if needed
Next y,x                     '' loop
?u,v                         '' print point counts