zdd-gf2.pl

    1:- module(gf2, [gf2_poly/3, gf2_poly_list/5,
    2				gf2_poly_card/3, gf2_card/4,
    3				gf2_poly_solve/3, gf2_poly_solve_list/3,
    4				gf2_solve/4, gf2_solve_zero/4,
    5				gf2_join/4, gf2_merge/4, gf2_join_list/4, gf2_merge_list/4,
    6				gf2_subst/5, gf2_insert/4,
    7				collect_atoms/2, collect_atoms_list/2, divide_at/4
    8			   ]).    9
   10:- use_module(library(apply)).   11:- use_module(library(apply_macros)).   12:- use_module(zdd('zdd-array')).   13:- use_module(util(math)).   14:- use_module(util(meta2)).   15:- use_module(pac(basic)).   16:- use_module(pac(meta)).   17:- use_module(util(misc)).   18:- use_module(zdd(zdd)).   19:- use_module(pac(op)).   20
   21		/*******************
   22		*     Benchmark    *
   23		*******************/
   24% ?- zdd X<<pow(:numlist(1, 100)).
   25%@ X = 101.
   26
   27% Power Polynomial for benchmark
   28% P1 = 1 + x1
   29% P(n+1) = Pn + x(n+1)*P(n)  (n>0)
   30% ?- power_poly(1, X).
   31% ?- power_poly(2, X).
   32power_poly(1, 1 + x(1)):-!.
   33power_poly(N, P+x(N)*P):- N0 is N-1, power_poly(N0, P).
   34
   35% ?- power_poly(3, P), (zdd gf2_poly(P, Q
   36
   37% ?- power_poly(1, _P), (zdd gf2_poly_solve_card(_P, X)).
   38% ?- power_poly(2, _P), (zdd gf2_poly_solve_card(_P, X)).
   39% ?- power_poly(3, _P), (zdd gf2_poly_solve_card(_P, X)).
   40% ?- power_poly(10, _P), (zdd gf2_poly_solve_card(_P, X)).
   41% ?- power_poly(20, _P), (zdd gf2_poly_solve_card(_P, X)).
   42% ?- power_poly(30, _P), (zdd gf2_poly_solve_card(_P, X)). % running for ever.
   43
   44% Read the term by zdd commands.
   45% ?- numlist(1, 1000, Ns),
   46%	time(((zdd zdd_power(Ns, P), gf2_poly_solve_card(P, C)))).
   47%@ % 45,060 inferences, 0.004 CPU in 0.004 seconds (86% CPU, 11780392 Lips)
   48%@ Ns = [1, 2, 3, 4, 5, 6, 7, 8, 9|...],
   49%@ P = 1001,
   50%@ C = 1.
   51
   52		/***************
   53		*     Notes    *
   54		***************/
   55
   56% For a, b in GF2, 1 and 2 hold.
   57% 1.  a  = 0  iff  1 + a = 1.
   58% 2.  ab = 1  iff  a = 1 and b = 1.
   59% 3.  Orthogonality of matrix A = [ aij ] ( n=< i, j =< n )
   60% is represented as a system of equations Eij = 0,
   61% where Eij is a polynomial with indeterminates aij obtained
   62% from the familiar multiplication of matricies.
   63
   64% It follows from 1, 2, 3 that the set Orth(n, GF2)
   65% and the intersection of sol(Eij=0) have the same cardinality.
   66% where Sol(E) is the set of solutions of E = 0.
   67
   68		/*****************
   69		*     Helpers    *
   70		*****************/

divide_at(+A, +L, -H, -T) is det: A : ZDD atom. L : sorted list. H,T are unified with lists such that append of reverse(H), [A], and T equals to the given L.

   78% ?- divide_at(a, [x, a, y], H, T).
   79% ?- divide_at(a, [x, y, a], H, T).
   80divide_at(A, L, H, T):- divide_at(A, L, [], H, T).
   81%
   82divide_at(A, [A|T], H, H, T):-!.
   83divide_at(A, [B|T], H, H0, T0):-
   84	divide_at(A, T, [B|H], H0, T0).
   85
   86		/**************************
   87		*     gf2_system_solver   *
   88		**************************/

C is unified with the number of solution of P = 1, where P is a polynomial over GF(2).

   94% ?- zdd gf2_poly_card(a, C).
   95% ?- zdd gf2_poly_card(a+b, C).
   96% ?- zdd gf2_poly_card(a+b+c, C).
   97% ?- zdd gf2_poly_card(a+b+c+d, C).
   98% ?- zdd gf2_poly_card(a, C).
   99% ?- zdd gf2_poly_card(a*b, C).
  100% ?- zdd gf2_poly_card(a*b*c, C).
  101% ?- zdd gf2_poly_card(a*b*c*d, C).
  102% ?- N=25, findall(a(I), between(1, N, I), A),
  103%	time((zdd gf2_poly_card(+(A), C))).
  104
  105gf2_poly_card(P, C, S):- collect_atoms(P, Vs),
  106	gf2_poly(P, X, S),
  107	gf2_card(Vs, X, C, S).
  108
  109% ?- zdd gf2_poly_solve_card(a, C).
  110% ?- zdd gf2_poly_solve_card(a+b, C).
  111% ?- zdd gf2_poly_solve_card(a+b+c, C).
  112% ?- zdd gf2_poly_solve_card(a+b+c+d, C).
  113% ?- zdd gf2_poly_solve_card(a, C).
  114% ?- zdd gf2_poly_solve_card(a*b, C).
  115% ?- zdd gf2_poly_solve_card(a*b*c, C).
  116% ?- zdd gf2_poly_solve_card(a*b*c*d, C).
  117%
  118gf2_poly_solve_card(P, C, S):- collect_atoms(P, Vs),
  119							gf2_poly(P, X, S),
  120							gf2_solve(Vs, X, Sols, S),
  121							card(Sols, C, S).
  122%
  123gf2_card(_,  0, 0, _):-!.
  124gf2_card(Vs, 1, N, _):-!, length(Vs, K), N is 2^K.
  125gf2_card(Vs, X, N, S):- memo(gf2_card(X, Vs)-N, S),
  126	(	nonvar(N) -> true
  127	; 	cofact(X, t(A, L, R), S),
  128		divide_at(A, Vs, H, T),
  129		gf2_card(T, L, I, S),
  130		gf2_join(L, R, U, S),
  131		gf2_card(T, U, J, S),
  132		N0 is I + J,
  133		length(H, K),
  134		N is N0*(2^K)
  135	).
  136
  137% ?- fetch_state(S), gf2_poly_solve(a+b+c, X, S), card(X, C, S).
  138% ?- zdd gf2_poly_solve(a+b*c, X), card(X, C).
  139% ?- findall(a(I), between(0, 1000, I), V),
  140%	(zdd gf2_poly_solve(+V, X), card(X, C)).
  141% ?- findall(a(I), between(0, 1000, I), V),
  142%	(zdd gf2_poly_solve(*(V), X), card(X, C)).
  143% ?- findall(a(I), between(1, 3, I), V),
  144%	(zdd gf2_poly_solve(+(V), X), card(X, C)).
  145% ?- zdd gf2_poly_solve(+[a,b, a, c,d], X), card(X, C).
  146% ?- zdd gf2_poly_solve(+[a,b, a, b, c, c, 1], X), card(X, C).
  147
  148gf2_poly_solve(X, Y, S):- collect_atoms(X, Vs),
  149					gf2_poly(X, X0, S),
  150					gf2_solve(Vs, X0, Y, S).
  151
  152
  153% ?- fetch_state(S), gf2_poly_solve_list([a], Y, S), card(Y, C, S).
  154% ?- fetch_state(S), gf2_poly_solve_list([a+b+c, a+b], Y, S), card(Y, C, S).
  155gf2_poly_solve_list(X, Y, S):- collect_atoms_list(X, A),
  156					gf2_poly_list(X, X0, S),
  157					zdd_power(A, Y0, S),
  158					gf2_solve_list(A, X0, Y0, Y, S).
  159%
  160gf2_solve_list(_, [], P, P, _):-!.
  161gf2_solve_list(Vs, [X|Xs], P, Q, S):-
  162	gf2_solve(Vs, X, P0, S),
  163	zdd_meet(P0, P, Q0, S),
  164	gf2_solve_list(Vs, Xs, Q0, Q, S).
  165
  166%
  167gf2_solve(_,  0, 0, _):-!.
  168gf2_solve(Vs, 1, P, S):-!, zdd_power(Vs, P, S).
  169gf2_solve(Vs, X, P, S):- memo(gf2_unit(X, Vs)-P, S),
  170	(	nonvar(P) -> true
  171	; 	cofact(X, t(A, L, R), S),
  172		divide_at(A, Vs, H, T),
  173		gf2_solve(T, R, R0, S),
  174		zdd_insert(A, R0, R1, S),
  175		(	L = 0 -> P0 = R1
  176		;	gf2_solve(T, L, L0, S),
  177			cofact(L1, t(A, L0, L0), S),
  178			zdd_xor(L1, R1, P0, S)
  179		),
  180		zdd_power_rev(H, P0, P, S)
  181	).
  182
  183% Without memo,  it goes easy into looping for "power polynomial"
  184% gf2_solve_without_memo(_,  0, 0, _):-!.
  185% gf2_solve_without_memo(Vs, 1, P, S):-!, zdd_power(Vs, P, S).
  186% gf2_solve_without_memo(Vs, X, P, S):-
  187% 	 	cofact(X, t(A, L, R), S),
  188% 		divide_at(A, Vs, H, T),
  189% 		gf2_solve_without_memo(T, R, R0, S),
  190% 		zdd_insert(A, R0, R1, S),
  191% 		(	L = 0 -> P0 = R1
  192% 		;	gf2_solve_without_memo(T, L, L0, S),
  193% 			cofact(L1, t(A, L0, L0), S),
  194% 			zdd_xor(L1, R1, P0, S)
  195% 		),
  196% 		zdd_power_rev(H, P0, P, S).

gf2_solve_zero(+Vs, +X, -P, +S) is det: P is unified with the set of assignments which makes X = 0.

  202gf2_solve_zero(Vs, X, P, S):-
  203	gf2_solve(Vs, X, U, S),
  204	zdd_power(Vs, W, S),
  205	zdd_subtr(W, U, P, S).
  206
  207%
  208zdd_power_rev([], P, P, _).
  209zdd_power_rev([X|Xs], P, Q, S):- cofact(P0, t(X, P, P), S),
  210						 zdd_power_rev(Xs, P0, Q, S).
  211
  212		/***********************************
  213		*     Read polynomial over GF(2).  *
  214		***********************************/
  215
  216% ?- collect_atoms(*[a,b,a], A).
  217collect_atoms(X, A):- collect_atoms(X, B, []), sort(B, A).
  218% ?- collect_atoms_list([a,b,a], A).
  219collect_atoms_list(X, A):- collect_atoms_list(X, B, []), sort(B, A).
  220
  221%
  222collect_atoms(I, A, A)	:- integer(I), !.
  223collect_atoms(-(X), A, B):-!, collect_atoms(X, A, B).
  224collect_atoms(P-Q, A, B):-!, collect_atoms(P+(-Q), A, B).
  225collect_atoms(X * Y, A, B):-!,
  226	collect_atoms(X, A, C),
  227	collect_atoms(Y, C, B).
  228collect_atoms(X + Y, A, B):-!,
  229	collect_atoms(X, A, C),
  230	collect_atoms(Y, C, B).
  231collect_atoms(+(V), A, B):-!, collect_atoms_list(V, A, B).
  232collect_atoms(*(V), A, B):-!, collect_atoms_list(V, A, B).
  233collect_atoms(X,[X|A], A).
  234
  235%
  236collect_atoms_list([], A, A):-!.
  237collect_atoms_list([X|Xs], A, B):-
  238	collect_atoms(X, A, C),
  239	collect_atoms_list(Xs, C, B).
  240
  241%
  242gf2_poly_list([], [], _).
  243gf2_poly_list([P|Ps], [Q|Qs], S):-
  244	gf2_poly(P, Q, S),
  245	gf2_poly_list(Ps, Qs, S).

gf2_poly(+E, -V, +S) is det: E: polynomial over GF2. V is unified with E that is read as a ZDD. GF2 rules applied: e.g. a + a = 0. a*(b+c)=a*b+a*c.

  252% ?- zdd gf2_poly(a, X), psa(X).
  253% ?- zdd gf2_poly(a+a+a+b, X), psa(X).
  254% ?- zdd gf2_poly(-(a+a+a+b), X), psa(X).
  255% ?- zdd gf2_poly(1+ -2*a+b+1, X), psa(X).
  256% ?- zdd gf2_poly(a*a, X), psa(X).
  257% ?- zdd gf2_poly((a+b+c)*(c+b+a), X), psa(X).
  258% ?- zdd gf2_poly(*[+[a]], X), psa(X).
  259% ?- zdd gf2_poly(*[+[a,b], +[c,d]], X), psa(X).
  260% ?- zdd gf2_poly(+[a, b] +  *[b, c, d], X), psa(X).
  261% ?- zdd gf2_poly(+[a,b], X), psa(X).
  262% ?- zdd gf2_poly(*[a,b], X), psa(X).
  263% ?- zdd gf2_poly(a*1*(1*a + 1*a), X), psa(X).
  264%
  265gf2_poly(I, M,  _)	:- integer(I), !, M is I mod 2.
  266gf2_poly(-(X), Y, S):-!, gf2_poly(X, Y, S).
  267gf2_poly(P-Q, Y, S):-!, gf2_poly(P+(-Q), Y, S).
  268gf2_poly(X * X, Y, S):-!, gf2_poly(X, Y, S).
  269gf2_poly(X * Y, Z, S):-!, % memo(gf2_poly(X*Y)-Z, S),
  270		gf2_poly(X, U, S),
  271		gf2_poly(Y, V, S),
  272		gf2_merge(U, V, Z, S).
  273gf2_poly(X + X, 0, _):-!.
  274gf2_poly(X + Y, Z, S):-!, % memo(gf2_poly(X+Y)-Z, S),
  275		gf2_poly(X, U, S),
  276		gf2_poly(Y, V, S),
  277		gf2_join(U, V, Z, S).
  278%
  279gf2_poly(+(V), Y, S):-!, gf2_poly_list(+, V, 0, Y, S).
  280gf2_poly(*(V), Y, S):-!, gf2_poly_list(*, V, 1, Y, S).
  281gf2_poly(X, Y, S)	:- zdd_singleton(X, Y, S).
  282%
  283gf2_poly_list(_, [], Y, Y, _):-!.
  284gf2_poly_list(F, [X|Xs], Y0, Y, S):-
  285	gf2_poly(X, A, S),
  286	(	F = (*) -> gf2_merge(A, Y0, Y1, S)
  287	;	gf2_join(A, Y0, Y1, S)
  288	),
  289	gf2_poly_list(F, Xs, Y1, Y, S).
  290
  291		/***********************************************
  292		*   Operations on polynomial over GF(2)    *
  293		***********************************************/

gf2_insert(+X, +Y, -Z, +S) is det: Insert a literal X into each set of literals in a ZDD Y, and unify Z with the result. For short, Z = X*Y.

  300gf2_insert(_, 0, 0, _):-!.
  301gf2_insert(A, 1, J, S):-!, zdd_singleton(A, J, S).
  302gf2_insert(A, I, J, S):- memo(gf2_insert(I, A)-J, S),
  303	(	nonvar(J)	-> true
  304	;	cofact(I, t(B, L, R), S),
  305		zdd_compare(C, A, B, S),
  306		( 	C = (=)	 ->  gf2_insert(A, L, L0, S),
  307						 zdd_cons(R0, A, R, S),
  308						 gf2_join(L0, R0, J, S)
  309		;	C = (<)  ->  zdd_cons(J, A, I, S)
  310		;   C = (>)  ->  gf2_insert(A, R, R0, S),
  311						 gf2_insert(A, L, L0, S),
  312						 cofact(J, t(B, L0, R0), S)
  313		)
  314	).

gf2_join(+X, +Y, -Z) is det: Unify Z with GF(2)-polynomial addition X +Y. Z = X + Y.

  320% ?- zdd X<< {[a]}, gf2_join(X, X, Z), psa(Z).
  321% ?- zdd X<< {[1]}, gf2_join(X, X, Z), psa(Z).
  322% ?- zdd X<<{[1,a]}, gf2_join(X, X, Z), psa(Z).
  323% ?- zdd X<< ({[b]}+{[a]}), gf2_join(X, X, Z), psa(Z).
  324% ?- zdd X<< {[1]} + {[b]}+ {[a]}, zdd_merge(X, X, Z), psa(Z).
  325% ?- zdd X<< ({[b]} + {[a]}), psa(X).
  326% ?- zdd X<< {[1], [b], [a]}, zdd_merge(X, X, Z), psa(X), psa(Z).

gf2_join(+X, +Y, -Z, +S) is det: X, Y: zdd of GF(2)-polyminals. Z is unified with GF(2)-polynomial addition X + Y.

  332gf2_join(X, X, 0, _):-!.
  333gf2_join(0, X, X, _):-!.
  334gf2_join(X, 0, X, _):-!.
  335gf2_join(1, X, Y, S):-!, gf2_join_1(X, Y, S).
  336gf2_join(X, 1, Y, S):-!, gf2_join_1(X, Y, S).
  337gf2_join(X, Y, Z, S):-
  338	(	X@<Y -> memo(gf2_join(X,Y)-Z, S)
  339	;	memo(gf2_join(Y, X)-Z, S)
  340	),
  341	(	nonvar(Z) -> true %, write(.)
  342	; 	cofact(X, t(A, L, R), S),
  343		cofact(Y, t(A0, L0, R0), S),
  344		zdd_compare(C, A, A0, S),
  345		(	C = (<) -> 	gf2_join(L, Y, U, S),
  346						cofact(Z, t(A, U, R), S)
  347		;	C = (=)	->	gf2_join(R, R0, U, S),
  348						gf2_join(L, L0, V, S),
  349						cofact(Z, t(A, V, U), S)
  350		; 	gf2_join(L0, X, U, S),
  351			cofact(Z, t(A0, U, R0), S)
  352		)
  353	).
  354
  355% Switch the left-most leaf (0 or 1).
  356gf2_join_1(0, 1, _):-!.
  357gf2_join_1(1, 0, _):-!.
  358gf2_join_1(X, Y, S):-!,	cofact(X, t(A, L, R), S),
  359	gf2_join_1(L, L0, S),
  360	cofact(Y, t(A, L0, R), S).

gf2_merge(+X, +Y, -Z, +S) is det: X, Y: GF(2)-polyniomials in zdd. Z is unified with GF(2)-polyniomial addition X * Y. Z = X*Y.

  367% ?- zdd X<< *[a,b], Y<< *[c,d], gf2_merge(X, Y, Z), psa(Z).
  368% ?- zdd X<< *[a,b], Y<< *[b,d], gf2_merge(X, Y, Z), psa(Z).
  369% This is faster. Strange.
  370gf2_merge(X, X, X, _):-!.
  371gf2_merge(0, _, 0, _):-!.
  372gf2_merge(_, 0, 0, _):-!.
  373gf2_merge(1, X, X, _):-!.
  374gf2_merge(X, 1, X, _):-!.
  375gf2_merge(X, Y, Z, S):-
  376	(	X@<Y -> memo(gf2_merge(X,Y)-Z, S)
  377	;	memo(gf2_merge(Y,X)-Z, S)
  378	),
  379	(	nonvar(Z) -> true
  380	;	cofact(X, t(A, L, R), S),
  381		gf2_merge(L, Y, L0, S),
  382		gf2_merge(R, Y, R1, S),
  383		gf2_insert(A, R1, R0, S),
  384		gf2_join(L0, R0, Z, S)
  385	).
  386%
  387gf2_merge_list([], X, X, _).
  388gf2_merge_list([A|As], X, Z, S):-
  389	gf2_merge(A, X, Y, S),
  390	gf2_merge_list(As, Y, Z, S).
  391%
  392gf2_join_list([], X, X, _).
  393gf2_join_list([A|As], X, Y, S):-
  394	gf2_join(A, X, Z, S),
  395	gf2_join_list(As, Y, Z, S).
  396
  397% ?- zdd U<<{[a]}, P<<pow([x,y]), gf2_subst(x, U, P, Q), psa(Q).
  398% ?- zdd U<<{[a], [b]}, P<<pow([x,y]), gf2_subst(x, U, P, Q), psa(Q).
  399% ?- zdd U<<{[a, b]}, P<<pow([a, b, c, x,y]), gf2_subst(x, U, P, Q), psa(Q).
  400gf2_subst(_, _, X, X, _):- X<2, !.
  401gf2_subst(X, U, P, Q, S):- memo(gf2_subst(X, U, P)-Q, S),
  402	(	nonvar(Q) -> true	%, write(.) % So so hits.
  403	;	cofact(P, t(Y, L, R), S),
  404		gf2_subst(X, U, L, L0, S),
  405		zdd_compare(C, X, Y, S),
  406		(	C = (=) 	-> 	gf2_merge(U, R, R0, S)
  407		;	C = (<)		-> 	zdd_cons(R0, Y, R, S)
  408		;	gf2_subst(X, U, R, R1, S),
  409			gf2_insert(Y, R1, R0, S)
  410		),
  411		gf2_join(L0, R0, Q, S)
  412	)