Counting Eisenstein primes

Julia, 66 62 60 bytes

!n=sum(isprime,[a<1<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² − ab + 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² − ab + 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<1<b%3.

(For the remaining green lines, we will never return a false positive: if a = b or b = 0 then a² − ab + b² = 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<1<b%3 – then the expression reduces to b. This doesn’t seem to be shorter, but it’s neat!