Sequence A102698 (* Mathematica Package *)

(* Created by the Wolfram Workbench Aug 16, 2007 *)

(* Translated by Rodrigo A. Obando from Maple code written by Eugen Ionascu *) 

BeginPackage["CountingEqTri`"]
(* Exported symbols added here with SymbolName::usage *) 
Sides::usage = "Sides[n] calculates the sides, l, of an equilateral triangle with integer coordinates through the map l^2/2 that will fit in {0,1,2,...,n} (via the characterization l=sqrt{2(m^2-mn+n^2)} and Gauss characterization of the numbers # m^2-mn+n^2 as those that do not contain in the prime factorization # odd powers of 2 or of primes of the form 6k-1.";
DKL::usage = "DKL[side] calculates the odd values of d that 'divide' a certain length.";
AddVector::usage = "AddVector[T, v] adds vector v to triangle T.";
ApplyPermutations::usage = "ApplyPermutations[p, mx] applies 48 permutations to p using mx.";
Orbit::usage = "Orbit[T] Finds the orbit of T.";
Transl::usage = "Transl[T] translates T";
Inters::usage = "Inters[T] intersects T";
Calc::usage = "Calc[T, n] checks the formula.";
InterSch::usage = "InterSch[T] Finding the intersection along translations along {0, 0, 1} and {0, 1, 1}.";
NoTriCn::usage = "NoTriCn[T, n] Calculates the number of triangles fitting in {0, 1, 2, ..., n}^3.";
CheckEq::usage = "CheckEq[T] Checks the lengths of the sides of the triangle T.";
ABCSol::usage = "ABCSol[d] finds the solutions of a^2+b^2+c^2=3*d^2 given d.";
FindPar::usage = "FindPar[a, b, c, m, n] finds parametrizations of a given solution (a, b, c).";
MinimalTr::usage = "MinimalTr[s, a, b, c, stopp] for finding the triangles of minimal nonnegative coordinates in terms of s, a, b, c and p.";
FindTotalEqTriangles::usage = "FindTotalEqTriangles[p, lastside, nuptols, detail] Gets the total number of equilateral triangles in a cube of size p (can be restarted in the middle). If detail = True it prints detailed output.";
FindTotalEquilateralTriangles::usage = "FindTotalEquilateralTriangles[dLength] Gets the total number of equilateral triangles in a cube of size dLength.";

Begin["`Private`"]
(* Implementation of the package *)
Sides[n_] := Module[{
	i, j, k, L, a, m, p, q, r, max},
		L = {};
		max = Floor[n^2];
   		For[i = 2, i &<= max, i++,
   			a = FactorInteger[i];
   			k = Length[a];
    		r = 0; 
    		For[j = 1, j &<= k, j++,
    			m = a[[j, 1]];
    			p = Mod[m, 6];
    			q = Mod[a[[j, 2]], 2];
    			If[r == 0 && (m == 2 || p == 5) && q == 1,
    				r = 1
    		]
    	]; 
    	If[r == 0,
    		L = Union[Join[L, {i}]]]
    	];
    	L = Union[Join[L, {1}]]
];

DKL[side_] := 
  Module[{
  	i, x, noft, div, y, y1, z},
  		x = Divisors[side]; 
  		noft = Length[x];
  		div = {}; 
  		For[i = 1, i &<= noft, i++,
  			z = Mod[x[[i]], 2]; 
    		If[z == 1,
    			y = side/x[[i]]^2;
    			y1 = Floor[y]; 
     			If[y == y1,
     				div = Union[Join[div, {x[[i]]}]]
     			]
     		]
     	];
     	div
];

