Counting solutions for a polynomial over GF(2)

General
, ,
1

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.

2

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, …

3

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.

4

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).
6

Thanks for clear explanation. GF2 is a boolean field !

8

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

10

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

12

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.

13

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)
	).
14

The fastest I have, from the year 2018 (the problem itself entered
my circuite of experimentation already some decade before, it
was designed to challenge some CLP(X) ideas):

% ?- dim(6,6,X), time((basis(X,E), expr_tree(E,T), expr_vars(T, N, _))).
% % Up 2,150 ms, GC 166 ms, Thread Cpu 1,954 ms (Current 02/02/18 22:48:40)
% 23040

Its based on BDD, basically a BDD tree that is extended by some fields,
so that it carries at the same time a bitset indicating the variables in the
tree and an integer, indicating the solution count.

Drawback it needs around 8GB, and the next dim(7,7,X) cannot compute
anymore. Possibly the only place where I was using popcount/1 (called bitcount/1
here, have renamed it later). This is the BDD tree construction with the

typical rule ITE(X,T,T) = T and the additional book keeping. ZDD trees
are differently constructed. Its also not the first time that Kuniaki Mukai has
picked up my ortho problem, and it seems to be still unresolved somehow?

tree_make(A, A, _, R) :- !,
   R = A.
tree_make(A, B, X, node3(X,R,L,A,B)) :-
   expr_vars(A, S, U),
   expr_vars(B, T, V),
   W is U \/ V,
   N is bitcount(W /\ (\U)),
   M is bitcount(W /\ (\V)),
   R is 2^N*S+2^M*T,
   L is (1<<X) \/ W.
15

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).

16

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.
17

Now you gave me a nice christmas exercise, redoing the 2018 benchmark.
Its 2022 now, and I have new generation machines, one 11th Gen Intel and
one AMD Ryzen for testing. Wonder what the timing will be.

Your timing 2.910 seconds is not very far away from my 1,954 ms. But maybe
new ideas spring up. And what about a multi-core algorithm and/or some
symmetry breaking, to squeeze the lemon further.

Edit 23.12.2022
Did @hakank ever try the problem with Picat? Picat has some SAT Solver offerings.

18

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.

19

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.

20

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.
21

A solution without foldnum/4, note the simple alternating factors:

h(1, X) :- !, X = 1.
h(N, X) :-
    M is N-1,
    h(M, Y),
    X is (2^M-N rem 2)*Y.

Seems to give the same results:

?- between(2, 12, N), h(N, X), write(ortho(N)=X), nl, fail; true.
ortho(2) = 2 /* = 1*2 */
ortho(3) = 6 /* = 2*(4-1) */
ortho(4) = 48 /* = 6*8 */
ortho(5) = 720 /* = 48*(16-1) */
ortho(6) = 23040 /* = 720*32 */
ortho(7) = 1451520 /* = 23040*(64-1) */
ortho(8) = 185794560 /* = 1451520*128 */
ortho(9) = 47377612800 /* = 85794560*(256-1) */
ortho(10) = 24257337753600 /* = 47377612800*512 */
ortho(11) = 24815256521932800 /* = 24257337753600*(1024-1) */
ortho(12) = 50821645356918374400 /* = 24815256521932800*2048 */
22

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

23

(Nothing of importance)

24

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.