I ported C++ implementation of Welch method for constructing Costas arrays to Mathematica.

Any feedback would be appreciated.

(* Welch method for generating Costas arrays *)

On["Assert"];

(* Find all primitive roots of a prime p *)
primitiveRoots[p_] := Module[{roots = {}},
  Do[
    If[MultiplicativeOrder[g, p] == p - 1, 
     AppendTo[roots, g]], 
    {g, 2, p - 1}];
  roots
]


(* Costas array verification function *)
(* Checks if a permutation is a valid Costas array using triangular difference table *)
isCostasArray[perm_] := Module[{n, diffTable, allDistinct},
  n = Length[perm];
  If[n <= 1, Return[True]]; (* trivial cases *)
  
  (* Generate triangular difference table *)
  diffTable = Table[
    Table[perm[[i + k]] - perm[[i]], {i, 1, n - k}],
    {k, 1, n - 1}
  ];
  
  (* Check if each row has distinct elements *)
  allDistinct = And @@ (DuplicateFreeQ[#] & /@ diffTable);
  
  (* Print difference table for debugging *)
  (*Print["Triangular difference table for ", perm, ":"];*)
  Do[
    Print["Row ", k, ": ", diffTable[[k]]] ,
    {k, 1, Length[diffTable]}
  ];
  Print["All rows have distinct elements: ", allDistinct];
  Print[""];
  
  allDistinct
]





(* Main computation *)
p = 13; (* Hard-coded prime *)

(* Verify p is prime *)
If[!PrimeQ[p] || p <= 2,
  Print[p, " is not a prime number greater than 2."];
  Abort[];
]

Print[p, " is a prime number."];

(* Find primitive roots *)
primitiveRootsList = primitiveRoots[p];
numPrimitiveRoots = Length[primitiveRootsList];

Print["Primitive roots:"];
Do[
  Print["Primitive root[", i - 1, "] = ", primitiveRootsList[[i]]],
  {i, 1, numPrimitiveRoots}
];

(* Generate Costas arrays *)
costasArrays = Table[
  Table[PowerMod[primitiveRootsList[[k]], j - 1, p], {j, 1, p - 1}],
  {k, 1, numPrimitiveRoots}
];

(* Print (p-1)th order Costas arrays *)
Print["(p-1)=", p - 1, "'th order Costas array:"];
Do[
  Print[costasArrays[[k]]],
  {k, 1, numPrimitiveRoots}
];

(* Generate (p-2)th order Costas arrays *)
(* Remove first element and subtract 1 from remaining elements *)
costas2Arrays = Table[
  Map[# - 1 &, Rest[costasArrays[[k]]]],
  {k, 1, numPrimitiveRoots}
];

Print["(p-2)=", p - 2, "'th order Costas array:"];
Do[
  Print[costas2Arrays[[k]]],
  {k, 1, numPrimitiveRoots}
];

(* Generate (p-3)th order Costas array if condition is met *)
(* Check if 2 is at position 2 in the first Costas array *)
If[Length[costasArrays] > 0 && Length[costasArrays[[1]]] > 1 && 
   costasArrays[[1]][[2]] == 2 && p > 3,
  (* Remove first two elements and subtract 2 from remaining *)
  costas3Array = Map[# - 2 &, Drop[costasArrays[[1]], 2]];
  Print["(p-3)=", p - 3, "'th order Costas array:"];
  Print[costas3Array];,
  (* Else *)
  Print["(p-3)'th order Costas array cannot be constructed."];
];











(* Verify that generated arrays are indeed Costas arrays *)
Print["\n=== VERIFICATION ==="];
Print["Verifying (p-1)th order Costas arrays:"];
Do[
  (*Print["Checking array ", k, ": ", costasArrays[[k]]];*)
  isValid = isCostasArray[costasArrays[[k]]];
  Print["Is valid Costas array: ", isValid];
  Print[""]; 
  Assert[isValid],
  {k, 1, numPrimitiveRoots}
];

Print["Verifying (p-2)th order Costas arrays:"];
Do[
  Print["Checking array ", k, ": ", costas2Arrays[[k]]];
  isValid = isCostasArray[costas2Arrays[[k]]];
  Print["Is valid Costas array: ", isValid];
  Print[""],
  {k, 1, numPrimitiveRoots}
];

(* Verify (p-3)th order array if it exists *)
If[Length[costasArrays] > 0 && Length[costasArrays[[1]]] > 1 && 
   costasArrays[[1]][[2]] == 2 && p > 3,
  costas3Array = Map[# - 2 &, Drop[costasArrays[[1]], 2]];
  Print["Verifying (p-3)th order Costas array:"];
  Print["Checking array: ", costas3Array];
  isValid = isCostasArray[costas3Array];
  (*Print["Is valid Costas array: ", isValid];
  Print[""];*)
  Assert[isValid];
];

(* Test with the example from the problem description *)
Print["=== EXAMPLE VERIFICATION ==="];
exampleArray = {1, 3, 4, 2, 5};
Print["Testing example array: ", exampleArray];
isCostasArray[exampleArray];

(* Display additional information *)
Print["\n=== SUMMARY ==="];
Print["Prime p = ", p];
Print["Number of primitive roots = ", numPrimitiveRoots];
Print["φ(p-1) = φ(", p - 1, ") = ", EulerPhi[p - 1]];
Print["Primitive roots: ", primitiveRootsList];

References

[1] J. H. van Lint and R. M. Wilson, "Latin squares which contain no repeated digrams," SIAM J. Appl. Math., vol. 35, no. 3, pp. 494–499, Nov. 1978, doi: 10.1137/1007035. [Online]. Available: https://www.sci-hub.ru/10.1137/1007035

[2] S. W. Golomb and H. Taylor, "Constructions and properties of Costas arrays," Proc. IEEE, vol. 72, no. 9, pp. 1143–1163, Sep. 1984, doi: 10.1109/PROC.1984.12994. [Online]. Available: https://ieeexplore.ieee.org/document/1457262

[3] A. Ardalani, "Contributions to the theory of Costas arrays," Ph.D. dissertation, Martin-Luther-Universität Halle-Wittenberg, Halle, Germany, 2023. [Online]. Available: https://opendata.uni-halle.de/bitstream/1981185920/110843/1/Ardalani_Ali_Dissertation_2023.pdf

[4] J. K. Beard, "A database of all known Costas arrays with 200 or fewer dots," [Online]. Available: https://web.archive.org/web/20210727010503/http://jameskbeard.com/jameskbeard/Costas_Arrays.html (accessed Jul. 3, 2025).