# Polygons in a cube

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:

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 challenge, so the shortest code wins.

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

• So the gist of the video is that we'll ever have only regular 3/4/6-gons made of 3D grid points, however large the grid is, right? Oct 23, 2020 at 0:34
• @Bubbler—that's right! Oct 23, 2020 at 0:35
• I conjecture that the squares only appear on the planes parallel to the cube's faces, and the triangles and hexagons only appear on the (1,1,1)-lattice equivalents. Is this true? Oct 23, 2020 at 0:57
• @Bubbler, I've edited the example of the square to provide a square that is not parallel to a face. Oct 23, 2020 at 1:06
• @Bubbler, for a triangle, see $(0,4,0),(1,0,1),(4,3,1)$. Oct 23, 2020 at 1:11

# 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

• You can get the coordinates a little shorter using pure numpy: Try it online!. See this StackOverflow answer for details.
– ovs
Oct 26, 2020 at 9:46
• And (p%d|(c<0)|(c>n)).any() saves a few more bytes
– ovs
Oct 26, 2020 at 9:53