
MATHEMATICA

(* This gives a polar function of a "k" sides polygon with side length "sidelength" and vertical rigthmost side *)
PolarPolygonSide[sidelength_, theta_, k_] := ((sidelength/2)/Tan[Pi/k])/Cos[Mod[theta  Pi/k , 2 Pi/k]  Pi/k];
(* uncomment next to generate and plot different polygons *)
(* Manipulate[PolarPlot[PolarPolygonSide[sidelength, theta + phase, sides], {theta, 0, 2 Pi}, PlotRange > sidelength, GridLines > {Range[sidelength, sidelength] + di, Range[sidelength, sidelength] + dj}], {sidelength, 1, 10, 1}, {sides, 3, 30, 1}, {phase, 0, 2 Pi/3, 2 Pi/300}, {dj, 0, 1/2, 0.01}, {di, 0, 1/2, 0.01}] *)
(* This function gives 1 if the point of coordinates (x, y) is strictly inside a polygon given by PolarPolygonSide[sidelength, theta, sides] rotated by "phase", and 0 otherwise *)
TruePointInsidePhase[x_, y_, sidelength_, phase_, sides_] :=
Module[{theta},
theta = ArcTan[x, y] + phase;
If[x^2 + y^2 == 0, 1,
If[x^2 + y^2  (PolarPolygonSide[sidelength, theta, sides]^2) <
0, 1, 0]] // Return];
sides = 3; (* number of sides of the polygon *)
(* The following step increments seem to be small enough for sidelengths up to 10 *)
dstep = 0.01; (* scanning step on x and y *)
phasestep = 2 Pi/3000; (* orientation angular increment step *)
seq = {};
Do[
npoints = {}; k = 0;
Do[Do[Do[
Do[Do[
k =
k + TruePointInsidePhase[i + di, j + dj, sidelength, phase,
sides]
, {i, sidelength  1, sidelength + 1, 1}], {j, sidelength 
1, sidelength + 1, 1}];
AppendTo[npoints, k];
k = 0;
, {dj, 0, 1/2, dstep}], {di, 0, 1/2, dstep}], {phase, 0, 2 Pi/3,
phasestep}] // Quiet;
temp = npoints // Min;
AppendTo[seq, temp];
, {sidelength, 0, 10, 1}]
seq
