(*Mathematica input file for the notebook accompanying the article
  "Motzkin nuimbers and related sequences modulo powers of 2"*)

(*Coefficient reduction modulo 2,4,8,16*)
Rmod2[Expr_] := If[NumberQ[Expr], Mod[Expr, 2], Expr]
Rmod4[Expr_] := If[NumberQ[Expr], Mod[Expr, 4], Expr]
Rmod8[Expr_] := If[NumberQ[Expr], Mod[Expr, 8], Expr]
Rmod16[Expr_] := If[NumberQ[Expr], Mod[Expr, 16], Expr]

(*Computation of the generating function for Motzkin numbers
  modulo 4*)
M4 := Module[{Erg, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, a1, 
   a2, a3, a4, a5, a6, a7},
  F[z_] := A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3;
  A0[z_] := 1/(1 - z) + z^2/(1 - z)^3 + 2 a0[z];
  A1[z_] := 2 a1[z];
  A2[z_] := (1 - z)/z^2 + 2 a2[z];
  A3[z_] := 2 a3[z];
  Erg = F[z];
  For[k = 1, k <= 1, k++,
   Gl = Expand[z ^2 Erg^2 + (z - 1) Erg + 1];
   Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] + z^4/(1 - z^4))^2, 
     f[z]];
   Gl = Table[
      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
       3}] /. X -> 1;
   Var = Solve[
     Table[Gl[[ii]] == 0, {ii, 1, 4}], {a0[z], a1[z], a2[z], a3[z]}];
   Erg = Erg /. {a0[z] -> 
       Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] + 2 a0[z], 
      a1[z] -> Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] + 2 a1[z],
       a2[z] -> 
       Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] + 2 a2[z], 
      a3[z] -> Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] + 2 a3[z]};
   ];
  Collect[
    Map[Rmod4, Expand[
      Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, a3[z] -> 0}], 2], 
    f[z]] /. X -> 1
  ]

(*Computation of the generating function for Motzkin numbers
  modulo 16*)
M16 := Module[{Erg, X, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, 
   a1, a2, a3, a4, a5, a6, a7},
  F[z_] := 
   A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3 + A4[z] f[z]^4 + 
    A5[z] f[z]^5 + A6[z] f[z]^6 + A7[z] f[z]^7;
  A0[z_] := 1/(1 - z) + z^2/(1 - z)^3 + z^6/(1 - z)^7 + 2 a0[z];
  A1[z_] := 2 a1[z];
  A2[z_] := 2 a2[z];
  A3[z_] := 2 a3[z];
  A4[z_] := (1 - z)/z^2 + 2 a4[z];
  A5[z_] := 2 a5[z];
  A6[z_] := 2 a6[z];
  A7[z_] := 2 a7[z];
  Erg = F[z];
  For[k = 1, k <= 3, k++,
   Gl = Expand[z ^2 Erg^2 + (z - 1) Erg + 1];
   Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] - z^4/(1 - z^4))^4, 
     f[z]];
   Gl = Table[
      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
       7}] /. X -> 1;
   Var = Solve[
     Table[Gl[[ii]] == 0, {ii, 1, 8}], {a0[z], a1[z], a2[z], a3[z], 
      a4[z], a5[z], a6[z], a7[z]}];
   Erg = Erg /. {a0[
        z] -> (Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a0[z], 
      a1[z] -> (Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a1[z], 
      a2[z] -> (Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a2[z], 
      a3[z] -> (Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a3[z],
      a4[z] -> (Factor[Factor[Var[[1, 5, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a4[z], 
      a5[z] -> (Factor[Factor[Var[[1, 6, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a5[z], 
      a6[z] -> (Factor[Factor[Var[[1, 7, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a6[z], 
      a7[z] -> (Factor[Factor[Var[[1, 8, 2]]], Modulus -> 2] /. 
          z -> -z) + 2 a7[z]};
   ];
  Collect[
    Map[Rmod16, 
     Expand[Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, 
        a3[z] -> 0}], 2], f[z]] /. X -> 1
  ]

(*Computation of the generating function for Motzkin prefix numbers
  modulo 4*)
MP4 := Module[{Erg, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, a1, 
		 a2, a3, a4, a5, a6, a7},
    F[z_] := A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3;
    A0[z_] := 1/(1 - z) + z/(1 - z)^2 + z^3/(1 - z)^4 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := 1/z + 2 a2[z];
    A3[z_] := 2 a3[z];
    Erg = F[z];
    For[k = 1, k <= 1, k++,
	Gl = Expand[z (1 - 3 z) Erg^2 - (3 z - 1) Erg - 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] + z^4/(1 - z^4))^2, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
					3}] /. X -> 1;
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 4}], {a0[z], a1[z], a2[z], a3[z]}];
	Erg = Erg /. {a0[z] -> Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] + 2 a0[z], 
		      a1[z] -> Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] + 2 a1[z],
                      a2[z] -> Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] + 2 a2[z], 
		      a3[z] -> Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] + 2 a3[z]};
	];
  Collect[Map[Rmod4, Expand[Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, a3[z] -> 0}], 2], 
    f[z]] /. X -> 1
  ]

(*Computation of the generating function for Motzkin prefix numbers
  modulo 16*)
MP16 := Module[{Erg, X, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, 
		  a1, a2, a3, a4, a5, a6, a7},
    F[z_] := 
   A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3 + A4[z] f[z]^4 + 
    A5[z] f[z]^5 + A6[z] f[z]^6 + A7[z] f[z]^7;
    A0[z_] := 
    1/(1 - z) + z/(1 - z)^2 + z^3/(1 - z)^4 + z^7/(1 - z)^8 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := 2 a2[z];
    A3[z_] := 2 a3[z];
    A4[z_] := 1/z + 2 a4[z];
    A5[z_] := 2 a5[z];
    A6[z_] := 2 a6[z];
    A7[z_] := 2 a7[z];
    Erg = F[z];
    For[k = 1, k <= 3, k++,
	Gl = Expand[z (1 - 3 z) Erg^2 - (3 z - 1) Erg - 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] - z^4/(1 - z^4))^4, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
				 7}] /. X -> 1;
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 8}], {a0[z], a1[z], a2[z], a3[z], 
						  a4[z], a5[z], a6[z], a7[z]}];
	Erg = Erg /. {a0[z] -> (Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] /.z -> -z) + 2 a0[z], 
		      a1[z] -> (Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] /.z -> -z) + 2 a1[z], 
		      a2[z] -> (Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] /.z -> -z) + 2 a2[z], 
		      a3[z] -> (Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] /.z -> -z) + 2 a3[z],
		      a4[z] -> (Factor[Factor[Var[[1, 5, 2]]], Modulus -> 2] /.z -> -z) + 2 a4[z], 
		      a5[z] -> (Factor[Factor[Var[[1, 6, 2]]], Modulus -> 2] /.z -> -z) + 2 a5[z], 
		      a6[z] -> (Factor[Factor[Var[[1, 7, 2]]], Modulus -> 2] /.z -> -z) + 2 a6[z], 
		      a7[z] -> (Factor[Factor[Var[[1, 8, 2]]], Modulus -> 2] /.z -> -z) + 2 a7[z]};
	];
  Collect[Map[Rmod16, Expand[Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, 
		  a3[z] -> 0}], 2], f[z]] /. X -> 1
  ]

