[Question]
  [
# Background
A **[triangular grid](http://mathworld.wolfram.com/TriangularGrid.html)** is a grid formed by tiling the plane regularly with equilateral triangles of side length 1. The picture below is an example of a triangular grid.
![](https://i.stack.imgur.com/pmooq.gif)
A **triangular lattice point** is a vertex of a triangle forming the triangular grid.
The **origin** is a fixed point on the plane, which is one of the triangular lattice points.
# Challenge
Given a non-negative integer `n`, find the number of triangular lattice points whose Euclidean distance from the origin is less than or equal to `n`.
# Example
The following figure is an example for `n = 7` (showing only 60-degree area for convenience, with point A being the origin):
![](https://i.stack.imgur.com/xw9yl.png)
# Test Cases
```
Input | Output
---------------
    0 |       1
    1 |       7
    2 |      19
    3 |      37
    4 |      61
    5 |      91
    6 |     127
    7 |     187
    8 |     241
    9 |     301
   10 |     367
   11 |     439
   12 |     517
   13 |     613
   14 |     721
   15 |     823
   16 |     931
   17 |    1045
   18 |    1165
   19 |    1303
   20 |    1459
   40 |    5815
   60 |   13057
   80 |   23233
  100 |   36295
  200 |  145051
  500 |  906901
 1000 | 3627559
```
**Hint**: This sequence is *not* [OEIS A003215](https://oeis.org/A003215).
# Rules
Standard rules for [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'") apply. The shortest submission wins.
Please include how you solved the challenge in your submission.
  
[Answer]
# [Python 2](https://docs.python.org/2/), 43 bytes
```
f=lambda n,a=1:n*n<a/3or n*n/a*6-f(n,a+a%3)
```
[Try it online!](https://tio.run/##FYzBCsIwEETv@Yq9CEmt1rbQQ7Efs2qiC80mbOKhXx/XgYGBx7x81E/iqbWw7RgfLwTucRtX7viOw5wEdA3YLZdglZzxNLtGMSepUI5itNfiq/jnVwol3ilStdNN44wJ/z8QgyC/vV1Gt2YhrqAu9TnXfg "Python 2 – Try It Online")
This is black magic. 
**[Offering 250 rep](https://codegolf.meta.stackexchange.com/a/16359/20260) for a written-up proof.** See [Lynn's answer](https://codegolf.stackexchange.com/a/164474/20260) for a proof and explanation. 
[Answer]
# [Haskell](https://www.haskell.org/), 48 bytes
```
f n=1+6*sum[(mod(i+1)3-1)*div(n^2)i|i<-[1..n^2]]
```
[Try it online!](https://tio.run/##y0gszk7Nyfn/P00hz9ZQ20yruDQ3WiM3P0UjU9tQ01jXUFMrJbNMIy/OSDOzJtNGN9pQTw/IiY39n5uYmadgq1BQlJlXoqCRpmBoYKD5HwA "Haskell – 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**](http://foldr.moe/blackmagic.pdf).
> 
> Long story short: we count [Eisenstein integers](https://en.wikipedia.org/wiki/Eisenstein_integer) of norm \$1 \le N \le n^2\$, by factoring \$N = (x+y\omega)(x+y\omega^\*)\$ into [Eisenstein primes](https://en.wikipedia.org/wiki/Eisenstein_prime) 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.
> 
> 
> 
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), ~~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!](https://tio.run/##y00syUjNTSzJTE78n2b7P7g0N9opPz8nNbpCo0K7UlO7Ms7IxlY5zihWp7pCR9dIWcdIuVanuhLGjFX775IfHVCUmVcSnaega6eQFp0XG6ujUJ2no2Cgo2BkUBv7HwA "Wolfram Language (Mathematica) – Try It Online")
## How?
Instead of thinking in this:
[![enter image description here](https://i.stack.imgur.com/C13c8m.png)](https://i.stack.imgur.com/C13c8m.png)
Think of it like this:
[![enter image description here](https://i.stack.imgur.com/U8KTVm.png)](https://i.stack.imgur.com/U8KTVm.png)
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](https://i.stack.imgur.com/1Spuu.png)](https://i.stack.imgur.com/1Spuu.png)
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 48 bytes
Based on OEIS [A004016](https://oeis.org/A004016).
```
1+6Sum[DivisorSum[i,#~JacobiSymbol~3&],{i,#^2}]&
```
[Try it online!](https://tio.run/##y00syUjNTSzJTE78n6Zg@99Q2yy4NDfaJbMsszi/CMTM1FGu80pMzk/KDK7MTcrPqTNWi9WpBorGGdXGqv13yY8OKMrMK4nOU9C1U0iLzouN1VGoztNRMNBRMDKojf0PAA "Wolfram Language (Mathematica) – Try It Online")
[Answer]
### 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](http://cjam.aditsu.net/#code=q'%7C-~%5D2%2F%7B%0A~%5C%0A%0A%7B_*_%2Cf%7B)_)3%25(%40%40%2F*%7D1b6*)%7D%0A%0A~%3Do%0A%7D%2F&input=%20%20%20%200%20%7C%20%20%20%20%20%20%201%0A%20%20%20%201%20%7C%20%20%20%20%20%20%207%0A%20%20%20%202%20%7C%20%20%20%20%20%2019%0A%20%20%20%203%20%7C%20%20%20%20%20%2037%0A%20%20%20%204%20%7C%20%20%20%20%20%2061%0A%20%20%20%205%20%7C%20%20%20%20%20%2091%0A%20%20%20%206%20%7C%20%20%20%20%20127%0A%20%20%20%207%20%7C%20%20%20%20%20187%0A%20%20%20%208%20%7C%20%20%20%20%20241%0A%20%20%20%209%20%7C%20%20%20%20%20301%0A%20%20%2010%20%7C%20%20%20%20%20367%0A%20%20%2011%20%7C%20%20%20%20%20439%0A%20%20%2012%20%7C%20%20%20%20%20517%0A%20%20%2013%20%7C%20%20%20%20%20613%0A%20%20%2014%20%7C%20%20%20%20%20721%0A%20%20%2015%20%7C%20%20%20%20%20823%0A%20%20%2016%20%7C%20%20%20%20%20931%0A%20%20%2017%20%7C%20%20%20%201045%0A%20%20%2018%20%7C%20%20%20%201165%0A%20%20%2019%20%7C%20%20%20%201303%0A%20%20%2020%20%7C%20%20%20%201459%0A%20%20%2040%20%7C%20%20%20%205815%0A%20%20%2060%20%7C%20%20%2013057%0A%20%20%2080%20%7C%20%20%2023233%0A%20%20100%20%7C%20%20%2036295%0A%20%20200%20%7C%20%20145051). Note that the two largest cases are too slow.
### Explanation
[alephalpha](/users/9288/alephalpha) noted in a comment on the question that
> 
> It is the sum of the first n^2+1 terms of [OEIS A004016](//oeis.org/A004016)
> 
> 
> 
and [xnor's answer](/a/163708/194) 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](http://cheddarmonk.org/papers/triangle-lattice.pdf) 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
}
```
[Answer]
# [APL (Dyalog Classic)](https://www.dyalog.com/), 23 bytes
```
{1+6×+/-/⌊⍵÷1+3⊥¨⍳⍵2}×⍨
```
[Try it online!](https://tio.run/##PY@xSsRAEIZ7nyJdihBuZ2d3k/VtwkkkGIjcXSNynQieRmzkKgsL4XoRwdK8yb5InNnsuEXYb/5/li/NdV9e3DT9cFmu@2a77dZzeH7thnD/ouaWvrdQuOlYrMpVeDqE8Wv6hgLD4eP3FMZPYr2fjmE8zfOOyzQID@8burZh/DnPh6s8PN7lbdP1OQ0o3uzPsixT2S5bDjDCP1aMWhA8IwpiTI2gi7tW0Ed0CUHHciVYR6wTahPLPiEqRhArdFwGsTLIGiBWFmIqVg6QUawqHZ8Sq1rHVKw8xjRZgTKWMVkBuIjJClDxrk5WYCxrmIS2Bi67Balr2apeUKNG2gW1IDrtqawXpIeUJQ27oFfO0@@H8Y3qPKF2Za3/Aw "APL (Dyalog Classic) – Try It Online")
tribute to [xnor's](https://codegolf.stackexchange.com/questions/163704/triangular-lattice-points-close-to-the-origin/163708#163708) and [lynn's](https://codegolf.stackexchange.com/questions/163704/triangular-lattice-points-close-to-the-origin/164474#164474) answers
the last test is commented because it needs more memory, e.g. `MAXWS=200M` in the env
[Answer]
# [J](http://jsoftware.com/), 27 bytes
```
[:+/@,*:>:(*++&*:)"{~@i:@+:
```
[Try it online!](https://tio.run/##y/qfVqxga6VgoADE/6OttPUddLSs7Kw0tLS11bSsNJWq6xwyrRy0rf5rcinpKaingRSrK@go1FoppBVzcaUmZ@QraOjopSkZaCpk6ikYGf4HAA "J – Try It Online")
Based on JungHwan Min's [method](https://codegolf.stackexchange.com/questions/163704/triangular-lattice-points-close-to-the-origin/163711#163711).
## 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
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 14 bytes
```
ḤŒR+²_×ʋþ`F½»ċ
```
Uses [@JungHwanMin's method](https://codegolf.stackexchange.com/a/163711/12012).
[Try it online!](https://tio.run/##ATwAw/9qZWxsef//4bikxZJSK8KyX8OXyovDvmBGwr3Cu8SL/zByMjA7NDAsNjAsODAsMTAwwrXFvMOH4oKsR/8 "Jelly – Try It Online")
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~15~~ 13 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
-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](https://codegolf.stackexchange.com/a/163708/53748), but uses iteration rather than recursion (and a little less calculation)
```
²:Ð€‘$Im3S×6C
```
A monadic link accepting a non-negative integer and returning a positive integer.
**[Try it online!](https://tio.run/##ASEA3v9qZWxsef//wrI6w5DigqzigJgkSW0zU8OXNkP///8xMDA "Jelly – Try It Online")** Or see the [test-suite](https://tio.run/##y0rNyan8///QJqvDEx41rXnUMEPFM9c4@PB0M@f/R/ccbgeKuQGx@///OZnFJRpFiXnpqRpGhpqa2tEmBjpmBjoWBjqGBgY6RkBsagBmG8QCAA "Jelly – Try It Online").
### How?
```
²:Ð€‘$Im3S×6C - Main Link: non-negative integer, n     e.g. 7
²             - 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
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E/wiki/Commands), ~~15~~ 14 [bytes](https://github.com/Adriandmen/05AB1E/wiki/Codepage)
```
nD>L÷¥3ιнO6*±Ì
```
Port of [*@JonathanAllan*s Jelly answer](https://codegolf.stackexchange.com/a/163729/52210), which in turn is a derivative from [@xnor's 'black magic' formula](https://codegolf.stackexchange.com/a/163708/52210).
[Try it online](https://tio.run/##AR8A4P9vc2FiaWX//25EPkzDt8KlM8650L1PNirCscOM//83) or [verify all test cases](https://tio.run/##yy9OTMpM/V9Waa@k8KhtkoKSfeX/PBc7n8PbDy01Prfzwl5/M61DGw/3/Nf5H22gY6hjpGOsY6JjqmOmY65joWOpYwgUNNQxNNIxNNYxNNExNNUxNNMxNNcxtNAxtNQxMtAxMdAxM9CxACozMADyDXRMDcBsg1gA).
**Explanation:**
```
n               # Square the (implicit) input-integer
 D>             # Duplicate it, and increase the copy by 1
   L            # Create a list in the range [1, input²+1]
    ÷           # Integer divide input² 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)
```
[Answer]
# JavaScript (ES6), 65 bytes
This is a port of [@JungHwanMin's solution](https://codegolf.stackexchange.com/a/163711/58563).
```
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!](https://tio.run/##HctBDsIwDETRq3hVJXGNoFti7lKVBoEqG7WocTZcPaTsRpr3X@M@btP6fH9I9D7XmthJX9g4s4TB863QN3eds2BooWAJJbZHPKYDEhNZpL/g3F@8r0lXJ8BwvoJAZBiOgehhUtl0mU@LPloKrW/6Bw "JavaScript (Node.js) – 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!](https://tio.run/##FcsxDsIwDEDRq3ihiuO6QmYk5i5V2yBQ5aAWoeT0Idn@8P57/s3ncrw@X7a0brVGdTYWzWpeBfVRyEiGwWWf/Y1c4XwR9F5CA4YUu2ZlzoGtM7VREGtMhzNQuN7BIChIDyKEJdmZ9m3a07Ot0P6m/w "JavaScript (Node.js) – Try It Online")
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 39 bytes
```
Length@Solve[x(x+y)+y^2<=#^2,Integers]&
```
[Try it online!](https://tio.run/##FctNC4JAFIXhv3JBkMIbeSczgwy3QYug5aAgMX5ATmBDKNFvn46Ls3qfM9SuM0Pt@kftG8r91djWdcX99fwYPa2maF5Hc6VOeVApvlhnWjO@y9Dfxt46HWzOjQ7KMqRtQd@YSZgU044pYdozpUwHpozpiLRkdAEQCAERGAESKAETOAWXYCmWLZ84/vk/ "Wolfram Language (Mathematica) – Try It Online")
Using [JungHwan Min's coordinate transformation](https://codegolf.stackexchange.com/a/163711/87058) and simply counting the solutions over the integers.
[Answer]
# 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](https://codegolf.stackexchange.com/a/163708/52210).
[Try it online.](https://tio.run/##LU9NT8MwDL3zK54mTUqWtWNMcFhX@AXssiNwMF2KsrXulKZDaOpvL05Asv387ecTXSk7Hc9T1VDf45Uc3@4Ax8H6miqLfQxTApWKlnUhmVFUZA9GiYmz51rxcq2L2FH/9S2jJX3zNgyewQve0Wrzst6Kt6LFUxZHyNB8o4txAi7DZ@Mq9IGCwLVzR7RCRx2Cd/z19hF3JS6RJFq5y/Y7BSpRAurOp9OuvC/gduXDo4AxOhWBw08fbJt3Q8gvsjI0rJyZbd/DzLQ55/Ke/v9tnH4B)
[Answer]
# [Pari/GP](http://pari.math.u-bordeaux.fr/), 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!](https://tio.run/##K0gsytRNL/ifpmD7P0/XzlDbSKssNbm4NFcjLDVZozCtKLVAI9pIx9DaUMcoVicvDsjU1NT8n5ZfpJGnYKtgoKNgBMQFRZl5JUABJQVdOyCRppEHUgQA "Pari/GP – Try It Online")
[Answer]
# [C# (Visual C# Interactive Compiler)](http://www.mono-project.com/docs/about-mono/releases/5.0.0/#csc), 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!](https://tio.run/##FYtBC8IgHMXvfYpdAl1ucww65LRDbKedIui8xC2h/oI6UKLPbnZ4jx8/3pOukk6ncQPZa/AkRyw8ARefjMWK/h3@uoiYi1CGPjbduT1lamJ5rFYUSDzEfYeZVX6zkD9AWsy@ie0G2N7Kzo@Xqq8zrApR0lJKcX0zk3Ye4Xo0dpjlEwUuLgacycu71V5NGhRaUMAYs/QD "C# (Visual C# Interactive Compiler) – 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. 
]