> 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.