Polygons in a cube

Question

Inspired in part by this Mathologer video on gorgeous visual "shrink" proofs, and my general interest in the topic, this challenge will have you count regular polygons with integer coordinates in 3D.

You'll be provided an input n, which is a non-negative integer. Your program should find the number of subsets of \$\{0, 1, \dots, n\}^3\$ such that the points are the vertices of a regular polygon. That is, the vertices should be 3D coordinates with nonnegative integers less than or equal to \$n\$.

Examples

For \$n = 4\$, there are \$2190\$ regular polygons: \$1264\$ equilateral triangles, \$810\$ squares, and \$116\$ regular hexagons. An example of each:

• Triangle: \$(1,0,1), (0,4,0), (4,3,1)\$
• Square: \$(1,0,0), (4,3,0), (3,4,4), (0,1,4)\$
• Hexagon: \$(1,1,0), (0,3,1), (1,4,3), (3,3,4), (4,1,3), (3,0,1)\$

The (zero-indexed) sequence begins:

0, 14, 138, 640, 2190, 6042, 13824, 28400, 53484, 94126, 156462, 248568, 380802, 564242, 813528, 1146472, 1581936, 2143878, 2857194, 3749240, 4854942, 6210442

Rules

To prevent the most naive and uninteresting kinds of brute-forcing, your program must be able to handle up to \$a(5) = 6042\$ on TIO.

This is a code-golf challenge, so the shortest code wins.

---

This is now on the On-Line Encyclopedia of Integer Sequences as A338323.

Answer 1

Python 2, 330 313 bytes

import numpy as N,itertools as I
L=list
D=N.dot
n=input()
def G((a,b,c)):
 A=a
 while 1:
	u=b-a;v=c-b;d=D(v,v)
	if D(u,u)-d:return 0
	p=2*D(u,v)*v
	a,b,c=b,c,c+a-b+p/d
	if(p%d|(c<0)|(c>n)).any()or L(a)<L(A):return 0
	if L(a)==L(A):return 1
print sum(map(G,I.permutations(N.indices((n+1,)*3).reshape(3,-1).T,3)))/2

Try it online!

TIO spent 36.893 s calculating a(5) = 6042, so this solution is just fast enough.

Given 3 consecutive points of a regular polygon a,b,c, the next point is d = c + a - b + 2*proj(b-a, c-b) (where proj is vector projection). My solution iterates through all triples of points and determines if they form a polygon using this formula.

Verified up to a(7) = 28400 on my machine.

-17 bytes from ovs