Mathematica implementation of Lempel method for constructing Costas arrays

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

How to fix my code? Thanks in advance.

(* Costas Arrays with Lempel Method in Mathematica *)
(* Converted from C++ code by Fatih GULEC *)

(* Helper function to find prime factorization p^k *)
PrimeFactor[q_] := Module[{factors, p, k},
  factors = FactorInteger[q];
  If[Length[factors] == 1,
    {p, k} = factors[[1]];
    {p, k},
    Print["Error: Input must be a prime power"];
    {0, 0}
  ]
]

(* Convert decimal representation to polynomial coefficients *)
DecimalToPolyCoeffs[dec_, p_, degree_] := Module[{coeffs, temp},
  coeffs = {};
  temp = dec;
  Do[
    AppendTo[coeffs, Mod[temp, p]];
    temp = Floor[temp/p];
  , {i, 1, degree + 1}];
  coeffs
]

(* Convert polynomial coefficients to Mathematica polynomial *)
CoeffListToPoly[coeffs_, x_] := Module[{},
  Sum[coeffs[[i]] x^(i-1), {i, 1, Length[coeffs]}]
]

(* Check if polynomial is primitive using manual implementation *)
IsPrimitivePolynomial[coeffs_, p_] := Module[{
  degree, poly, x, isIrreducible, isPrimitive, order, testPoly, remainder
},
  degree = Length[coeffs] - 1;
  poly = CoeffListToPoly[coeffs, x];
  
  (* First check if it's irreducible by trying to factor *)
  isIrreducible = IrreduciblePolynomialQ[poly, Modulus -> p];
  
  If[!isIrreducible, Return[False]];
  
  (* Check if it's primitive by verifying the multiplicative order *)
  order = p^degree - 1;
  testPoly = x;
  
  (* Check if x^order ≡ 1 (mod poly) and x^d ≢ 1 for proper divisors d *)
  Do[
    If[order > 1 && Mod[order, d] == 0,
      remainder = PolynomialMod[testPoly^(order/d), poly, Modulus -> p];
      If[remainder == 1, Return[False]];
    ];
  , {d, 2, Floor[Sqrt[order]]}];
  
  (* Final check: x^order ≡ 1 (mod poly) *)
  remainder = PolynomialMod[testPoly^order, poly, Modulus -> p];
  remainder == 1
]

(* Find primitive polynomials of degree k over GF(p) *)
FindPrimitivePolynomials[p_, k_] := Module[{
  primitives, testDec, testEnd, coeffs, isValid
},
  primitives = {};
  testEnd = 2*p^k - 1;
  testDec = p^k + 1;
  
  While[testDec <= testEnd,
    If[Mod[testDec, p] != 0, (* Not divisible by X *)
      coeffs = DecimalToPolyCoeffs[testDec, p, k];
      
      (* Check if leading coefficient is non-zero *)
      If[coeffs[[-1]] != 0,
        isValid = IsPrimitivePolynomial[coeffs, p];
        If[isValid,
          AppendTo[primitives, coeffs];
        ];
      ];
    ];
    testDec++;
  ];
  
  Print["Number of primitive polynomials found: ", Length[primitives]];
  primitives
]

(* Generate field elements using primitive polynomial *)
GenerateFieldElements[p_, k_, primPolyCoeffs_] := Module[{
  elements, alpha, temp, i, j, q
},
  q = p^k;
  elements = Table[Table[0, {k}], {q}];
  
  (* Element 0: [0,0,...,0] *)
  elements[[1]] = Table[0, {k}];
  
  (* Element 1: [1,0,...,0] *)
  elements[[2]] = Table[0, {k}];
  elements[[2, 1]] = 1;
  
  (* Elements α, α^2, ..., α^(k-1) *)
  For[i = 3, i <= k + 1, i++,
    elements[[i]] = Table[0, {k}];
    If[i - 2 <= k, elements[[i, i - 1]] = 1];
  ];
  
  (* Generate higher powers α^k, α^(k+1), ... *)
  For[i = k + 2, i <= q, i++,
    (* Multiply previous element by α *)
    temp = Table[0, {k}];
    
    (* Shift left (multiply by α) *)
    For[j = 1, j < k, j++,
      temp[[j + 1]] = elements[[i - 1, j]];
    ];
    
    (* Handle overflow using primitive polynomial *)
    If[elements[[i - 1, k]] != 0,
      For[j = 1, j <= k, j++,
        temp[[j]] = Mod[temp[[j]] - elements[[i - 1, k]] * primPolyCoeffs[[j]], p];
      ];
    ];
    
    elements[[i]] = temp;
  ];
  
  elements
]

(* Generate Costas Array using Lempel method *)
GenerateCostasArray[primPolyCoeffs_, p_, k_] := Module[{
  q, elements, costasArray, i, j, oneMinusAi, found, match
},
  q = p^k;
  elements = GenerateFieldElements[p, k, primPolyCoeffs];
  costasArray = {};
  
  Print["Generating Costas array for primitive polynomial coefficients: ", primPolyCoeffs];
  Print["Field size q = ", q];
  
  (* Lempel algorithm: find j such that α^i + α^j = 1 *)
  For[i = 2, i < q, i++,
    (* Calculate 1 - α^i *)
    oneMinusAi = Table[0, {k}];
    oneMinusAi[[1]] = 1; (* Start with 1 *)
    For[j = 1, j <= k, j++,
      oneMinusAi[[j]] = Mod[oneMinusAi[[j]] - elements[[i, j]], p];
    ];
    
    (* Find matching element *)
    found = False;
    For[j = 2, j < q && !found, j++,
      match = True;
      For[l = 1, l <= k && match, l++,
        If[elements[[j, l]] != oneMinusAi[[l]], match = False];
      ];
      If[match,
        AppendTo[costasArray, j - 1];
        found = True;
      ];
    ];
    
    If[!found,
      Print["Warning: Could not find matching element for i = ", i];
    ];
  ];
  
  costasArray
]

(* Check if (q-3)rd order Costas array can be constructed *)
CheckQMinus3CostasArray[costasArray_, q_, k_] := Module[{},
  If[k == 1 && Length[costasArray] >= q - 2 && Last[costasArray] == q - 2,
    Print["(q-3)rd order Costas Array: ", Take[costasArray, q - 3]];
    Take[costasArray, q - 3],
    Print["(q-3)rd order Costas Array cannot be constructed."];
    {}
  ]
]

(* Main function to find Costas Arrays *)
FindCostasArrays[q_] := Module[{p, k, numPrim, primitives, costasArrays, i, costasArray},
  Print["Finding Costas Arrays for q = ", q];
  
  {p, k} = PrimeFactor[q];
  If[p == 0, Return[{}]];
  
  Print["p = ", p, ", k = ", k];
  
  numPrim = EulerPhi[q - 1]/k;
  Print["Expected number of primitive elements: ", numPrim];
  
  primitives = FindPrimitivePolynomials[p, k];
  costasArrays = {};
  
  Print["\nGenerating Costas Arrays:"];
  Print["========================"];
  
  For[i = 1, i <= Length[primitives], i++,
    Print["\nPrimitive polynomial ", i, " coefficients: ", primitives[[i]]];
    costasArray = GenerateCostasArray[primitives[[i]], p, k];
    
    If[Length[costasArray] > 0,
      Print["(q-2)nd order Costas Array: ", costasArray];
      AppendTo[costasArrays, costasArray];
      
      (* Check for (q-3)rd order *)
      CheckQMinus3CostasArray[costasArray, q, k];
    ];
  ];
  
  costasArrays
]

(* Visualize Costas Array *)
VisualizeCostasArray[costasArray_] := Module[{n, positions, plot},
  n = Length[costasArray];
  positions = Table[{i, costasArray[[i]]}, {i, 1, n}];
  
  plot = ListPlot[positions, 
    PlotStyle -> {Red, PointSize[0.02]},
    GridLines -> {Range[0, n + 1], Range[0, Max[costasArray] + 1]},
    GridLinesStyle -> LightGray,
    PlotRange -> {{0, n + 1}, {0, Max[costasArray] + 1}},
    AspectRatio -> 1,
    Frame -> True,
    FrameLabel -> {"Column", "Row"},
    PlotLabel -> "Costas Array Visualization"
  ];
  
  Print["Costas Array: ", costasArray];
  plot
]

(* Verify Costas property using triangular difference table *)
VerifyCostasProperty[costasArray_] := Module[{n, diffTable, allDistinct, k, i},
  n = Length[costasArray];
  If[n <= 1, Return[True]]; (* trivial cases *)
  
  Print["Verifying Costas array: ", costasArray];
  
  (* Generate triangular difference table *)
  diffTable = Table[
    Table[costasArray[[i + k]] - costasArray[[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:"];
  Do[
    Print["Row ", k, " (separation ", k, "): ", diffTable[[k]]] ,
    {k, 1, Length[diffTable]}
  ];
  Print["All rows have distinct elements: ", allDistinct];
  Print[""];
  
  If[allDistinct,
    Print["✓ Costas property verified successfully!"];
  , 
    Print["✗ Costas property violated!"];
  ];
  
  allDistinct
]

(* Additional verification function - checks if permutation is valid Costas array *)
IsCostasArray[perm_] := Module[{n, diffTable, allDistinct},
  n = Length[perm];
  If[n <= 1, Return[True]]; (* trivial cases *)
  
  (* Check if it's a valid permutation first *)
  If[Sort[perm] != Range[Max[perm]],
    Print["Error: Input is not a valid permutation"];
    Return[False];
  ];
  
  (* 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);
  
  allDistinct
]

(* Example usage *)
Print["Costas Arrays Generator using Lempel Method"];
Print["==========================================="];

(* Example: Find Costas arrays for q = 5 *)
q = 5;
costasArrays = FindCostasArrays[q];

Print["\nSummary:"];
Print["Found ", Length[costasArrays], " Costas arrays for q = ", q];

Print["\nAdditional functions available:"];
Print["- VisualizeCostasArray[array] : Plot the Costas array"];
Print["- VerifyCostasProperty[array] : Verify using triangular difference table"];
Print["- IsCostasArray[perm] : Check if permutation is a valid Costas array"];
Print["- FindCostasArrays[q] : Find all Costas arrays for given q"];

(* Example verification of known small Costas arrays *)
Print["\n=== Testing with Known Costas Arrays ==="];


(* Example visualization if arrays were found *)
If[Length[costasArrays] > 0,
  Print["\nVisualizing and verifying first generated Costas array:"];
  VisualizeCostasArray[costasArrays[[1]]];
  VerifyCostasProperty[costasArrays[[1]]];
]