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[]