Ellipse Lattice Point Counter

Question

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 \2ドルa\$ and height \2ドルb\$ 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.

ellipse

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

• qwr: Is there a closed form?
• vengy: 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.
• qwr: 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: 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.
• vengy: Oh, I'll have to take a closer look and hopefully try to understand it. ;) Thanks.

Answer 1

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 \2ドルa \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.

Answer 2

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.

• Giuseppe: You can save another 2 bytes by reversing the order of output and a little manipulation; Attempt it online

Answer 3

Jelly, 13 bytes

rNŒp2÷2§’ṠĠẈṖ

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
 2 | Squared
 ÷2 | 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Œp2Uḋ2_2P$ṠĠẈṖ

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).

• Nick Kennedy: @JonathanAllan thanks. The test cases were confusing me!

Answer 4

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}))) (削除ここまで)

Answer 5

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.

Answer 6

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) } 
)

• z..: -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})))
• doubleunary: @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.)

Answer 7

05AB1E, (削除) 23 (削除ここまで) 22 bytes

Ýí€ûnR**`âOIPns-.±{Åγ¦

Outputs the pair in reversed order, so \$[boundaryCount, interiorCount]\$.

Try it online or verify all test cases.

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.

Answer 8

Charcoal, 46 bytes

×ばつθη2+ικI⟦ΣEυ›ι0Noυ0

Try it online! Link is to verbose version of code. Explanation: Similar approach to @vengy's C answer.

NθNη

Input a and b.

×ばつη...·±θθ2

Loop x from -a to a inclusive.

×ばつθ...·±ηη2

Loop y from -b to b inclusive.

×ばつθη2+ικ

Calculate the squared distance in from the perimeter, multiplied by a2b2.

I⟦ΣEυ›ι0Noυ0

Output the number of points for where this is positive and the number for where this is zero.

Answer 9

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}
]

Answer 10

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

• vengy: Oops...my inefficient answer considered four quarters.
• enzo: -4 bytes?
• l4m2: @vengy 4x slower doesn't matter at all. Seems yours is shorter 4 quarter

Answer 11

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!

Answer 12

APL+WIN, 65 bytes

Prompts for a followed by b and outputs boundary the interior:

((×ばつa)×ばつ⌊n)×ばつ ̄1++/0=1|n×ばつ(1-(((⌽⍳a),0,⍳a)*2)÷(a←⎕)*2)*.5

Try it online! Thanks to Dyalog Classic

Answer 13

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.

Answer 14

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 

Answer 15

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

Answer 16

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!

• Neil: 121 bytes
• AZTECCO: 117
• l4m2: @AZTECCO place l?i+=l>0:B++ into for and 0 to 115
• ceilingcat: Building on @l4m2 113 bytes

Answer 17

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./+n×ばつ2..⊂

Try It

∩/+<1:=1./+n×ばつ2..⊂
 ×ばつ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
 n2 # 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.

• vengy: Try It
• Maxim Khanov: Oh, completely missed that feature, I'll edit my submission with your link!
• vengy: Also, most of the Uiua code has a section that details how it works. For example
• Maxim Khanov: Removing the trailing fix (¤), since join doesn't need array args. -2 bytes.

Answer 18

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.

Answer 19

APL(NARS), 40 chars

+/ ̈0(>,=)⊂{+/ ̄1,2*⍨⍵÷k} ̈,(⊂k+1)×ばつk←⎕

How to use and test:

 +/ ̈0(>,=)⊂{+/ ̄1,2*⍨⍵÷k} ̈,(⊂k+1)×ばつk←⎕
⎕:
 5 3
41 4 
 {+/ ̈0(>,=)⊂{+/ ̄1,2*⍨⍵÷k} ̈,(⊂k+1)×ばつ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 
 

Answer 20

64-bit Aarch64 assembly, 84 bytes / 21 insns

f:
 mov x0, xzr
 mov x1, xzr
 fneg d2, d0 // x = -a
 fneg d3, d1 // y = -b
 fmov d4, #1.0
 fmul d5, d0, d0 // a^2
 fmul d6, d1, d1 // b^2
 1:
 fmul d7, d3, d3 // y^2
 fdiv d7, d7, d6 // y^2/b^2
 fdiv d16, d2, d5 // x /a^2
 fmadd d16, d16, d2, d7 // x/a^2 * x + y^2/b^2
 fcmp d16, d4
 cinc x0, lt
 cinc x1, eq
 fadd d2, d2, d4
 fcmp d2, d0
 b.le 1b
 fneg d2, d0
 fadd d3, d3, d4
 fcmp d3, d1
 b.le 1b
 ret

This is a crude, brute-force search from (-a, -b) to (a, b) checking whether the lattice points match. It's not well-optimised, but I don't really see a better way out of it. It expects a, b as doubles in d0, d1 and (hopefully*) returns the number of points on the interior in x0 and on the perimiter in x1. Its signature would look something like struct answer { uint64_t interior; uint64_t perimiter }; extern answer f(double, double);

It should* work for a reasonable range of input numbers.

Notes:

• *I haven't tested it, but the compiler emits what looks like closely enough similar code at a glance, but it wastes instructions setting up temporary registers for the return values, and doesn't perform a fmadd fuse instead choosing to just multiply and add separately (why? this is a conceptually pretty trivial optimisation that both GCC and Clang fail to do even on -Os)