> with(numtheory):
# This procedure is shown to give the correct number of equilateral
# triangles whose coordinates are integers 
# between 0 and n as long as the Diophantine equation 3d^2=a^2+b^2+c^2
# for d an odd number between 1 and n
# does not have a degenerate solution, (i.e. gcd(a,b,c)=1 but
# min[gcd(a,d),gcd(b,d),gcd(c,d)]>1).  
# The existence of degenerate solutions is known but the smallest one
# seems to be far away (d>400).
# The number of such equilateral triangles for n=400 can be estimates
# to about (n+1)^5.5=207631886300000.
# We would like to thank Ray J. Chandler who helped correcting an
# earlier version of this code.
# Generating the sides less than n^2. 
> sides:=proc(n)
> local i,j,k,L,a,m,p,q,r,max;
> L:={};max:=n^2;
> for i from 2 to max do
> a:=ifactors(i);
> k:=nops(a[2]);r:=0;
>     for j from 1 to k do 
>       m:=a[2][j][1]; p:=m mod 3; q:=a[2][j][2] mod 2;
>       if  r=0 and p=2 and q=1 then r:=1 fi;
>     od;
> if r=0 then L:=L union {i};fi;
> od;
> L:=L union {1};
> L:=convert(L,list);
> end:
# Adding a vector to the coordinates of a triangle T.
> addvect:=proc(T,v)
> local Q,a,b,c;
> Q:=convert(T,list);a:=v[1];b:=v[2];c:=v[3];
> {[Q[1][1]+a,Q[1][2]+b,Q[1][3]+c],[Q[2][1]+a,Q[2][2]+b,Q[2][3]+c],[Q[3]
> [1]+a,Q[3][2]+b,Q[3][3]+c]};
> end:

# Finding the orbit length under the cube group action 
> orbit1:=proc(T)
> local
> i,k,T1,a,b,c,x,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17
> ,T18,
> T19,T20,T21,T22,T23,T24,S,Q,d,a1,b1,c1;
> Q:=convert(T,list);
> a:=max(Q[1][1],Q[2][1],Q[3][1]);a1:=min(Q[1][1],Q[2][1],Q[3][1]);
> b:=max(Q[1][2],Q[2][2],Q[3][2]);b1:=min(Q[1][2],Q[2][2],Q[3][2]);
> c:=max(Q[1][3],Q[2][3],Q[3][3]);c1:=min(Q[1][3],Q[2][3],Q[3][3]);
> d:=max(a,b,c);
> T1:=T;
> T2:={[Q[1][2],Q[1][3],Q[1][1]],[Q[2][2],Q[2][3],Q[2][1]],[Q[3][2],Q[3]
> [3],Q[3][1]]};
> T3:={[Q[1][1],Q[1][3],Q[1][2]],[Q[2][1],Q[2][3],Q[2][2]],[Q[3][1],Q[3]
> [3],Q[3][2]]};

> T4:={[Q[1][1],Q[1][2],d-Q[1][3]],[Q[2][1],Q[2][2],d-Q[2][3]],[Q[3][1],
> Q[3][2],d-Q[3][3]]};
> T5:={[Q[1][2],Q[1][3],d-Q[1][1]],[Q[2][2],Q[2][3],d-Q[2][1]],[Q[3][2],
> Q[3][3],d-Q[3][1]]};
> T6:={[Q[1][1],Q[1][3],d-Q[1][2]],[Q[2][1],Q[2][3],d-Q[2][2]],[Q[3][1],
> Q[3][3],d-Q[3][2]]};

> T7:={[Q[1][1],d-Q[1][2],Q[1][3]],[Q[2][1],d-Q[2][2],Q[2][3]],[Q[3][1],
> d-Q[3][2],Q[3][3]]};
> T8:={[Q[1][2],d-Q[1][3],Q[1][1]],[Q[2][2],d-Q[2][3],Q[2][1]],[Q[3][2],
> d-Q[3][3],Q[3][1]]};
> T9:={[Q[1][1],d-Q[1][3],Q[1][2]],[Q[2][1],d-Q[2][3],Q[2][2]],[Q[3][1],
> d-Q[3][3],Q[3][2]]};

