The answer is that I am using BooleanConvert, but needed an alternative version:

ANFtoArith[fun_] := With[{fun2 = Replace[ReplaceRepeated[
      fun, {x_ && y_ :> x*y, ! x_ :> x + 1,
       x_ \[Xor] y_ :> x + y, False -> 0, True -> 1}],
     {Function[x_] :> Mod[Expand[x], 2]}]}, Function[fun2]]

f2sBool = BooleanConvert[BooleanFunction[#], "ANF"] & /@ Tuples[{0, 1}, 2^2];
f2sArith = ANFtoArith /@ f2sBool
Out[]={0 &, Mod[1 + #1 + #2 + #1 #2, 2] &, Mod[#2 + #1 #2, 2] &, 
 Mod[1 + #1, 2] &, Mod[#1 + #1 #2, 2] &, Mod[1 + #2, 2] &, 
 Mod[#1 + #2, 2] &, Mod[1 + #1 #2, 2] &, Mod[#1 #2, 2] &, 
 Mod[1 + #1 + #2, 2] &, Mod[#2, 2] &, Mod[1 + #1 + #1 #2, 2] &, 
 Mod[#1, 2] &, Mod[1 + #2 + #1 #2, 2] &, Mod[#1 + #2 + #1 #2, 2] &, 1 &}

The full story is kind of funny actually... I got into metamath relatively late, actually at the very end. The assumption seemed to be "oh this guy is a physicist, what could he possibly know". Well, for one, "circuit design" and Boolean logic. They had a proof of De Morgan's theorem on the record I think it was taking hundreds of steps and I said, "No, it's just one step, Expand[...]":

Mod[Expand[ANFtoArith[And[Not[#1], Not[#2]]] @@ {a, b}], 2]
ANFtoArith[BooleanConvert[BooleanFunction[{0, 0, 0, 1}], "ANF"]] @@ {a, b}
Out[47]= Mod[1 + a + b + a b, 2]
Out[48]= Mod[1 + a + b + a b, 2] 

There was probably another implication that I was too small minded to fully understand the notion of how to count operations, but I also did some Bernstein-style analysis breaking down Expand to more elemental logic or arthmetic functions and counting repeated usages up to forty or fifty at max, not hundreds.

SW asked for this to be developed for WFR, then someone else demoted it to a Hobby project in his official notes. If the disagreement gets any bigger I guess we will have to ask Thomas Hales.

For now the priority is relatively low with hopes to have some more done by July, but generalizing ANF to other radix looks to be really not that easy, perhaps explaining why it isn't in NKS (but should be). If we figure it out, in the end, we could have a more powerful, more general version of BooleanConvert, which arguably could follow into the kernel.

Yea I agree find patterns in the pictures could be a way forward, but another issue I'm worried about is possible multiple solutions. Is there exactly one way to "functionalize a truth table", or do we need to worry about canonical choosing?

Hi Todd, nice typing, but did you mean "important columns" rather than "other stuff"... I'm glad I checked because I may not have noticed the similarity to my recent proof finding the Algebraic Normal Form for certain radix cellular automata:

TermsUpTo[vars_ : {__Slot}, nMax_
  ] := TermsUpTo[vars /. Slot -> S, nMax] /. S -> Slot

TermsUpTo[vars_ /; ! MemberQ[Head /@ vars, Slot], nMax_Integer
  ] := TermsUpTo[vars, Table[nMax, {Length[vars]}]]

TermsUpTo[vars_ /; ! MemberQ[Head /@ vars, Slot], nMax_ : {__Integer}
  ] := Reverse[SortBy[Union[Flatten[
     Outer[Times @@ DeleteCases[Times[#1, #2], 0] &,
      Flatten[Outer[MapThread[#1^#2 &,
          {vars, {##}}] &,
        Sequence @@ Range /@ nMax,
        1], Length[nMax] - 1],
      Tuples[{0, 1}, Length[nMax]], 1]]],
   Flatten[{#} /. {Power[x_, y_] :> Table[x, {y}],
       Times -> List}] &]]

CAtoANF[color_, radius_] := Function[Evaluate[
    Mod[TermsUpTo[Slot /@ Range[2 radius + 1], color - 1], color]]
   ] @@@ Tuples[Reverse[Range[0, color - 1]], 2 radius + 1]

With[{CAtoANFPrimes = CAtoANF[Prime[#], 1] & /@ Range[4]},
 Column[MapIndexed[
   Row[Show[#, ImageSize -> 400] & /@ {ArrayPlot[#1],
       ArrayPlot[Inverse[#1, Modulus -> Prime[#2[[1]]] ]]}, 
     Spacer[5]] &, CAtoANFPrimes]]]

Effectively a constructive proof of existence of ANF for first few primes, which is expected to extend to all primes. If you try this with composites mixed in, matrix inversion breaks, and it starts to look like Galois Field Theory will be necessary.

However, we can go immediately to four colors, because I fixed the proof in that one particular case:

CAtoANF4 = Function[Evaluate[
     Mod[TermsUpTo[Slot /@ Range[6], 1], 4]]
    ] @@@ Tuples[{1, 0}, 6];

Row[Show[#, ImageSize -> 400] & /@ {ArrayPlot[CAtoANF4],
   ArrayPlot[Inverse[CAtoANF4, Modulus -> 4]]}, Spacer[5]]

Notice similarity to your BitReverseOrder, but it's not exactly the same.

Anyways, I tried to go to radix six, but couldn't figure it out. Any ideas?

Do you know a reference for a proof that even says ANF exists for any radix?

I might lecture this summer how useful ANF is for complexity ordering search spaces... When you can cut something down from say $3^{27}$ to $3^{10}$ it becomes a lot easier to be thorough and actually get some results. Even if you stay at $3^{27}$ or higher, extra info about functional complexity can help building ocean maps for deep dive searches.