Solving for a unique set of solutions

Question

I have the following code:

In[1]:=n = 6561;
Solve[{n == a^2 + b^2 + c^2, n == d^2 + e^2 + f^2, 
  n == g^2 + h^2 + i^2, n == a^2 + d^2 + g^2, n == b^2 + e^2 + h^2, 
  n == c^2 + f^2 + i^2, 
  1 <= a <= n && 1 <= b <= n && 1 <= c <= n && 1 <= d <= n && 
   1 <= e <= n && 1 <= f <= n && 1 <= g <= n && 1 <= h <= n && 
   1 <= i <= n}, {a, b, c, d, e, f, g, h, i}, PositiveIntegers]

And I want to do two things to this code:

1. How can I edit the code that it only gives a unique set of solutions, such that: $$a\ne b\ne c\ne d\ne e\ne f\ne g\ne h\ne i$$
2. How can I make a loop that tests different values of $n$ for a unique solution?

Thanks for any advice.

Answer 1

n = 6561;

Using PowersRepresentations and filtering:

(powersreps = Select[DuplicateFreeQ] @ PowersRepresentations[n, 3, 2];
 sPR = Sort[Join @@@ Join @@ (Permutations /@ 
   Select[SubsetQ[powersreps, Sort /@ Transpose[#]] &] @
   Select[DuplicateFreeQ @* Flatten] @
   Subsets[Join @@ (Permutations /@ powersreps), {3}])];) // AbsoluteTiming // First

2.21796 

Note: We could also use IntegerPartitions instead of PowersRepresentations:

 Sort[Sort /@ Select[DuplicateFreeQ] @ Sqrt @ 
    IntegerPartitions[n, {3}, Range[Sqrt @ n]^2]] == powersreps

True

Compare with using Solve and filtering:

vars = {a, b, c, d, e, f, g, h, i};
eqns1 = {n == a^2 + b^2 + c^2, n == d^2 + e^2 + f^2, n == g^2 + h^2 + i^2};
eqns2 = {n == a^2 + d^2 + g^2, n == b^2 + e^2 + h^2, n == c^2 + f^2 + i^2};
rangeconstraints = {And @@ Thread[1 <= vars < n]};

sOP = Select[DuplicateFreeQ] @  Values @ 
    Solve[Join[eqns1, eqns2, rangeconstraints], vars, Integers]; // 
   AbsoluteTiming // First

28.9619 

sPR == sOP

True

Update: A much faster method using FindInstance to find a single instance and processing it to generate the desired solution:

constraints = {And @@ Join[Thread[1 <= vars < n], {a < b < c, a != h}]};

(singleinstance = Partition[First @ Values @
    FindInstance[Join[eqns1, eqns2, constraints], vars, Integers], 3];
 sFI = Sort[Join @@@ Join @@ Transpose @* Map[Permutations] /@ 
    (Join @@ (Permutations /@ {#,  Transpose @ #} & @ singleinstance))];) // 
   AbsoluteTiming // First

0.183416 

 sFI == sOP

True

Answer 2

Your code gives 600 solutions.

n = 6561;
sol = {a, b, c, d, e, f, g, h, i} /.Solve[{n == a^2 + b^2 + c^2, n == d^2 + e^2 + f^2, 
n == g^2 + h^2 + i^2, n == a^2 + d^2 + g^2, n == b^2 + e^2 + h^2, 
n == c^2 + f^2 + i^2, 
1 <= a <= n && 1 <= b <= n && 1 <= c <= n && 1 <= d <= n && 
1 <= e <= n && 1 <= f <= n && 1 <= g <= n && 1 <= h <= n && 
1 <= i <= n}, {a, b, c, d, e, f, g, h, i}, PositiveIntegers];
Length[sol](*600*)

Reduce to unique solutions

Select[sol, Length[Union[#]] == 9 &]

evaluates 72!

Answer 3

Select[sol, DuplicateFreeQ]

Or

Pick[sol, DuplicateFreeQ /@ sol]