The Caged Circles

PowerShell, 317 313 305 298 281 274 231 221 210 207 195 181 136 126 111 109 108 bytes

-an enormous 94 bytes thanks to mazzy!

I must say, I never expected to beat Charcoal in this challenge; and when I posted that 317 byte answer, I certainly never expected to get it all the way down to the low 100s. I'm hopeful I can shave those final 10 to get to sub-100, but this feels so very optimized!

Takes the x and y offsets (center of the circle) as inputs. Generates all points in the range (-9..9,-9..9), then adds points for potential border intersections, which handles offset circles, and handles doubling point values on x,0 or 0,y coordinates when the circle is centered on the origin. Groups all of those points by their squared distance from the origin, then for each group outputs the sum of the group's count and the counts of all groups with distance smaller than that group. Multiplies all of the squared distances by 1e7 then casts to Integer to avoid issues with floating points.

Because it groups the distances, it does output more than the first 16 terms, but as the challenger stated that this is acceptable, the TIO link truncates the output to the first 16 terms to make it more clear which tests succeed.

$x,$y=$args|%{,(-9..9+($p=$_%1)|%{($t=$_-$p)*$t*1e7})}
$x|%{$n=$_;$y|%{[Int]($n+$_)}}|group|% c*|%{($s+=$_)}

Try it online!

Ungolfed

At this point, the ungolfed solution is more a representation of the logic behind my solution than it is an ungolfed version of the solution - as the iterations have been made, I haven't been able to keep up with the changes on the ungolfed version. I am leaving it in here, despite the discrepancy, in hopes that it will inspires to answer in their language of choice (or port it to their language)

#Take the x-offset, y-offset, and starting radius as a parameter
#radius starts at 0, so that the first element will be the count of points on the origin of the circle
param($xOffset,$yOffset,$radius=0)
#set the x and y offsets to be in the range [0,1)
$xOffset=$xOffset % 1
$yOffset=$yOffset % 1
#do 16 times (for the required 16 elements)
1..16|%{
	#make an array of points
	$points = @()
	#set our limits to the edge of the bounding box of the circle
	#this will normally be within 1 of the actual circle, but offsets of .5
	#will increase the search space due to .NET's banker's rounding
	$upperSearchLimitX = [Int]($xOffset+$radius+1)
	$lowerSearchLimitX = -$upperSearchLimitX 
	$upperSearchLimitY = [Int]($yOffset+$radius+1)
	$lowerSearchLimitY = -$upperSearchLimitY 
	#generate x and y for every combination of values in those ranges
	$lowerSearchLimitX..$upperSearchLimitX | %{
		$x = $_
		$lowerSearchLimitY..$upperSearchLimitY | %{
			$y = $_
			#add the point to our list
			$points+=@{x=$x;y=$y}
			#If x is zero, also add a special point that represents a border intersection to the list
			#this will also make x,0 points count for two when the center is at 0,0
			if (!$x) {
				$points+=@{x=$x+$xOffset;y=$y}
			}
			#If y is zero, also add a special point that represents a border intersection to the list
			#this will also make 0,y points count for two when the center is at 0,0
			if (!$y) {
				$points+=@{x=$x;y=$y+$yOffset}
			}
		}
	}
	#add the center point to the list
	$points += @{x=$xOffset;y=$yOffset}
	#set the minimum point as the closest point to our current circle that's outside of the bounds of the circle
	$minimum =	$points|?{ 
				($_.distance=[Math]::Sqrt(($_.x-$xOffset)*($_.x-$xOffset)+($_.y-$yOffset)*($_.y-$yOffset))) -gt $radius
			}|sort distance -t 1
	#set our new radius to be that point's distance from the x offset
	$radius = $minimum.distance
	#count the number of points that lie inside (but not on the boundary of) this new circle, output implicitly
	@($points|?{$_.distance -lt $radius}).Count
}

Try it online!

I'm not the best at making images; if anyone can recommend a good tool for things like these, let me know. But here's a rough sketch of all the points that get checked along with the order we iterate to them in for x=3/8, y=1/4

Charcoal, 117 bytes