AddVector[T_, v_] := Map[v + # &, T];

ApplyPermutations[p_, mx_] := Flatten[Map[Flatten[Outer[(Outer[(#1[[#2]]) &, {#1}, #2, 1]) &, {#}, Permutations[{1, 2, 3}], 1], 2] &, Map[#[[{1, 2, 3}]] &, Map[Map[If[Length[#] &> 1, -#, #] &, #] &, Map[#.{a, b, c, -1} &, Map[Transpose[{{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, #}] &, Map[# &, Table[Join[d IntegerDigits[n, 2, 3], {1}], {n, 0, 7}]]]]]]], 1] /. a -&> p[[1]] /. b -&> p[[2]] /. c -&> p[[3]] /. d -&> mx;

Orbit[p_] := Union[Map[Sort, Transpose[Map[ApplyPermutations[#, Max[Flatten[p]]] &, p]],{1}]];

Transl[T_] := Module[{
	mx, par},
	mx = Max[Flatten[T]] - Map[Max, Outer[T[[#2]][[#1]] &, {1, 2, 3}, {1, 2, 3}]];
	par = Map[Table[n, {n, 0, #}] &, mx];
	Union[Map[Sort,Union[Orbit[T], Flatten[Map[Orbit[AddVector[T, #]] &, Flatten[Outer[{#1,#2,#3} &, par[[1]], par[[2]], par[[3]]], 2]], 1]],{1}]]
 ];
 
Inters[T_] := Module[{S1, S2},
   S2 = Transl[T];
   S1 = Union[Map[AddVector[#, {0, 0, 1}] &, S2]];
   Intersection[S1, S2]
];
   
Calc[T_, n_] := Module[{
	S1, mx, par},
	mx = Max[Flatten[T]];
	par = Map[Table[m, {m, 0, #}] &, Table[n - mx, {p, 1, 3}]];
	S1 = Transl[T];
		{Length[Union[Flatten[Outer[AddVector[#1, #2] &, S1, Flatten[Outer[{#1, #2, #3} &, par[[1]], par[[2]], par[[3]]], 2], 1], 1]]], n - mx}
];

InterSch[T_] := Length[Intersection[Union[Map[AddVector[#, {0, 0, 1}] &, Transl[T]]], Union[Map[AddVector[#, {0, 1, 0}] &, Transl[T]]]]];
   
NoTriCn[T_, n_] := Module[{
	d, x, y, w},
	d = Max[Flatten[T]];
	x = Length[Transl[T]];
	y = Length[Inters[T]];
	w = InterSch[T];
	If[n &> d, (1 - d + n) (3 (d - n)^2 w + (1 - d + n)^2 x + 3 (d - n) (1 - d + n) y), x]
];

CheckEq[T_] := Union[Map[Sqrt[(#[[1, 1]] - #[[2, 1]])^2 + (#[[1, 2]] - #[[2, 2]])^2 + (#[[1, 3]] - #[[2, 3]])^2] &, Rest[Union[Map[If[#[[1]] != #[[2]], #] &, Union[Map[Sort, Tuples[T, 2]]]]]]]];

ABCSol[d_] := Module[{
	u, k, cd, i, j, x, y, sol = {}},
	For[i = 1, i &<= d, i++,
 		u = Reduce[3*d^2 - i^2 == x^2 + y^2, Integers];
 		k = Length[u];
 		For[j = 1, j &<= k, j++,
  			If[u[[j, 1, 2]] &>= i && u[[j, 2, 2]] &>= i,
   				cd = GCD[GCD[i, u[[j, 1, 2]]], u[[j, 2, 2]]];
   				If[cd == 1,
    				sol = Union[sol, {Sort[{i, u[[j, 1, 2]], u[[j, 2, 2]]}]}]
    			]
   			]
  		];
 	];
	sol   
];

FindPar[a_, b_, c_, m_, n_] := Module[{
	sol, q, x, y, ns, ef, i, r, s, d, mz, nz, mw, nw, mx, nx, my, ny, mu, nu, mv, nv, u, v, w, z, om={}},
	q = a^2 + b^2;
	sol = Reduce[2*q == x^2 + 3*y^2, {x, y}, Integers];
	ns = Length[sol];
	d = Sqrt[(a^2 + b^2 + c^2)/3];
	ef = 0;
	For[i = 1, i &<= ns, i++,
		If[ef == 0,
			s = sol[[i, 1, 2]];
			r = sol[[i, 2, 2]];
			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] && nx == Floor[nx] && my == Floor[my] && ny == Floor[ny] && mu == Floor[mu] && nu == Floor[nu] && mv == Floor[mv] && nv == Floor[nv],
				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
				]
			]
		];
		om
	];

MinimalTr[s_, a_, b_, c_, stopp_] := Module[{
	tr, out, d, z, u, q, r, nt, i, T, alpha, beta, gamma, k, tri, tria, noft, orb, avb},
	d = Sqrt[(a^2 + b^2 + c^2)/3];
	z = s/d^2;
	u = Reduce[z == q^2 - q*r + r^2, {q, r}, Integers];
	nt = Length[u];
	tr = Table[{},{k,1,nt}];
	out = Table[{},{k,1,nt}];
	For[i = 1, i &<= nt, i++,
		T = FindPar[a, b, c, u[[i, 1, 2]], u[[i, 2, 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]]}]
	];
	(*L = Sort[out];*)
	tri = {};
	For[i = 1, i &<= nt, i++,
		If[out[[i]] &<= stopp, tri = Union[tri, {tr[[i]]}]
		]
	];
	tri = Union[Map[Sort,tri,{1}]];
	tria = {};
	If[Length[tri] &> 0,
		noft = Length[tri];
		tria = {tri[[1]]};
		orb = Union[Map[Sort,Transl[tri[[1]]],{1}]];
		For[i = 1, i &<= noft, i++,
			avb = MemberQ[orb, tri[[i]]];
			If[avb == False,
				orb = Union[orb, Union[Map[Sort,Transl[tri[[i]]],{1}]]];
				tria = Union[tria, {tri[[i]]}];
			]
		]
	];
	Union[Map[Sort, tria, {1}]]
];
FindTotalEqTriangles[p_, lastside_, nuptols_, detail_: True] := Module[{
	netr, s, nos, i, div, nod, j, sol, nop, k, x, noft, z, l},
	netr = nuptols;
	s = Sides[p];
	nos = Length[s];
	If[detail,Print[s]];
	For[i = lastside, i &<= nos, i++,
 		div = DKL[s[[i]]];
 		nod = Length[div];
 		For[j = 1, j &<= nod, j++,
  			sol = ABCSol[div[[j]]];
  			nop = Length[sol];
  			For[k = 1, k &<= nop, k++,
   				x = MinimalTr[s[[i]], sol[[k, 1]], sol[[k, 2]], sol[[k, 3]], p];
   				noft = Length[x];
   				If[detail,Print["Number of Minimal Triangles: "&<&>ToString[noft]]];
   				If[noft &>= 1,
    				For[l = 1, l &<= noft, l++,
     					z = NoTriCn[x[[l]], p];
     					netr = netr + z;
     					If[detail,
     						Print[ToString[s[[i]]] &<&> ", " &<&> ToString[div] &<&> ", " &<&> ToString[sol[[k]]] &<&> ", " &<&> ToString[x[[l]]] &<&> ", " &<&> ToString[z] &<&> ", " &<&> ToString[netr] &<&> ", " &<&> ToString[i]]]],
     				If[detail,Print["No Minimal triangles when Doing " &<&> ToString[s[[i]]] &<&> " with div = " &<&> ToString[div] &<&> " and sol = " &<&> ToString[sol]]];
    			]
   			]
  		]
 	];
 	Print["n = "&<&>ToString[p]&<&>", eqt = "&<&>ToString[netr]];
 	netr
];

FindTotalEquilateralTriangles[dLength_] := Apply[Plus, 
 Flatten[Map[Map[NoTriCn[#, dLength] &, #] &, 
   Flatten[Map[
     Outer[MinimalTr[#1, #2[[1]], #2[[2]], #2[[3]], 
         dLength] &, {#[[1]]}, Flatten[ABCSol[#[[2]]], {1}], 1] &, 
     Flatten[Map[Flatten[Outer[{#1, #2} &, {#}, DKL[#]], {1, 2}] &, 
       Sides[dLength]], 1]], 2]]]];
       	
End[]

EndPackage[]