Counting solutions for a polynomial over GF(2)

I have written a codes gf2_sat_count/3 which counts all solutions of a given polynomioal over the well known GF2 (Galoi field {0, 1}). I am not familiar with this important field, but my interest was applying a similar method for the boolean algebra {true, false}, which, I think, works fine.

gf2_sat_count/3 works polynomials in ZDD.

I got this table, in which C is the number of orthogonal matricies of size n x n over GF(2).

n x n.      C (the number ortho(GF2,n)
------     --------------------------
1 x 1		1
2 x 2		2
3 x 3		6
4 x 4		48
5 x 5		720
6 x 6		23040

It took less thatn one hour for n = 6, most of which was for reading a polynomial condition of the orthogonality ( M * transpos M = the unit matrix ) in ZDD. The time for counting was negligible.

However, I remember I saw a different result table for the same problem. So I am not sure my table is correct. Anyway, this table was obtained as naive and direct application of the gf2_sat_count.

In addition, I have just noticed there seems to be a different SAT counting method for GF2 which may reduce much time for the time consuming reading polynomials. I have to check if
it really makes possible to read the orthogonal condition as a system of equations, avoiding unnecessary expansion to one huge polynomial.

These numbers seems to be correct.

From OEIS “A003053 Order of orthogonal group O(n, GF(2)).” (A003053 - OEIS). The sequence continues with 1, 2, 6, 48, 720, 23040, 1451520, 185794560, 47377612800, …

Thanks for your kindness. In fact, I was sure only fifty-to-fifty. Now I would like to try on the boolean ring for next, expecting similarity to GF2, but after making sure on the definition of that ring.

I am learning pldoc for the time in the future when I might show codes for one who asks me.

As sample, below is codes with %! comments. I am not sure for now if it fits the pldoc recommendation. The code is on gf2-solver for systems of polynomial equations over GF2. It is
an initial version, but I have no idea and plan to revise it. It seems difficult for me. I like more easier prolog programming.

EDIT: 2022-12-11
This code seems to have a bug. Applying orth(GF2, n), it works for n = 1, 2, 3 but not for n = 4, 5, 6. I have to debug it. However, for n = 6, it runs about 100 times faster than the correct version, which, I think, the idea of solving in a systrem of equations may work as I expected.
END

gf2_system_solver
		/**************************
		*     gf2_system_solver   *
		**************************/

%!	gf2_solve_system(+Vs, +Eqs, -Sols, +S) is det.
%	Vs : a list varialbles (indeterminates) of Eqs.
%	Eqs: a list [P1, P2, ..., Pn] as a system of equations
%		P1 = 0, P2 = 0, ...,, Pn = 0  of polynomials over GF2.
%	Sols is unified with a zdd of a family of all subsets U of Vs
%   such that U satisfies all Pi = 0.

% ?- #S, X<< (a+b), Y<< (x+y),
%	gf2_solve_system([a,b,x,y], [X, Y], Z, S), psa(Z, S), card(Z, C, S).

gf2_solve_system(Ls, Eqs, Sols, S):- zdd_power(Ls, Ini, S),
	gf2_solve_system(Ls, Eqs, Ini, Sols, S).
%
gf2_solve_system(_, [], Y, Y, _):-!.
gf2_solve_system(_, _, 0, 0, _):-!.
gf2_solve_system(Ls, [X|Xs], Y, Z, S):-
	solve(Ls, X, Sols, S),
	zdd_meet(Sols, Y, Y0, S),
	gf2_solve_system(Ls, Xs, Y0, Z, S).

%!	gf2_solve(+Vs, +X, -Sols, +S) is det.
%	Vs: indeterminates(variables) of a polynomial over GF2.
%	X : a zdd of the polynomial.
%	Sols is unified with a zdd of all solutions of the polynomial.

% ?- A = [], #S, X<<pow(A), gf2_solve([], 0, Y, S), psa(Y, S).
% ?- A = [], #S, X<<pow(A), gf2_solve([], 1, Y, S), psa(Y, S).
% ?- A = [c], #S, zdd_singleton(c, C, S), gf2_solve(A, C, Y, S), psa(Y, S).
% ?- A = [a, b], #S, gf2_solve(A, 0, Y, S), psa(Y, S).
% ?- A = [], #S, X<<pow(A), gf2_solve([], X, Y, S), psa(Y, S).
% ?- A = [a], #S, X<<pow(A), gf2_solve(A, X, Y, S), psa(Y, S).
% ?- #S, X<< +[a, b], psa(X, S), gf2_solve([a,b], X, Y, S), psa(Y, S).
% ?- #S, X<< *[a, b], psa(X, S), gf2_solve([a,b], X, Y, S), psa(Y, S).
% ?- #S, X<< +[a], psa(X, S), gf2_solve([a], X, Y, S), psa(Y, S).
% ?- #S, X<< *[a,b], psa(X, S), gf2_solve([a,b], X, Y, S), psa(Y, S).
% ?- #S, X<< (*a), psa(X, S), gf2_solve([a], X, Y, S), psa(Y, S).
% ?- #S, gf2_solve([a], 0, Y, S), psa(Y, S).
% ?- #S, X<< (a*b + 1), psa(X, S), gf2_solve([a,b], X, Y, S), psa(Y, S), card(Y, C, S).
% ?- #S, X<< (a*b + x*y), psa(X, S), gf2_solve([a,b], X, Y, S), psa(Y, S), card(Y, C, S).
% ?- #S, X<< (a*b + x*y), psa(X, S), gf2_solve([a,b,x,y], X, Y, S), psa(Y, S), card(Y, C, S).
% ?- #S, X<< (a+b+x+y), psa(X, S), gf2_solve([a,b,x,y], X, Y, S), psa(Y, S), card(Y, C, S).

gf2_solve(Ls, 0, Sols, S):-!, zdd_power(Ls, Sols, S).
gf2_solve(Ls, X, Sols, S):- select_constant(X, C, X0, S),
	(	C = 0 -> gf2_solve_eq0(Ls, X0, Sols, S)
	;	gf2_solve_eq0(Ls, X0, Sols0, S),
		zdd_power(Ls, All, S),
		zdd_subtr(All, Sols0, Sols, S)
	).
%
gf2_solve_eq0(Ls, 0, X, S):-!, zdd_power(Ls, X, S).
gf2_solve_eq0(_, 1, 0, _):-!.
gf2_solve_eq0(Ls, P, X, S):- memo(gf2_solve(P, Ls)-X, S),
	(	nonvar(X) -> true	%, write(.)	% no hits. Useless here ?
	;	cofact(P, t(A, L, R), S),
		divide_at(A, Ls, H, T),
		zdd_power(T, PowT, S),
		cofact(PowAT, t(A, PowT, PowT), S),
		(	gf2_solve_eq0(T, R, R0, S),
			cofact(AR_sols_zero, t(A, PowT, R0), S),
			gf2_solve_eq0(T, L, L0, S),
			cofact(L_sols_zero, t(A, L0, L0), S),
			zdd_meet(AR_sols_zero, L_sols_zero, Sols_zero, S),
			zdd_subtr(PowAT, AR_sols_zero, AR_sols_unit, S),
			zdd_subtr(PowAT, L_sols_zero, L_sols_unit, S),
			zdd_meet(AR_sols_unit, L_sols_unit, Sols_unit, S),
			zdd_join(Sols_zero, Sols_unit, Sols0, S)
		),
		prefix_pow_rev(H, Sols0, X, S)
	).
Helpers
		/*****************
		*     Helpers    *
		*****************/
% ?- #S, select_constant(0, X, Y, S).
% ?- #S, select_constant(1, X, Y, S).

% ?- #S, X<< (*[a]+1), psa(X, S), select_constant(X, C, Y, S), psa(Y, S).
% ?- #S, X<< (*[a]), select_constant(X, C, Y, S), psa(Y, S).

select_constant(X, 1, Z, S):- zdd_has_1(X, S), !,
	zdd_subtr(X, 1, Z, S).
select_constant(X, 0, X, _).

% ?- divide_at(a, [x,a, y], H, T).
% ?- divide_at(a, [x,y,a], H, T).
divide_at(A, L, H, T):- divide_at(A, L, [], H, T).
%
divide_at(A, [A|T], H, H, T):-!.
divide_at(A, [B|T], H, H0, T0):-
	divide_at(A, T, [B|H], H0, T0).

% ?- #S, prefix_pow_rev([b, a], 1, X, S), psa(X, S).
% ?- #S, prefix_pow_rev([b], 1, X, S), prefix_pow_rev([a], X, Y, S).
prefix_pow_rev([], X, X, _).
prefix_pow_rev([A|As], X, Y, S):-
	cofact(X0, t(A, X, X), S),
	prefix_pow_rev(As, X0, Y, S).

Thanks for clear explanation. GF2 is a boolean field !

Luckily, I think, I did exactly as your comment.

Also luckily, I did so, being aware of 1 on diagonal.

Before asking help, I have to debug for myself to make points clear.
There are so many things I am not sure about my codes.

Fixing several careless bugs in codes, I could run the query which counts all O(6, G2) by a system of equations. Performance gets better but not so much as I expected, which was a little bit disappointing.

Below is predicate gf2_solve, which is a main codes, which I like because it uses mainly set operation in ZDD union/intersection/subtraction/powerset, though the codes polishing is necessary.

% ?- #S, time(gf2_count_orth(6, C, S)).
%@ % 21,529,533,838 inferences, 2558.382 CPU in 2606.041 seconds (98% CPU, 8415293 Lips)
%@ S = ..,
%@ C = 23040.
Main codes gf2_solve.
gf2_solve(Ls, P, X, S):- memo(gf2_solve(P, Ls)-X, S),  % P>1.
	(	nonvar(X) -> true	%, write(.)	% no hits. Useless here ?
	;	cofact(P, t(A, L, R), S),
		divide_at(A, Ls, H, T),
		zdd_power(T, PowT, S),
		cofact(PowAT, t(A, PowT, PowT), S),
		(	R = 1 -> R0 = 0
		;	gf2_solve(T, R, R0, S)
		),
		cofact(AR_sols_zero, t(A, PowT, R0), S),
		(	L = 0 -> X0 = AR_sols_zero
		;	(	L = 1 -> L0 = PowT
			;	gf2_solve(T, L, L0, S)
			),
			L_sols_zero = L0,
			zdd_meet(AR_sols_zero, L_sols_zero, Sols_zero, S),
			zdd_subtr(PowAT, AR_sols_zero, AR_sols_unit, S),
			zdd_subtr(PowAT, L_sols_zero, L_sols_unit, S),
			zdd_meet(AR_sols_unit, L_sols_unit, Sols_unit, S),
			zdd_join(Sols_zero, Sols_unit, X0, S)
		),
		prefix_pow_rev(H, X0, X, S)
	).

It’s good to see this. I regret I never thought to use bits operations.

A little bit sounds strange. Dim(7,7) is seems not feasible for also my codes, as I watched the actitivity monitor of MacOS. By the way, I wonder how to compute for dim(n, n) n> 6. I guess they use some generating function for O(GF2, n).

There was a careless but fatal bug in the intersection operation zdd_meet/4. Fixing the one-letter bug, now my algorithm works as I expected.

% ?- time(count_orth(6, C)).
%@ % 29,124,413 inferences, 2.866 CPU in 2.910 seconds (98% CPU, 10162503 Lips)
%@ C = 23040.
% ?- call_with_time_limit(600, time(count_orth(7, C))).
%@ % 958,026,253 inferences, 147.126 CPU in 148.643 seconds (99% CPU, 6511622 Lips)
%@ C = 1451520.

I thought union/intersection/subtraction/multiplication operations
of my zdd library are stable and efficient, but intersection had a careless bug as I said.
I found the bug during testing a system of equations a + b = 0, a = 0, b = 0.
which surprised me.

My approach for counting orth(GF2, N) is extremely simple-minded.
First, read a given polynomial(s) P in ZDD form, then enumerate all GF2-value assignments which make value of P is 0, 1 respectively. This means polynomials and assginements coexist in ZDD at the
same time. In particular they use 0, 1 for both structures in a very simlar way, which confuses me often.
The other things are straightforward and easy routines.

Anyway, counting orth(GF2, N) of yours makes me more familiar with ZDD programming. Thanks.

Taking more than two hours cpu time and 80GB memory, it terminates with
cardinality for orth(GF2, 8) , which is more than I had expected.

iMac Intel 128GB, Ventura, SWI-Prolog (threaded, 64 bits, version 9.1.2-7-gcb02d296b),

% ?- time(count_orth(8, C)).
%@ % 77,269,642,066 inferences, 8239.463 CPU in 8300.403 seconds (99% CPU, 9377995 Lips)
%@ C = 185794560.

It would be nice if this result supports that even naive implementation of ZDD in SWI-Prolog
is useful for a certain class of computation.

EDIT.

Visting the site A003053 - OEIS, I have found the following reference

Jessie MacWilliams. “Orthogonal Matrices Over Finite Fields.”
The American Mathematical Monthly, vol. 76, no. 2, 1969, pp. 152–64. JSTOR,
https://doi.org/10.2307/2317262

in which a simple recursive formula for counting orth(GF(q), n). orth(GF(2), n) is a case with q = 2.

EDIT. 24/12/2022 Power of math !!

% ?- time(forall(between(2, 12, N), ( card_orth(2, N, V), writeln(orth(N)= V)))).
%@ orth(2)=2
%@ orth(3)=6
%@ orth(4)=48
%@ orth(5)=720
%@ orth(6)=23040
%@ orth(7)=1451520
%@ orth(8)=185794560
%@ orth(9)=47377612800
%@ orth(10)=24257337753600
%@ orth(11)=24815256521932800
%@ orth(12)=50821645356918374400
%@ % 215 inferences, 0.000 CPU in 0.000 seconds (78% CPU, 3467742 Lips)
%@ true.

% Using foldnum of pack pac.
card_orth(Q, N, V):-
	(	0 =:=  N mod 2 -> card_even_orth(Q, N, V)
	;	card_odd_orth(Q, N, V)
	).

card_orth_help(Q, T, I, U, V):- V is U* (Q^(2*T)-Q^(2*I)).

card_even_orth(Q, E, V):- T is E//2,
	foldnum(card_orth_help(Q, T), 1-(T-1), 1, U),
	V is Q^T* U.

card_odd_orth(Q, O, V):- E is O-1,
	card_even_orth(Q, E, U),
	V is (Q^E-1)*U.

Mine was a dead copy of the formula in the article without thought. Your code
is surprising me.

On my way of recalling my old codes on counting solutions of polynomials over GF(2), which is still unsuccessful, I have found a drastically simple way, and I wrote a prototype codes.

gf2_poly_card(Poly, C) counts the number of solutions
Poly = 1 in GF(2). Here is a simple test on the prototype. I am not sure whether good or bad, and it may have bugs. I am interested in its simplicity, which uses heavily join operations of ZDD with being aware with additional property x+x=0 from GF(2).

Unfortunately my interest is narrow, and no idea and plan to develop details further of the codes in the future. I really hope your currently developing method is successful.

% ?- N=25, findall(a(I), between(1, N, I), A),
%	time((zdd gf2_poly_card(+(A), C))).
%@ % 889,217,468 inferences, 65.727 CPU in 65.883 seconds (100% CPU, 13528918 Lips)
%@ N = 25,
%@ A = [a(1), a(2), a(3), a(4), a(5), a(6), a(7), a(8), a(...)|...],
%@ C = 16777216.
% ?- 16777216 =:= 2^24.
%@ true.

Since some point of time, I believe that GF(2) is
exactly the same as the Boolean ring of order 2. So your
comment sounds for me that GF(2) is not extremely useful … . This disappointed me a little, but I have to believe you who are rich of experience as for XOR. Anyway I will keep your comment in mind. Thanks.

My sat uses DNF as normal form. So the “nasty” test
rather sounds for me that my sat works efficiently as I expected.

BTW, my dnf/cnf assumes input form is ground.
Prolog variables are converted with numbervars and
atomic propositions are ground prolog terms. I never
think my sat is used as substitute for clpb.

Here is the nasty test in addition with counting.
It also works “efficiently”

% ?- between(10,20,N), nasty(N,X), time(sat(X=:=X)), 
%	sat_count(C), writeln(C), fail; true.
%@ % 274,225 inferences, 1.584 CPU in 2.891 seconds (55% CPU, 173102 Lips)
%@ 2048
%@ % 349,802 inferences, 0.299 CPU in 0.518 seconds (58% CPU, 1168211 Lips)
%@ 4096
%@ % 469,595 inferences, 0.322 CPU in 0.550 seconds (59% CPU, 1459013 Lips)
%@ 8192
%@ % 734,666 inferences, 0.369 CPU in 0.596 seconds (62% CPU, 1988836 Lips)
%@ 16384
%@ % 1,118,548 inferences, 0.461 CPU in 0.700 seconds (66% CPU, 2427178 Lips)
%@ 32768
%@ % 1,822,181 inferences, 0.648 CPU in 0.912 seconds (71% CPU, 2812898 Lips)
%@ 65536
%@ % 3,183,806 inferences, 1.007 CPU in 1.315 seconds (77% CPU, 3161699 Lips)
%@ 131072
%@ % 5,859,967 inferences, 1.732 CPU in 2.127 seconds (81% CPU, 3382993 Lips)
%@ 262144
%@ % 11,280,442 inferences, 2.713 CPU in 3.268 seconds (83% CPU, 4158314 Lips)
%@ 524288
%@ % 21,852,124 inferences, 5.585 CPU in 6.498 seconds (86% CPU, 3912755 Lips)
%@ 1048576
%@ % 42,888,461 inferences, 12.466 CPU in 15.622 seconds (80% CPU, 3440310 Lips)
%@ 2097152
%@ true.