ＮθＮηＮζＮε≧×εθ≧×ηζ≧×ηε≔⁰ηＷ‹Ｌυ¹⁶«≔⟦⟧δＦ…·±ηηＦ…·±ηηＦ¬›ΣＸ⟦κλ⟧²ηＦ…·÷⁺θκε±÷⁺θκ±εＦ…·÷⁺ζλε±÷⁺ζλ±εＦ¬№δ⟦νξ⟧⊞δ⟦νξ⟧Ｆ¬№υＬδ⊞υＬδ≦⊕η»Ｉυ

Try it online! Link is to verbose version of code. Takes input as numerator denominator numerator denominator. Explanation:

ＮθＮηＮζＮε

Input the co-ordinates.

≧×εθ≧×ηζ≧×ηε

Increase the scale of the grid so that the gridlines are drawn at intervals given by the product of the denominators.

≔⁰η

Start enumerating circles with radius \$ \sqrt h + ε \$.

Ｗ‹Ｌυ¹⁶«

Repeat until 16 box counts have been found.

≔⟦⟧δ

Start with an empty list of boxes.

Ｆ…·±ηηＦ…·±ηη

Search points with a distance of up to and including \$ h \$ both horizontally and vertically.

Ｆ¬›ΣＸ⟦κλ⟧²η

Skip points whose distance from the desired point is greater than \$ h \$.

Ｆ…·÷⁺θκε±÷⁺θκ±εＦ…·÷⁺ζλε±÷⁺ζλ±ε

Consider all boxes touching this point.

Ｆ¬№δ⟦νξ⟧⊞δ⟦νξ⟧

Add the box to the list of boxes if it's not already present.

Ｆ¬№υＬδ⊞υＬδ

Add the count of boxes to the list of box counts if it's not already present.

≦⊕η

Move on to the next radius.

»Ｉυ

Print all of the box counts found.

MATL, 51 50 bytes

p-9:9*tJ*!+1ei1GP*1Gp\tFF5Mth3$v-!Z}J*!+u-|X<S&Y'Ys

The input consists of an array containing the denominators followed by an array containing the numerators. The output may contain more than 16 terms.

Try it online! Or verify all test cases (with output trimmed to 16 terms for convenience)

Explanation

Let (x, y) be the input values modulo 1. The code checks the distance from the (x, y) to each square with lower left corner (a, b), where a and b are integers. To do that, 9 points of the square are tested: (a+m,b+n) with m = 0, x, 1 and n = 0, y, 1, and the minimum of the 9 distances is kept.

A grid of 19×19 squares with a, b = −9, −8, ..., 9 is generated. The distance from (x, y) to each square is computed as described in the preceding paragraph. The 361 distances are sorted and run-length encoded, and the cumulative sum of the run lengths is the desired output.

To keep floating-point issues as controlled as possible, instead of using the actual values for x, y, a, b, m, n, the code uses those values multiplied by the product of the two input denominators.

p       % Implicit input: array with the denominators. Product
-9:9    % Push [-9, -8, ... 9]
*       % Multiply, element-wise
tJ*     % Duplicate, multiply by j (imaginary unit) element-wise
!+      % Transpose, add with broadcast. This generates the 19×19 grid of
        % complex values
1e      % Linearize into a row vector
i       % Input: array with the numerators
1GP     % Push first input again. Flip it
*       % Multiply, element-wise
1Gp     % Push first input again. Product (*)
\       % Modulo, element-wise (**)
t       % Duplicate
FF      % Push [0 0]
5Mth    % Push (*) again. Duplicate, concatenate with itself horizontally
3$v     % Concatenate three elements vertically. Gives a 3×2 matrix
-       % Subtract (*) with broadcast
!       % Transpose. The result is a 2×3 matrix
Z}      % Split into two rows with 3 elements
J*      % Multiply by j (imaginary unit) element-wise
!+      % Transpose, add with broadcast. This gives a 3×3 matrix with the
        % 9 points to be considered in each square
u       % Unique. This gives a column vector with the (unique) values
-|      % Subtract with broadcast, absolute value. Gives a 9×361 matrix
X<      % Minimum of each column. This gives the 361 distances
S       % Sort
&Y'     % Run-length encoding, producing only the run lengths
Ys      % Cumulative sum. Implicit display