Counting Eisenstein primes Introduction

Eisenstein integers are complex numbers of the form

a+bω

Where a,b are integers, and 

ω = e^(2πi/3)

The Eisenstein integers form a triangular lattice in the complex plane:



We say that an Eisenstein integer z=a+bω is prime if it cannot be written as the product of two non-unit (not 1,-1,ω,-ω,ω^2, or -ω^2) Eisenstein integers

Program

Input: A natural number n.

Output: The number of Eisenstein primes which are of the form a+bω for which a,b are natural numbers (including zero) less than or equal to n

Test Cases

0 → 0

1 → 0

2 → 5

3 → 9

4 → 13

5 → 20

Scoring

This is code-golf, so least amount of bytes wins
Julia, 66 62 60 bytes

!n=sum(isprime,[a&lt;1&lt;b%3?b:a^2-a*b+b^2for a=[0;0;0:n],b=0:n])


Try it online!

Explanation

We’re interested in the primes in this parallelogram on the complex plane (example for n = 4):



We can split them into primes on the green lines, and on the gray lines.

Wikipedia tells me an Eisenstein number z is a green line Eisenstein prime iff |z| is a natural prime equal to 2 mod 3.

It also says z is a gray line Eisenstein prime iff |z|² = a²&nbsp;&minus;&nbsp;ab&nbsp;+&nbsp;b² is a natural prime.



So we loop over a = 0…n and b = 0…n, and check:


If (a = 0 or b = 0 or a = b) and max(a, b) % 3 = 2, then count whether max(a, b) is a prime.
Else, count whether a²&nbsp;&minus;&nbsp;ab&nbsp;+&nbsp;b² is a prime.


However, we can abuse the distribution’s symmetry. Instead of counting each green line once, we can just count one green line thrice! That is, only check a = 0 and increment the counter by three when we find a green line prime. The a=[0;0;0:n] achieves exactly this.

Since we know we’re only considering the green line a = 0, we can replace max(a, b) by b.

The “green line condition” is nicely expressed in Julia using operator chaining: a&lt;1&lt;b%3.

(For the remaining green lines, we will never return a false positive: if a = b or b = 0 then a²&nbsp;&minus;&nbsp;ab&nbsp;+&nbsp;b²&nbsp;=&nbsp;a², which can’t be prime.)

Ideas

Maybe, instead of writing a^2-a*b+b^2, I can conditionally replace the exponent at b by 1 when a&lt;1&lt;b%3 – then the expression reduces to b. This doesn’t seem to be shorter, but it’s neat!