> T10:={[d-Q[1][1],Q[1][2],Q[1][3]],[d-Q[2][1],Q[2][2],Q[2][3]],[d-Q[3][
> 1],Q[3][2],Q[3][3]]};
> T11:={[d-Q[1][2],Q[1][3],Q[1][1]],[d-Q[2][2],Q[2][3],Q[2][1]],[d-Q[3][
> 2],Q[3][3],Q[3][1]]};
> T12:={[d-Q[1][1],Q[1][3],Q[1][2]],[d-Q[2][1],Q[2][3],Q[2][2]],[d-Q[3][
> 1],Q[3][3],Q[3][2]]};

> T13:={[Q[1][1],d-Q[1][2],d-Q[1][3]],[Q[2][1],d-Q[2][2],d-Q[2][3]],[Q[3
> ][1],d-Q[3][2],d-Q[3][3]]};
> T14:={[Q[1][2],d-Q[1][3],d-Q[1][1]],[Q[2][2],d-Q[2][3],d-Q[2][1]],[Q[3
> ][2],d-Q[3][3],d-Q[3][1]]};
> T15:={[Q[1][1],d-Q[1][3],d-Q[1][2]],[Q[2][1],d-Q[2][3],d-Q[2][2]],[Q[3
> ][1],d-Q[3][3],d-Q[3][2]]};

> T16:={[d-Q[1][1],d-Q[1][2],Q[1][3]],[d-Q[2][1],d-Q[2][2],Q[2][3]],[d-Q
> [3][1],d-Q[3][2],Q[3][3]]};
> T17:={[d-Q[1][2],d-Q[1][3],Q[1][1]],[d-Q[2][2],d-Q[2][3],Q[2][1]],[d-Q
> [3][2],d-Q[3][3],Q[3][1]]};
> T18:={[d-Q[1][1],d-Q[1][3],Q[1][2]],[d-Q[2][1],d-Q[2][3],Q[2][2]],[d-Q
> [3][1],d-Q[3][3],Q[3][2]]};

> T19:={[d-Q[1][1],Q[1][2],d-Q[1][3]],[d-Q[2][1],Q[2][2],d-Q[2][3]],[d-Q
> [3][1],Q[3][2],d-Q[3][3]]};
> T20:={[d-Q[1][2],Q[1][3],d-Q[1][1]],[d-Q[2][2],Q[2][3],d-Q[2][1]],[d-Q
> [3][2],Q[3][3],d-Q[3][1]]};
> T21:={[d-Q[1][1],Q[1][3],d-Q[1][2]],[d-Q[2][1],Q[2][3],d-Q[2][2]],[d-Q
> [3][1],Q[3][3],d-Q[3][2]]};

