Lattice Points of the Hexagonal Lattice in a Disc by a Convex Hull Walk: the Sequences A038589, A014201 and A307109

Let N(M) be the number of integer pairs (x,y) with x^2+xy+y^2 <= M, representing the number of points of the planar hexagonal lattice A_2 inside a closed disc of squared radius M. This arithmetic function determines the OEIS sequences A038589, A014201, and A307109, which are traditionally computed using double loops or character sum enumerations requiring O(M) or O(sqrt(M)) operations. We prove a rigorous 12-fold sector decomposition identity that reduces the evaluation of N(M) to an isolated monotone sum S over a single unique arc of the underlying ellipse where the absolute slope is bounded between 1/2 and 1. We derive an exact, deterministic sublinear algorithm to compute S by executing a fast walk along the upper convex hull of the lattice points. The Stern-Brocot descent step is tracked analytically using a convex quadratic characterization of the principal norm form x^2+xy+y^2, bounding intermediate steps via an exact integer square root. The developed algorithm achieves the exact count in O(M^(1/3) log M) operations, requiring only 0.78 * M^(1/3) iterations, which was verified numerically up to M=10^20. In Python 3.12, computing N(10^20 + 12345) takes 11.2 s, while smaller scales such as M=10^16 take 0.46 s, improving upon the vectorized Dirichlet character sum baseline by several orders of magnitude. The full source code and cross-verification logs are provided.