
(**
 * function[n, i] is the i-th n-ary function {0,1}^n -> {0,1}.
 * The index i may range from 1 to 2^(2^n).
 **)
function[n_Integer, i_Integer][vars__Integer] := Module[ {d, k},

  (* give the variables, construct the digit position *)
  d = Sum[{vars}[[k]] * 2^(k - 1), {k, 1, n}];

  (* return the $d$th digit in the binary representation of $i$. *)
  IntegerDigits[i, 2, 2^n][[d + 1]]
];

(**
 * Checks whether two boolean functions of three variables are equivalent.
 **)
equivQ[3, i_Integer, j_Integer] := Module[ {a, check},

  (* check whether there is one permutation/flip of arguments where
     the functions agree *)

  a = binstrings[3];

  check[1][x_, y_, z_] := Module[ {b, k},
    b = True;
    For[ k = 1, k <= Length[a], k++,
      If[ function[3, i][a[[k, 1]], a[[k, 2]], a[[k, 3]]] =!=
          function[3, j][a[[k, x]], a[[k, y]], a[[k, z]]], 
        b = False;
        Break[];
      ]
    ];
    Return[b]
  ];

  check[2][x_, y_, z_] := Module[ {b, k},
    b = True;
    For[ k = 1, k <= Length[a], k++,
      If[ function[3, i][a[[k, 1]], a[[k, 2]], a[[k, 3]]] =!=
          function[3, j][1-a[[k, x]], 1-a[[k, y]], 1-a[[k, z]]], 
          b = False;
          Break[];
      ]
    ];
    Return[b]
  ];

  check[x_,y_,z_] := check[1][x,y,z] || check[2][x,y,z];

  check[1,2,3] || check[1,3,2] || check[2,1,3] ||
  check[2,3,1] || check[3,1,2] || check[3,2,1] 
];

(**
 * Checks whether two boolean functions of two variables are equivalent.
 **)
equivQ[2, i_Integer, j_Integer] := Module[ {a, b, check, k},

  (* check whether there is one permutation/flip of arguments where
     the functions agree *)

  a = binstrings[2];

  check[1][x_, y_] := Module[ {b},
    b = True;
    For[ k = 1, k <= Length[a], k++,
      If[ function[2, i][a[[k, 1]], a[[k, 2]]] =!= 
          function[2, j][a[[k, x]], a[[k, y]]], 
          b = False;
      ]
    ];
    Return[b]
  ];

  check[2][x_, y_] := Module[ {b},
    b = True;
    For[ k = 1, k <= Length[a], k++,
      If[ function[2, i][a[[k, 1]], a[[k, 2]]] =!= 
          function[2, j][1 - a[[k, x]], 1 - a[[k, y]]], 
          b = False;
      ]
    ];
    Return[b]
  ];

  check[x_, y_] := check[1][x, y] || check[2][x, y];

  check[1,2] || check[2,1]
];

(**
 * Computes all equivalence classes of n-ary boolean functions.
 * In the resulting list, functions are represented by their indicees.
 **)
equivclasses[n_] := Module[ {list, i, done, j, k},

  list = {{1}};

  For[ i = 2, i <= 2^(2^n), i++,

    (* check if function[i] belongs to one of the equivalence classes 
    encountered so far. *)
    done = False;
    For[ j = 1, j <= Length[list], j++,
       If[ equivQ[n, i, First[list[[j]]]],
         AppendTo[list[[j]], i];
         done = True;
       ]
    ];
    If[ !done, AppendTo[list, {i}]];

  ];

  list
];

(* All binary strings of a given length *)
binstrings[n_] := Module[ {L},

   If[ n===1, 

      Return[{{0}, {1}}],

      L = binstrings[n - 1];
      Return[Join[Map[Prepend[#, 0]&, L], Map[Prepend[#, 1]&, L]]]
   ]
];

