(* To compute the equivalence classes for quartic MRS functions in n variables under permutations which preserve rotation symmetry enter the value of n below and run the program. Note that the program lines below must

be entered into Mathematica to produce a .nb file which can then be run *)

Clear[equivclasses]

(* This is the permutation in Th. 1.7 on p.196 of T. W. Cusick, Y. Cheon, Information Sciences 259 (2014), 192-211. This paper is \referenced as [TC14] below *)

mu[x_, sigtwo_] := (x - 1) (sigtwo - 1) + 1;

funcs[n_] := Cases[Tuples[{{1}, Range[2, Floor[n/4] + 1], Range[3, n - 2], Range[4, n - 1]}], x_ /;

((2*x[[2]] - 1 <= x[[3]] & x[[2]] - 1 <= x[[4]] - x[[3]] & x[[4]] <= n - x[[2]] + 1) & If[x[[2]] - 1 == x[[4]] - x[[3]], x[[3]] - x[[2]] <= n - x[[4]] + 1, True] || (Divisible[n, 4] & x[[2]] == n/4 + 1 & x[[3]] == n/2 + 1 & x[[4]] == 3 n/4 + 1))];

freduce[x_, n_] := Sort[

Mod[{Sort[x], {1, 2 + n - Sort[x][[4]], Sort[x][[2]] - Sort[x][[4]] + n + 1, Sort[x][[3]] - Sort[x][[4]] + n + 1}, {1, Sort[x][[4]] - Sort[x][[3]] + n + 1, 2 + n - Sort[x][[3]], Sort[x][[2]] - Sort[x][[3]] + n + 1}, {1, Sort[x][[3]] - Sort[x][[2]] + n + 1, Sort[x][[4]] - Sort[x][[2]] + n + 1, 2 + n - Sort[x][[2]]}}, n, 1]][[1]];

n = 27;

fun = funcs[n];

(* Initially, all functions are given their own equivalence class index in a hash map. The number of distinct functions is given in [Lemma 1.3, p. 5071, TC14] *)

For[i = 1, i <= Length[fun], i++, equivclasses[fun[[i]]] = i ];

For[sig2 = 3, sig2 <= n, sig2++,

If[GCD[sig2 - 1, n] == 1,

For[i = 1, i <= Length[fun], i++,

newfun = freduce[Mod[mu[fun[[i]], sig2] - 1, n] + 1, n];

(*Print["For \[Sigma](2)=",sig2,", ",fun[[i]]," ---> ", newfun];*)

(* Combine the two equivalence classes if they are different. Assign the larger numbered function to the equivalence class of the lower numbered function. We do not have to check the other functions because of the group structure [Theorem 1.8, p. 197, TC14] *)

newval = Min[{equivclasses[fun[[i]]], equivclasses[newfun]}]; equivclasses[fun[[i]]] = equivclasses[newfun] = newval;

]

]

]

(* Rebuild information the other way. We now have a mapping from a function to the number of its equivalence class but we want a mapping from the equivalence class index to the set of functions. *)

classes \ = ConstantArray[{}, Length[fun]]; (* Preallocation *)

For[i = 1, i <= Length[fun], i++, AppendTo[classes[[equivclasses[fun[[i]]]]], fun[[i]]];

]

classes = Cases[classes, x_ /; x != {}]; (* Remove unused indices *)

For[i = 1, i <= Length[classes], i++, Print["Class ", i, ", Size: ", Length[classes[[i]]], ": ", classes[[i]]]

]

(* Now the program lists the classes in the format shown below. The classes appear with their initial elements in lexicographic order *)

Class 1, Size: 9:

{{1,2,3,4},{1,2,14,15},{1,3,5,7},{1,4,11,21},{1,4,12,20},{1,5,9,13},{1,6,11,16},{1,6,12,17},{1,7,14,21}}

Class 44, Size: 2: {{1,4,10,19},{1,4,13,22}}