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)

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

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 
 

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.

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.

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 

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.

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.

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

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

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.

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

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

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.

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!

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.

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

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

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.