> T22:={[d-Q[1][1],d-Q[1][2],d-Q[1][3]],[d-Q[2][1],d-Q[2][2],d-Q[2][3]],
> [d-Q[3][1],d-Q[3][2],d-Q[3][3]]};
> T23:={[d-Q[1][2],d-Q[1][3],d-Q[1][1]],[d-Q[2][2],d-Q[2][3],d-Q[2][1]],
> [d-Q[3][2],d-Q[3][3],d-Q[3][1]]};
> T24:={[d-Q[1][1],d-Q[1][3],d-Q[1][2]],[d-Q[2][1],d-Q[2][3],d-Q[2][2]],
> [d-Q[3][1],d-Q[3][3],d-Q[3][2]]};
> S:={T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17,T18,T19
> ,T20,T21,T22,T23,T24};
> S;
> end:
> orbit:=proc(T)
> local S,Q,T1;
> Q:=convert(T,list);
> T1:={[Q[1][3],Q[1][2],Q[1][1]],[Q[2][3],Q[2][2],Q[2][1]],[Q[3][3],Q[3]
> [2],Q[3][1]]};
> S:=orbit1(T) union orbit1(T1);
> S;
> end:
> transl:=proc(T)
> local S,Q,i,j,k,a2,b2,c2,a,b,c,d;
> Q:=convert(T,list);
> a:=max(Q[1][1],Q[2][1],Q[3][1]);
> b:=max(Q[1][2],Q[2][2],Q[3][2]);
> c:=max(Q[1][3],Q[2][3],Q[3][3]);
> d:=max(a,b,c);
> a2:=d-a;b2:=d-b;c2:=d-c;
> S:=orbit(T);
>  for i from 0 to a2 do 
>    for j from 0 to b2 do
>     for k from 0 to c2 do
>         S:=S union orbit(addvect(T,[i,j,k]));
>       od;
>    od;
>  od;
> S;
> end:
> inters:=proc(T)
> local a,b,c,Q,d,S,m,i,S1,S2;
> Q:=convert(T,list);
> a:=max(Q[1][1],Q[2][1],Q[3][1]);
> b:=max(Q[1][2],Q[2][2],Q[3][2]);
> c:=max(Q[1][3],Q[2][3],Q[3][3]);
> d:=max(a,b,c);S2:=transl(T);S:=convert(S2,list);m:=nops(S);S1:={};
> for i from 1 to m do 
> S1:=S1 union {addvect(S[i],[0,0,1])}; 
> od;
> S2 intersect S1;
> end:
# Finding the intersection along translations along [0,0,1] and [0,1,1]
> intersch:=proc(T)
> local a,b,c,Q,d,S,m,i,S1,S2,S3,S4;
> Q:=convert(T,list);
> S2:=transl(T);S:=convert(S2,list);m:=nops(S);S1:={};
> for i from 1 to m do 
> S1:=S1 union {addvect(S[i],[0,0,1])}; 
> od;
> S3:={};
> for i from 1 to m do 
> S3:=S3 union {addvect(S[i],[0,1,0])}; 
> od;
> nops(S1 intersect S3);
> end:
# Calculating the number of triangles fitting in {0,1,2,...,n}^3.
> f:=(n,d,x,y,w)->(n-d+1)^3*x-(12+3*(n-d-1)^3+24*(n-d-1)+15*(n-d-1)^2)*y
> +3*w*(n-d+1)*(n-d)^2:
>  notrincn:=proc(T,n)
> local Q,a,b,c,x,a2,b2,c2,d,y,z,w;
> Q:=convert(T,list);
> a2:=max(Q[1][1],Q[2][1],Q[3][1]);
> b2:=max(Q[1][2],Q[2][2],Q[3][2]);
> c2:=max(Q[1][3],Q[2][3],Q[3][3]);
> d:=max(a2,b2,c2);
> x:=nops(transl(T));y:=nops(inters(T));w:=intersch(T);
> if n>d then 
> z:=f(n,d,x,y,w);else z:=x;
> fi;
> z;
> end:

