Python 2, 43 bytes

f=lambda n,a=1:n*n<a/3or n*n/a*6-f(n,a+a%3)

Try it online!

This is black magic.

(削除) Offering 250 rep for a written-up proof. (削除ここまで) See Lynn's answer for a proof and explanation.

• 7
  \$\begingroup\$ How does this work? I've been wondering for a good 30 minutes... It looks so simple but I can't find a relationship between that recursion and circles... \$\endgroup\$

  JungHwan Min

  – JungHwan Min

  2018年04月29日 09:24:26 +00:00

  Commented Apr 29, 2018 at 9:24
• 7
  \$\begingroup\$ @JungHwanMin My proof is an epic journey through plane geometry, Eisenstein integers, factorization over number fields, quadratic reciprocity, arithmetic progressions, and interchanging summations -- all for such a simple expression. Writing it all would be a major undertaking that I don't have time for now, so I'm hoping someone else will give a proof, likely a simpler one that mine that makes the connection clearer. \$\endgroup\$

  xnor

  – xnor

  2018年05月01日 00:44:32 +00:00

  Commented May 1, 2018 at 0:44
• 17
  \$\begingroup\$ Proof. This is longer than Lynn's but more self-contained: it makes no use of unproven assertions about factorisation over the Eisenstein integers. \$\endgroup\$

  Peter Taylor

  – Peter Taylor

  2018年05月12日 17:14:44 +00:00

  Commented May 12, 2018 at 17:14
• 2
  \$\begingroup\$ @PeterTaylor Cheddar Monk? As in Darths & Droids? \$\endgroup\$

  Neil

  – Neil

  2019年02月07日 20:29:34 +00:00

  Commented Feb 7, 2019 at 20:29
• 3
  \$\begingroup\$ @Neil, congratulations on being the first person ever to ask! I registered the domain to use it as a bargaining chip for Negotiation, Level 1 in the Academy. \$\endgroup\$

  Peter Taylor

  – Peter Taylor

  2019年02月07日 20:44:39 +00:00

  Commented Feb 7, 2019 at 20:44

Haskell, 48 bytes

f n=1+6*sum[(mod(i+1)3-1)*div(n^2)i|i<-[1..n^2]]

Try it online!

Uses xnor's "black magic" formula:

$$f(n)=1+6\sum_{a=0}^\infty \left\lfloor \frac{n^2}{3a+1}\right\rfloor - \left\lfloor \frac{n^2}{3a+2}\right\rfloor$$

A proof of its correctness, and an explanation of how xnor managed to express it in 43 bytes of Python, can be found here.

Long story short: we count Eisenstein integers of norm \1ドル \le N \le n^2\$, by factoring \$N = (x+y\omega)(x+y\omega^*)\$ into Eisenstein primes and counting how many solutions for \$(x,y)\$ come out of the factorization. We recognize the number of solutions as being equal to