(*Computation of the generating function for Riordan numbers
  modulo 4*)
R4 := Module[{Erg, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, a1, 
		 a2, a3, a4, a5, a6, a7},
    F[z_] := A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3;
    A0[z_] := 1/(1 - z) + z/(1 - z)^2 + z^3/(1 - z)^4 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := 1/z + 2 a2[z];
    A3[z_] := 2 a3[z];
    Erg = F[z];
    For[k = 1, k <= 1, k++,
	Gl = Expand[z (1 + z) Erg^2 - (z + 1) Erg + 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] + z^4/(1 - z^4))^2, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
								    3}] /. X -> 1;
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 4}], {a0[z], a1[z], a2[z], a3[z]}];
	Erg = Erg /. {a0[z] -> 
		      Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] + 2 a0[z], 
		      a1[z] -> Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] + 2 a1[z],
       a2[z] -> 
		      Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] + 2 a2[z], 
		      a3[z] -> Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] + 2 a3[z]};
	];
  Collect[
	  Map[Rmod4, Expand[
			 Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, a3[z] -> 0}], 2], 
    f[z]] /. X -> 1
  ]

(*Computation of the generating function for Riordan numbers
  modulo 16*)
R16 := Module[{Erg, X, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, 
		  a1, a2, a3, a4, a5, a6, a7},
    F[z_] := 
   A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3 + A4[z] f[z]^4 + 
    A5[z] f[z]^5 + A6[z] f[z]^6 + A7[z] f[z]^7;
    A0[z_] := 
    1/(1 - z) + z/(1 - z)^2 + z^3/(1 - z)^4 + z^7/(1 - z)^8 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := 2 a2[z];
    A3[z_] := 2 a3[z];
    A4[z_] := 1/z + 2 a4[z];
    A5[z_] := 2 a5[z];
    A6[z_] := 2 a6[z];
    A7[z_] := 2 a7[z];
    Erg = F[z];
    For[k = 1, k <= 3, k++,
	Gl = Expand[z (1 + z) Erg^2 - (z + 1) Erg + 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] - z^4/(1 - z^4))^4, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
								    7}] /. X -> 1;
   
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 8}], {a0[z], a1[z], a2[z], a3[z], 
						  a4[z], a5[z], a6[z], a7[z]}];
	Erg = Erg /. {a0[
			 z] -> (Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a0[z], 
		      a1[z] -> (Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a1[z], 
		      a2[z] -> (Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a2[z], 
		      a3[z] -> (Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a3[z],
		      a4[z] -> (Factor[Factor[Var[[1, 5, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a4[z], 
		      a5[z] -> (Factor[Factor[Var[[1, 6, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a5[z], 
		      a6[z] -> (Factor[Factor[Var[[1, 7, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a6[z], 
		      a7[z] -> (Factor[Factor[Var[[1, 8, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a7[z]};
	];
  Collect[
	  Map[Rmod16, 
	      Expand[Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, 
		  a3[z] -> 0}], 2], f[z]] /. X -> 1
  ]

(*Computation of the generating function for hex tree numbers
  modulo 4*)
H4 := Module[{Erg, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, a1, 
		 a2, a3, a4, a5, a6, a7},
    F[z_] := A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3;
    A0[z_] := 1/(1 - z) + z^2/(1 - z)^3 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := (1 - z)/z^2 + 2 a2[z];
    A3[z_] := 2 a3[z];
    Erg = F[z];
    For[k = 1, k <= 1, k++,
	Gl = Expand[z ^2 Erg^2 + (3 z - 1) Erg + 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] + z^4/(1 - z^4))^2, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
								    3}] /. X -> 1;
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 4}], {a0[z], a1[z], a2[z], a3[z]}];
	Erg = Erg /. {a0[z] -> 
		      Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] + 2 a0[z], 
		      a1[z] -> Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] + 2 a1[z],
       a2[z] -> 
		      Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] + 2 a2[z], 
		      a3[z] -> Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] + 2 a3[z]};
	];
  Collect[
	  Map[Rmod4, Expand[
			 Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, a3[z] -> 0}], 2], 
    f[z]] /. X -> 1
  ]

(*Computation of the generating function for hex tree numbers
  modulo 16*)
H16 := Module[{Erg, X, Var, F, A0, A1, A2, A3, A4, A5, A6, A7, a0, 
		  a1, a2, a3, a4, a5, a6, a7},
    F[z_] := 
   A0[z] + A1[z] f[z] + A2[z] f[z]^2 + A3[z] f[z]^3 + A4[z] f[z]^4 + 
    A5[z] f[z]^5 + A6[z] f[z]^6 + A7[z] f[z]^7;
    A0[z_] := 1/(1 - z) + z^2/(1 - z)^3 + z^6/(1 - z)^7 + 2 a0[z];
    A1[z_] := 2 a1[z];
    A2[z_] := 2 a2[z];
    A3[z_] := 2 a3[z];
    A4[z_] := (1 - z)/z^2 + 2 a4[z];
    A5[z_] := 2 a5[z];
    A6[z_] := 2 a6[z];
    A7[z_] := 2 a7[z];
    Erg = F[z];
    For[k = 1, k <= 3, k++,
	Gl = Expand[z^2 Erg^2 + (3 z - 1) Erg + 1];
	Gl = PolynomialRemainder[Gl, (f[z]^2 + f[z] - z^4/(1 - z^4))^4, 
				 f[z]];
   Gl = Table[
	      Map[Rmod2, Coefficient[Expand[Gl/2^(k)*X], f[z], i], 2], {i, 0, 
								    7}] /. X -> 1;
   Var = Solve[
	       Table[Gl[[ii]] == 0, {ii, 1, 8}], {a0[z], a1[z], a2[z], a3[z], 
						  a4[z], a5[z], a6[z], a7[z]}];
	Erg = Erg /. {a0[
			 z] -> (Factor[Factor[Var[[1, 1, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a0[z], 
		      a1[z] -> (Factor[Factor[Var[[1, 2, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a1[z], 
		      a2[z] -> (Factor[Factor[Var[[1, 3, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a2[z], 
		      a3[z] -> (Factor[Factor[Var[[1, 4, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a3[z],
		      a4[z] -> (Factor[Factor[Var[[1, 5, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a4[z], 
		      a5[z] -> (Factor[Factor[Var[[1, 6, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a5[z], 
		      a6[z] -> (Factor[Factor[Var[[1, 7, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a6[z], 
		      a7[z] -> (Factor[Factor[Var[[1, 8, 2]]], Modulus -> 2] /. 
				z -> -z) + 2 a7[z]};
	];
  Collect[
	  Map[Rmod16, 
	      Expand[Erg*X /. {a0[z] -> 0, a1[z] -> 0, a2[z] -> 0, 
		  a3[z] -> 0}], 2], f[z]] /. X -> 1
  ]


