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 \$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

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 \$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.

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.

Answer 3

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

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

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

ＮθＮηＦＸ×η…·±θθ²ＦＸ×θ…·±ηη²⊞υ⁻Ｘ×θη²⁺ικＩ⟦ΣＥυ›ι⁰№υ⁰

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.

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

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

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 13

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 14

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!

Answer 15

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

Answer 16

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.

Answer 17

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 18

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 
 

Answer 19

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