$6ドル \times ((\text{# of divisors of }N \equiv 1\space(\text{mod }3)) - (\text{# of divisors of }N \equiv 2\space(\text{mod }3)))$$

and apply a clever trick to make that really easy to compute for all integers between \1ドル\$ and \$n^2\$ at once. This yields the formula above. Finally, we apply some Python golf magic to end up with the really tiny solution xnor found.

• 5
  \$\begingroup\$ I certainly didn't expect this when xnor said "there's some deep mathematical insights behind golfing the problem". \$\endgroup\$

  Bubbler

  – Bubbler

  2018年05月10日 23:18:15 +00:00

  Commented May 10, 2018 at 23:18
• \$\begingroup\$ In this question I've just said "But if you find yourself considering computations like matrix division or using Eisenstein integers or Euclid's algorithm in the complex plane please save that for the follow-up code-challenge question" \$\endgroup\$

  uhoh

  – uhoh

  2020年10月20日 11:27:25 +00:00

  Commented Oct 20, 2020 at 11:27

Wolfram Language (Mathematica), (削除) 53 (削除ここまで) (削除) 51 (削除ここまで) 50 bytes

^{-1 byte thanks to @miles}

Sum[Boole[x(x+y)+y^2<=#^2],{x,-2#,2#},{y,-2#,2#}]&

Try it online!

How?

Instead of thinking in this:

enter image description here

Think of it like this:

enter image description here

So we apply the tranformation matrix [[sqrt(3)/2, 0], [1/2, 1]] to transform the second figure to the first one.

Then, we must find the circle in the triangular grid in terms of Cartesian coordinates.

(sqrt(3)/2 x)^2 + (1/2 x + y)^2 = x^2 + x y + y^2

So we find lattice points x, y such that x^2 + x y + y^2 <= r^2

For example, with r = 3:

enter image description here

• 1
  \$\begingroup\$ FYI, the formula x^2+x y+y^2 can also be derived from the Law of Cosines with 120 degrees. \$\endgroup\$

  Bubbler

  – Bubbler

  2018年04月29日 16:38:24 +00:00

  Commented Apr 29, 2018 at 16:38
• 3
  \$\begingroup\$ x^2+x y+y^2 -> x(x+y)+y^2 saves a byte \$\endgroup\$

  miles

  – miles

  2018年04月30日 07:31:59 +00:00

  Commented Apr 30, 2018 at 7:31
• \$\begingroup\$ The formula x^2 + xy + y^2 can also be derived from the norm of an Eistenstein integer, which is a^2 - ab + b^2. Note that the sign of a and b is irrelevant except in the term ab so it has the same amount of solutions. \$\endgroup\$

  orlp

  – orlp

  2018年05月01日 13:07:44 +00:00

  Commented May 1, 2018 at 13:07

Wolfram Language (Mathematica), 48 bytes

Based on OEIS A004016.

1+6Sum[DivisorSum[i,#~JacobiSymbol~3&],{i,#^2}]&

Try it online!

CJam (24 bytes)

{_*_,f{)_)3%(@@/*}1b6*)}

This is an anonymous block (function) which takes one argument on the stack and leaves the result on the stack. Online test suite. Note that the two largest cases are too slow.

Explanation

alephalpha noted in a comment on the question that

It is the sum of the first n^2+1 terms of OEIS A004016

and xnor's answer implements this sum (although I'm not sure whether their unposted proof uses it explicitly) as $$f(n) = 1 + 6 \sum_{a=0}^\infty \left\lfloor\frac{n^2}{3a+1}\right\rfloor - \left\lfloor\frac{n^2}{3a+2}\right\rfloor$$

My proof of correctness of that formula is based on some information gleaned from alephalpha's OEIS link:

G.f.: 1 + 6*Sum_{n>=1} x^(3*n-2)/(1-x^(3*n-2)) - x^(3*n-1)/(1-x^(3*n-1)). - Paul D. Hanna, Jul 03 2011

for which the relevant reference is the paper by Hirschhorn. An elementary proof is possible using nothing more than a basic understanding of complex numbers (cube roots of unity, magnitude), the concept of generating functions, the derivative of \$x^a\,ドル and the chain rule of differentiation. In summary, we first prove from first principles the Jacobi triple-product identity $$\prod_{k=0}^\infty (1-q^{k+1})(1 + xq^{k+1})(1 + x^{-1}q^k) = \sum_{k\in \mathbb{Z}} q^{k(k+1)/2}x^k$$ That then bootstraps a proof that $$\sum_{m,n \in \mathbb{Z}} \omega^{m-n} q^{m^2+mn+n^2} = \prod_{k=1}^\infty \frac{(1-q^k)^3}{1-q^{3k}}$$ where \$\omega\$ is a primitive cube root of unity. The final big step is to use this to show that $$\sum_{m,n \in \mathbb{Z}} q^{m^2+mn+n^2} = 1 + 6 \sum_{k \ge 0} \left(\frac{q^{3k+1}}{1-q^{3k+1}} - \frac{q^{3k+2}}{1-q^{3k+2}} \right)$$

Code dissection

{ e# Define a block. Stack: ... r
 _* e# Square it
 _,f{ e# Map with parameter: invokes block for (r^2, 0), (r^2, 1), ... (r^2, r^2-1)
 ) e# Increment second parameter. Stack: ... r^2 x with 1 <= x <= r^2
 _)3%( e# Duplicate x and map to whichever of 0, 1, -1 is equal to it (mod 3)
 @@/* e# Evaluate (r^2 / x) * (x mod 3)
 }
 1b6* e# Sum and multiply by 6
 ) e# Increment to count the point at the origin
}

APL (Dyalog Classic), 23 bytes

×ばつ+/-/⌊⍵÷1+3⊥ ×ばつ⍨

Try it online!

tribute to xnor's and lynn's answers

the last test is commented because it needs more memory, e.g. MAXWS=200M in the env

J, 27 bytes

[:+/@,*:>:(*++&*:)"{~@i:@+:

Try it online!

Based on JungHwan Min's method.

Explanation

[:+/@,*:>:(*++&*:)"{~@i:@+: Input: n
 +: Double
 i: Range [-2n .. 2n]
 "{~ For each pair (x, y)
 *: Square both x and y
 + Add x^2 and y^2
 + Plus
 * Product of x and y
 >: Less than or equal to
 *: Square of n
 , Flatten
 +/ Reduce by addition

Jelly, 14 bytes

ḤŒR+2_×ばつʋþ`F1⁄2»ċ

Uses @JungHwanMin's method.

Try it online!

Jelly, (削除) 15 (削除ここまで) 13 bytes

-2 thanks to Dennis (just increment the square to avoid concatenation of a zero; avoid head by using a post-difference modulo-slice rather than a pre-difference slice)

^{Uses the "black magic" method of honing in on the answer that was exposed by xnor in their Python answer, but uses iteration rather than recursion (and a little less calculation)}

×ばつ6C

A monadic link accepting a non-negative integer and returning a positive integer.

Try it online! Or see the test-suite.

How?

×ばつ6C - Main Link: non-negative integer, n e.g. 7
2 - square 49
 $ - last two links as a monad:
 ‘ - increment 50
 Ѐ - map across (implicit range of) right with:
 : - integer division [49,24,16,12,9,8,7,6,5,4,4,4,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0]
 I - incremental differences [-25,-8,-4,-3,-1,-1,-1,-1,-1,0,0,-1,0,0,0,-1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1]
 m3 - every third element [-25,-3,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,-1]
 S - sum (vectorises) -31
 ×ばつ6 - multiply by six -186
 C - complement (1-X) 187

05AB1E, (削除) 15 (削除ここまで) 14 bytes

nD>L÷3円ιнO6*±Ì

Port of @JonathanAllans Jelly answer, which in turn is a derivative from @xnor's 'black magic' formula.

Try it online or verify all test cases.

Explanation:

n # Square the (implicit) input-integer
 D> # Duplicate it, and increase the copy by 1
 L # Create a list in the range [1, input2+1]
 ÷ # Integer divide input2 by each of these integers
 \ # Take the deltas / forward differences
 3ι # Uninterleave this list into 3 parts
 # i.e. [a,b,c,d,e] → [[a,d],[b,e],[c]]
 н # Pop and leave the first inner list
 O # Take the sum of this list
 6* # Multiply it by 6
 ± # Take the bitwise NOT (-n-1)
 Ì # And increase it by 2
 # (after which the result is output implicitly)

JavaScript (ES6), 65 bytes

This is a port of @JungHwanMin's solution.

f=(n,y=x=w=n*2)=>y-~w&&(x*x+x*y+y*y<=n*n)+f(n,y-=--x<-w&&(x=w,1))

Try it online!

Original answer (ES7), 70 bytes

Simply walks through the grid and counts the matching points.

f=(n,y=x=n*=2)=>y+n+2&&(x*x*3+(y-x%2)**2<=n*n)+f(n,y-=--x<-n&&(x=n,2))

Try it online!

• 1
  \$\begingroup\$ Porting xnor's answer is shorter: 42 bytes (outputs true instead of 1; 46 if we also integer-divide it). And I don't know JavaScript well enough to golf the integer-divisions ~~(a/b), but I'm sure there is a shorter way for those as well.. \$\endgroup\$

  Kevin Cruijssen

  – Kevin Cruijssen

  2019年08月30日 12:09:44 +00:00

  Commented Aug 30, 2019 at 12:09

Wolfram Language (Mathematica), 39 bytes

Length@Solve[x(x+y)+y^2<=#^2,Integers]&

Try it online!

Using JungHwan Min's coordinate transformation and simply counting the solutions over the integers.

Java 8, 65 bytes

n->f(n,1)int f(int n,int a){return n*n<a/3?1:n*n/a*6-f(n,a+a%3);}

Port of @xnor's Python 2 answer.

Try it online.

Pari/GP, 42 bytes

Using the built-in qfrep.

n->1+2*vecsum(Vec(qfrep([2,1;1,2],n^2,1)))

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.

Try it online!

C# (Visual C# Interactive Compiler), 68 bytes

n=>{int g(int x,int y)=>x*x<y/3?1:x*x/y*6-g(x,y+y%3);return g(n,1);}

Try it online!

Same as everyone else, unfortunately. I know there's probably a better way of writing this, but declaring and calling a lambda at the same time in c# is not exactly something I do, well, ever. Though in my defense, I can't think of a good reason (outside code golf, of course) to do so. Still, if someone knows how you can do this, let me know and/or steal the credit, I guess.

• code-golf
• math
• geometry
• hexagonal-grid