# Checking the side lengths
> checkeq:=proc(T)
> local d1,d2,d3,Q,D;
> Q:=convert(T,list);
> d1:=sqrt((Q[1][1]-Q[2][1])^2+(Q[1][2]-Q[2][2])^2+(Q[1][3]-Q[2][3])^2);
> d2:=sqrt((Q[1][1]-Q[3][1])^2+(Q[1][2]-Q[3][2])^2+(Q[1][3]-Q[3][3])^2);
> d3:=sqrt((Q[3][1]-Q[2][1])^2+(Q[3][2]-Q[2][2])^2+(Q[3][3]-Q[2][3])^2);
> D:={d1,d2,d3};
> end:
# Finding the solutions of a^2+b^2+c^2=3*d^2 given d.
> abcsol:=proc(d)
> local i,j,k,m,u,x,y,sol,cd;sol:={};
> for i from 1 to d do
>  u:=[isolve(3*d^2-i^2=x^2+y^2)];k:=nops(u);
>    for j from 1 to k do 
>      if rhs(u[j][1])>=i and rhs(u[j][2])>=i then
>         cd:=gcd(gcd(i,rhs(u[j][1])),rhs(u[j][2]));
>       if cd=1 then sol:=sol union
> {sort([i,rhs(u[j][1]),rhs(u[j][2])])};fi;
>      fi;
>    od;
>  od;
> convert(sol,list);
> end:
# Find the divisors d that work for a side s
> dkl:=proc(s)
> local i,x,noft,div,y,y1,z;
> x:=convert(divisors(s),list); noft:=nops(x); div:={};
> for i from 1 to noft do
> y:=s/x[i]^2; y1:=floor(y);z:=x[i] mod 2;
> if z=1 and y=y1 then div:=div union {x[i]}; fi;
> od;
> convert(div,list);
> end:
# Finding parametrizations of a given solution [a,b,c].
> findpar:=proc(a,b,c,m,n)
> local i,j,sol,mx,nx,r,s,my,ny,q,d,u,v,w,x,y,z,mu,nu,mv,nv,ef,ns,mz,nz
> ,mw,nw,om,l,uu,t;q:=a^2+b^2;
> sol:=convert({isolve(2*q=x^2+3*y^2)},list);ns:=nops(sol);d:=sqrt((a^2+
> b^2+c^2)/3);
> ef:=0;
> for i from 1 to ns do 
>  if ef=0 then
>    r:=rhs(sol[i][1]); s:=rhs(sol[i][2]);
>    uu:=(s^2+3*r^2-2*q)^2;   if uu>0 then t:=s;s:=r;r:=t;fi; 
>   mz:=(r-s)/2;nz:=r;mw:=r;nw:=(r+s)/2;
>  
> mx:=-(d*b*(3*r+s)+a*c*(r-s))/(2*q);nx:=-(r*a*c+d*b*s)/q;my:=(d*a*(3*r+
> s)-b*c*(r-s))/(2*q);ny:=-(r*b*c-d*a*s)/q;
>   mu:=-(r*a*c+d*b*s)/q;nu:=-(d*b*(s-3*r)+a*c*(r+s))/(2*q);
>   mv:=(d*a*s-r*b*c)/q;nv:=-(d*a*(3*r-s)+b*c*(r+s))/(2*q);
>   
>  if mx=floor(mx) and nx=floor(nx) and my=floor(my) and ny=floor(ny)
>  and mu=floor(mu) and nu=floor(nu) and mv=floor(mv) and nv=floor(nv)
>  then  
> u:=mu*m-nu*n;v:=mv*m-nv*n;w:=mw*m-nw*n;
> x:=mx*m-nx*n;y:=my*m-ny*n;z:=mz*m-nz*n;
> om:=[[u,v,w],[x,y,z]]; ef:=1; fi;fi;
> od;
> om;
> end:
# Finding the triangles of minimal nonnegative coordinates in terms of 
# s,a,b,c and p and the main program 
> minimaltr:=proc(s,a,b,c,stopp)
> local
> i,z,u,nt,d,m,n,T,alpha,beta,gamma,tr,out,L,tri,noft,tria,orb,avb,lengt
> h,lengthn;
> d:=sqrt((a^2+b^2+c^2)/3);
> z:=s/d^2;u:=convert({isolve(z=q^2-q*r+r^2)},list);nt:=nops(u);
> for i from 1 to nt do
> T:=findpar(a,b,c,rhs(u[i][1]),rhs(u[i][2]));
> alpha:=min(T[1][1],T[2][1],0);
> beta:=min(T[1][2],T[2][2],0);
> gamma:=min(T[1][3],T[2][3],0);
> tr[i]:={[T[1][1]-alpha,T[1][2]-beta,T[1][3]-gamma],
> [T[2][1]-alpha,T[2][2]-beta,T[2][3]-gamma],[-alpha,-beta,-gamma]};
> out[i]:=max(tr[i][1][1],tr[i][1][2],tr[i][1][3],
> tr[i][2][1],tr[i][2][2],tr[i][2][3],> tr[i][3][1],tr[i][3][2],tr[i][3][3]);
> od;
> L:=sort([seq(out[i],i=1..nt)]);
> tri:={};
> for i from 1 to nt do
> if out[i]<= stopp then tri:=tri union {tr[i]};fi;
> od;
> tri:=convert(tri,list);
> tria:={};
> if nops(tri)>0 then
> noft:=nops(tri);tria:={tri[1]};
> orb:=transl(tri[1]);
>  for i from 1 to noft do
>       avb:=evalb(tri[i] in orb);
>   if avb=false  then
>       orb:=orb union transl(tri[i]);
>       tria:=tria union {tri[i]};
>     fi;
>   od;
> fi;
> tria;
> end:
> main:=proc(p)
> local i,j,k,s,nos,div,nod,nop,sol,x,netr,noft,l,z;
> netr:=0;
>     s:=sides(p);nos:=nops(s);
>  for i from 1 to nos do
>     div:=dkl(s[i]);nod:=nops(div);
>        for j from 1 to nod do
>          sol:=abcsol(div[j]);nop:=nops(sol);        
>         for k from 1 to nop do  
>          x:=minimaltr(s[i],sol[k][1],sol[k][2],sol[k][3],p); 
>          noft:=nops(x); if noft>=1 then
>           for l from 1 to noft do
>              z:=notrincn(x[l],p);
>              netr:=netr+z; print(checkeq(x[l]),sol[k],x[l],z,netr);
>           od;fi;
>      od;
>     od;
> od;
> netr;
> end:
# Example
> main(9);
 / (1/2)\                                                           
{ 2      }, [1, 1, 1], {[0, 1, 1], [1, 0, 1], [1, 1, 0]}, 5832, 5832
 \      /                                                           
 / (1/2)\                                                           
{ 6      }, [1, 1, 1], {[0, 1, 2], [2, 0, 1], [1, 2, 0]}, 4096, 9928
 \      /                                                           
 /   (1/2)\                                                            
{ 2 2      }, [1, 1, 1], {[2, 0, 2], [0, 2, 2], [2, 2, 0]}, 4096, 14024
 \        /                                                            
 /  (1/2)\                                                            
{ 14      }, [1, 1, 1], {[3, 1, 0], [0, 3, 1], [1, 0, 3]}, 5488, 19512
 \       /                                                            
 /   (1/2)\                                                            
{ 3 2      }, [1, 1, 1], {[0, 3, 3], [3, 0, 3], [3, 3, 0]}, 2744, 22256
 \        /                                                            
 /   (1/2)\                                                            
{ 3 2      }, [1, 1, 5], {[0, 3, 1], [3, 0, 1], [4, 4, 0]}, 7776, 30032
 \        /                                                            
 /   (1/2)\                                                            
{ 2 6      }, [1, 1, 1], {[2, 0, 4], [4, 2, 0], [0, 4, 2]}, 1728, 31760
 \        /                                                            
 /  (1/2)\                                                            
{ 26      }, [1, 1, 1], {[0, 4, 1], [1, 0, 4], [4, 1, 0]}, 3456, 35216
 \       /                                                            
 /   (1/2)\                                                            
{ 4 2      }, [1, 1, 1], {[0, 0, 4], [4, 0, 0], [0, 4, 0]}, 1728, 36944
 \        /                                                            
 /  (1/2)\                                                            
{ 38      }, [1, 1, 1], {[3, 0, 5], [0, 5, 3], [5, 3, 0]}, 2000, 38944
 \       /                                                            
 /  (1/2)\                                                            
{ 42      }, [1, 1, 1], {[5, 0, 4], [4, 5, 0], [0, 4, 5]}, 2000, 40944
 \       /                                                            
 /   (1/2)\                                                            
{ 5 2      }, [1, 1, 1], {[5, 0, 0], [0, 5, 0], [0, 0, 5]}, 1000, 41944
 \        /                                                            
 /   (1/2)\                                                            
{ 5 2      }, [1, 5, 7], {[4, 0, 4], [7, 5, 0], [0, 5, 1]}, 4320, 46264
 \        /                                                            
 /   (1/2)\                                                           
{ 3 6      }, [1, 1, 1], {[3, 6, 0], [6, 0, 3], [0, 3, 6]}, 512, 46776
 \        /                                                           
 /   (1/2)\                                                            
{ 3 6      }, [1, 1, 5], {[5, 7, 0], [7, 0, 1], [0, 2, 2]}, 1728, 48504
 \        /                                                            
 /    (1/2)\                                                            
{ 2 14      }, [1, 1, 1], {[6, 0, 2], [0, 2, 6], [2, 6, 0]}, 1024, 49528
 \         /                                                            
 /  (1/2)\                                                            
{ 62      }, [1, 1, 1], {[6, 1, 0], [0, 6, 1], [1, 0, 6]}, 1024, 50552
 \       /                                                            
 /   (1/2)\                                                           
{ 6 2      }, [1, 1, 1], {[0, 0, 6], [6, 0, 0], [0, 6, 0]}, 512, 51064
 \        /                                                           
 /   (1/2)\                                                           
{ 6 2      }, [1, 1, 5], {[0, 6, 2], [8, 8, 0], [6, 0, 2]}, 768, 51832
 \        /                                                           
 /  (1/2)\                                                           
{ 74      }, [1, 1, 1], {[4, 7, 0], [0, 4, 7], [7, 0, 4]}, 432, 52264
 \       /                                                           
 /  (1/2)\                                                           
{ 78      }, [1, 1, 1], {[7, 0, 2], [0, 2, 7], [2, 7, 0]}, 432, 52696
 \       /                                                           
 /  (1/2)\                                                           
{ 86      }, [1, 1, 1], {[7, 0, 1], [1, 7, 0], [0, 1, 7]}, 432, 53128
 \       /                                                           
 /   (1/2)\                                                          
{ 4 6      }, [1, 1, 1], {[0, 4, 8], [8, 0, 4], [4, 8, 0]}, 64, 53192
 \        /                                                          
 /   (1/2)\                                                           
{ 7 2      }, [1, 1, 1], {[8, 3, 0], [3, 0, 8], [0, 8, 3]}, 128, 53320
 \        /                                                           
 /   (1/2)\                                                           
{ 7 2      }, [1, 1, 1], {[7, 0, 7], [0, 7, 7], [7, 7, 0]}, 216, 53536
 \        /                                                           
 /   (1/2)\                                                            
{ 7 2      }, [1, 5, 11], {[0, 9, 0], [9, 5, 1], [1, 0, 4]}, 288, 53824
 \        /                                                            
 /    (1/2)\                                                           
{ 2 26      }, [1, 1, 1], {[8, 6, 0], [0, 8, 6], [6, 0, 8]}, 128, 53952
 \         /                                                           
 /   (1/2)\                                                           
{ 114      }, [1, 1, 1], {[8, 7, 0], [7, 0, 8], [0, 8, 7]}, 128, 54080
 \        /                                                           
 /   (1/2)\                                                          
{ 122      }, [1, 1, 1], {[4, 9, 0], [9, 0, 4], [0, 4, 9]}, 16, 54096
 \        /                                                          
 /    (1/2)\                                                          
{ 3 14      }, [1, 1, 1], {[0, 3, 9], [3, 9, 0], [9, 0, 3]}, 16, 54112
 \         /                                                          
 /   (1/2)\                                                          
{ 8 2      }, [1, 1, 1], {[8, 0, 8], [0, 8, 8], [8, 8, 0]}, 64, 54176
 \        /                                                          
 /   (1/2)\                                                          
{ 134      }, [1, 1, 1], {[0, 9, 2], [9, 2, 0], [2, 0, 9]}, 16, 54192
 \        /                                                          
 /   (1/2)\                                                          
{ 146      }, [1, 1, 1], {[1, 0, 9], [9, 1, 0], [0, 9, 1]}, 16, 54208
 \        /                                                          
 /   (1/2)\                                                         
{ 9 2      }, [1, 1, 1], {[0, 9, 0], [9, 0, 0], [0, 0, 9]}, 8, 54216
 \        /                                                         
                             54216

# The program lists the possible side lengths, some normal vector to a
# plane containing equilateral triangles given by a parameterization
# formula based on the normal vector, a so called minimal generating
# triangle, the number of equilateral triangles that this one generates
# under all cube symmetries and the tally number up to that particular
# side length.