Fast Term Factorization via Paige-Tarjan: Bypassing ==/2 and compare/2

General
1

I am experimenting with an alternative approach to term factorization based on the Paige-Tarjan minimization method. From a theoretical standpoint, I doubt the current implementation of term_factorized/3—which relies on compare/2 and ==/2—has a strict theoretical basis for evaluating cyclic term equivalence, and it inherently risks structural traversal overhead.

To address this, I implemented the initial extraction phase in C using the Foreign Language Interface. This completely bypasses Prolog-level structural checks and ==/2. Instead, it uses PL_compare to enforce strict O(1) physical memory address identity. It safely flattens cyclic structures into a minimal coalgebra vector, which is then passed back to Prolog for Paige-Tarjan minimization.

While this is “hot tested” and not yet exhaustively verified, the performance on massive cyclic structures is promising. Here is a test using a deeply nested cyclic graph (hydra/3):

% ?- hydra(3, H, A), A=H.
hydra(0, A, A).
hydra(N, h(X, X), A):- N>0, N0 is N-1, hydra(N0, X, A).

pt_term_min_coa(Xs, Is, Mvec) :-
    c_extract_coa(Xs, Js, Tvec),
    pt_min_coa(Tvec, Mvec, Qmap),
    apply_quotient_map(Js, Qmap, Is).

On iMac intell 2020 macOS Tahoe 26.62
SWI-Prolog version 10.1.11 for x86_64-darwin

% ?- time((hydra(2000, _H, _A), _A=_H,  pt_term_min_coa([_H], Is, Mvec))).
%@ % 4,079 inferences, 0.011 CPU in 0.011 seconds (99% CPU, 376778 Lips)
%@ Is = [1],
%@ Mvec = u(h(1, 1)) .

% ?- time((hydra(10000, _H, _A), _A=_H,  pt_term_min_coa([_H], Is, Mvec))).
%@ % 20,078 inferences, 0.267 CPU in 0.268 seconds (100% CPU, 75065 Lips)
%@ Is = [1],
%@ Mvec = u(h(1, 1)) .

The C-level extraction isolates the physical pointers efficiently before applying the quotient map, completely avoiding infinite traversal or stack overflows without relying on ==/2.

2

Interesting indeed! But, PL_compare() is even more expensive than ==/2 as it walks both terms depth-first like ==/2, but when it finds a difference it must decide on the order rather than simply failing. There is also PL_same_compound() that compares two terms by address rather than content.

Anyway, a faster and more reliable term factorization is surely welcome, so please keep working on it.

3

Thanks for comment, Jan

Replacing ‘PL_compare’ with PL_same_compound, and adding some hashing,
It gets roughly 3 times faster, though not sure whether my hydra is suitable for benchmark test.
they always shrink drastically into one.

% ?- time((hydra(10000, _H, _A), _A=_H,  pt_term_coa([_H], Is, Mvec))).
%@ % 20,005 inferences, 0.106 CPU in 0.172 seconds (61% CPU, 189244 Lips)
%@ Is = [1],
%@ Mvec = u(h(1, 1)) .

% ?- time((generate_hydras(1000, _Xs), pt_term_coa(_Xs, Is, _Mvec))).
%@ % 1,006,003 inferences, 0.293 CPU in 0.401 seconds (73% CPU, 3436213 Lips)
%@ Is = [1, 1, 1, 1, 1, 1, 1, 1, 1|…] .

pt_term_coa(Xs, Is, Mvec) :-
    c_extract_coa(Xs, Js, Tvec),
    c_paige_tarjan(Tvec, Mvec, Qmap),
    apply_quotient_map(Js, Qmap, Is).

% ?- hydra(3, H, A), A=H.
hydra(0, A, A).
hydra(N, h(X, X), A):- N>0, N0 is N-1, hydra(N0, X, A).

generate_hydras(0, []) :- !.
generate_hydras(N, [H|Rest]) :-
    N > 0,
    hydra(N, H, A), A=H,
    N1 is N - 1,
    generate_hydras(N1, Rest).
4

If you play around with the current git version, there is a new library(random_terms) that provides random_term/2 that lets you create terms with a lot of parameters. In addition, xpce got a new class xdot in the Prolog library that does roughly the same as the Python xdot program: display GraphViz output with pan and zoom. You embed that in a lot of stuff, handle clicks on nodes and edges or add a menu to them. swipl-win has the menu GUI/GUI demo programs that provides a demo of the xdot viewer as well as a demo for random_terms/2 that shows the generated term graphically.

That should give you plenty of randomness to test your implementation :slight_smile:

5

How to use the random_term/2 ?
No hit by apropos and help despite of use_module(library(random_terms).

6

The library is not included in the documentation. It has PlDoc comments inside though. Basically it is just

  ?- random_term(Term, []).

The option list defines domains and weights that define the distribution of the shape of the term and the types of the leaf nodes. Suggestions for further are welcome. We may add predicates to generate atoms/strings with certain properties as well as numbers with certain properties and distributions.

I wonder whether we should move this from the core to the plunit package?

7

Thanks. ‘random_term/2’ looks like random/1 on numbers, and seems works for all types of terms, which is impressive.

% ?-  N = 100000, time(forall(between(1, N, _), 
%	(random_term(_T, []), pt_term_coa([_T], _Is, _Mvec)))).
%@ % 22,118,826 inferences, 0.942 CPU in 0.954 seconds (99% CPU, 23480557 Lips)
%@ N = 100000.
8

For the fun of it, here is an example from the GUI demo. Note that setting w_compound and the depth_decay higher gives you deeper terms (on average; still limited by depth)

image
image1504×935 84.9 KB

9

test(L, C) is a test predicate defined below. C is the cardinality of the set of distinct subterms of members of S, where S is a certain list of compound terms and L is the length of S.

This query repeats 10 times running test/2 with L = 10000 diplaying C for
each generated list S.

% ?- N = 10, Len=10000,  forall(between(1, N, _), (test(Len, C), writeln(C))).
%@ 16994
%@ 17098
%@ 17249
%@ 16763
%@ 17307
%@ 17841
%@ 16861
%@ 17209
%@ 16721
%@ 17015
%@ N = 10,
%@ Len = 10000.

For now I am not sure this is normal, and it seems difficult for me to judge,
perhaps also in the future.

%
test(Len, C):-random_terms(Len, Rs),
			pt_term_coa(Rs, Is, Mvec),
			count_compound(Mvec, C).
%
count_compound(V, C):- count_compound(V, 1, 0, C).
%
count_compound(V, I, N, N0):- arg(I, V, A), !,
	(	compound(A) -> U is N + 1
	;	U = N
	),
	J is I + 1,
	count_compound(V, J, U, N0).
count_compound(_, _, N, N).

%
random_terms(N, Rs):- length(Rs, N), 
	  maplist(random_term,  Rs).
%
random_term(R):- random_term(R, []).
%
compound_random_terms(N, Rs):- length(Rs, N),
	  maplist(compound_random_term,  Rs).
%
compound_random_term(R) :-
	(	random_term(R, []) -> true
	;	compound_random_term(R)
	).
10

Hi Jan,

On trying to see relation between term_factorized/3' and my pt_term_coa/3’,
I happened to find the following error, which is reproducible.
The code compoun_random_term/2 below is to drop 'compound arity zero` term
for simplicity.

% ?- N = 100, forall(between(1, N, _), (compound_random_term(R), term_factorized(R, _, L))).

%@ ERROR: Domain error: `compound_non_zero_arity' expected, found `
%@     [[ EXCEPTION while printing message ' expected, found `~p\''
%@        with arguments [h()]:
%@        raised: domain_error(compound_non_zero_arity,h())
%@     ]]
%@ 
%@ ERROR: In:
%@ ERROR:   [17] functor(
%@     [[ EXCEPTION while printing message '~p'
%@        with arguments [functor(h(),_354,_356)]:
%@        raised: domain_error(compound_non_zero_arity,h())
%@     ]]
%@ 
%@ ERROR:   [16] terms:mk_subst([... - _404,...],[_420=_422|_416],t(black('',_434,_436,''),black(...,...,_446,...))) at /Users/cantor/lib/swipl/library/terms.pl:309
%@ ERROR:   [13] forall('<garbage_collected>',pt_mincoa:(...,...)) at /Users/cantor/lib/swipl/boot/apply.pl:52
%@ ERROR:   [12] '<meta-call>'('<garbage_collected>') <foreign>
%@ ERROR:   [11] toplevel_call('<garbage_collected>') at /Users/cantor/lib/swipl/boot/toplevel.pl:1529
%@ ERROR: 
%@ ERROR: Note: some frames are missing due to last-call optimization.
%@ ERROR: Re-run your program in debug mode (:- debug.) to get more detail.

% ?-forall(between(1, 100000, _), compound_random_term(_)).
compound_random_term(R) :- random_term(R, []),
	compound(R),
	compound_name_arity(R, _, K),
	K > 0,
	!.
compound_random_term(R) :- compound_random_term(R).
11

For me, it just says:

105 ?- N = 100, forall(between(1, N, _), (compound_random_term(R), term_factorized(R, _, L))).
ERROR: Domain error: `compound_non_zero_arity' expected, found `g()'
ERROR: In:
ERROR:   [21] functor(g(),_5988,_5990)
ERROR:   [20] terms:insert_vars(g(),_6024,t(black('',_6040,_6042,''),black(...,.

Maybe you defined a portray/1 rule that does not handle h()?

P.s. If you want random terms without f(), use [p_zero_arity(0)] in the option list. That sets the probability to generate zero-arity compounds.

12

p_zero_arity(0) option is exactly what I need. Thanks. Compound arity zero object
should be treated as atomic (not compound) in state minimization context.

As the output of pt_term_coa/3 is a minimized deterministic automata,
there should be clear relation between them, which is what I am trying to show for now.

13

Due to the option ‘p_zero_arity(0)’, it is easy to find structure shared nodes
for given lists of terms generated by random_term. The following sample query
repeats N=10 times to generate a list of length Len = 10000 of compound terms
with ‘random_term’, and find equivalence classes, i.e., shared node
in the minimized deterministic automaton. In the output, for example, [29685,29312,10857],
means the nodes 29685,29312,10857 are all shared nodes which came from
given random terms. By first glance on the output, they seems shorter than expected,
So far every test so far results the similar, though not enough. I could not see reasons.

% ?- N = 10, Len = 10000,
%	time(forall(between(1, N, _), (find_shared_nodes(Len, Cs), writeln(Cs)))).
%@ []
%@ [29685,29312,10857]
%@ [17337,14004,11004]
%@ []
%@ [28419,20988,13619]
%@ [23141,19414]
%@ [28937,15642]
%@ [25465,12925,9392,9246]
%@ []
%@ [23378]
%@ % 26,925,581 inferences, 3.385 CPU in 3.462 seconds (98% CPU, 7954828 Lips)
%@ N = 10,
%@ Len = 10000.

In mathemacs, as it simply finds non-singleton equivance classes of the quotient set
induced by a surjection, there should be short codes. It is pity that I could not
find clear relationship between term_factorized/3 and the Paige-Tarjan menthod, though
there should be clear way to bridge between the two.

find_shared_nodes(Len, Cs):-compound_random_terms(Len, Rs),
			   pt_term_coa(Rs, _, Mvec, Qmap),
			   non_singleton_equiv_classes(Qmap, Mvec, Cs).

% ?-forall(between(1, 100000, _), compound_random_term(_)).
compound_random_term(R) :- random_term(R, [p_zero_arity(0)]).
% ?-  compound_random_terms(1000, _Rs).
compound_random_terms(N, Rs):- length(Rs, N),
	  maplist(compound_random_term,  Rs).
%
random_term(R):- random_term(R, []).
%
random_terms(N, Rs):- length(Rs, N),  maplist(random_term,  Rs).
%
rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
%
rel_to_fun([], X, X).
rel_to_fun([A-B|L], [A-U|V], R):-!,
	rel_to_fun(L, [A-[B|U]|V], R).
rel_to_fun([A-B|L], U, R):-!,
	rel_to_fun(L, [A-[B]|U], R).
%
non_singleton_equiv_classes(Qmap, Mvec, Cs):- Qmap=..[_|Js],
	functor(Qmap, _, N),
	zip(Js, Is, R),
	numlist(1, N, Is),
	sort(R, R0),
	rel_to_fun(R0, Fun),
	findall(X, (member(X-V, Fun),
				V=[_,_|_],
				arg(X, Mvec, A),
				compound(A)
			   ),
			Cs).

pt_term_coa(Xs, Is, Mvec) :- pt_term_coa(Xs, Is, Mvec, _).
%
pt_term_coa(Xs, Is, Mvec, Qmap) :-
    user:extract_coa(Xs, Js, Tvec),
    user:paige_tarjan(Tvec, Mvec, Qmap),
    apply_quotient_map(Js, Qmap, Is).

%
zip([], [], []).
zip([X|Xs], [Y|Ys], [X-Y|XYs]):- zip(Xs, Ys, XYs).
%
14

I dropped p_zero_arity(0) option from random_term/2 with additional portray action.
and modified pt_term_coa in foreign C inteferface works well as expected.

From curiosity, I compared term_factorized/3 and the pt_term_coa/4 on time/1
for a massive hydra. The result is much more than I expected so that
I have to suspect some bug there in my pt_term_coa. I do not think hydra is not approapriate
benchmark at all.

% ?- N = 100000, hydra(N, _H, _R), _H = _R,
%	time(term_factorized(_H, _A, _L)),
%	time(pt_term_coa([_H], _I, _M, _Q)).
%@ % 67 inferences, 74.050 CPU in 74.364 seconds (100% CPU, 1 Lips)
%@ % 4 inferences, 0.007 CPU in 0.009 seconds (74% CPU, 601 Lips)
%@ N = 100000 .
15

I’d first setup a test using the random term generator and comparing the answers of the two implementations.

16

Hi Jan,

Using random_term/2 and variant/2, I have run the pt_term_coa/3 100000 times
to find a possible bug, but found no error, which of course does not mean it is correct, though.

You are surely busy for many other tasks of higher priority. I am just posting this since
I have managed to find a prolog codes ‘expand_coa_node/3’ below which converts a coa as a system of equations of flat terms on right hand side to terms.

I think the result of this benchmark may not be bad and enough to push me to polish it further.

% ?- time(random_coa_bench(100000)).
%@ % 27,228,220 inferences, 2.782 CPU in 2.785 seconds (100% CPU, 9787536 Lips)
%@ true.

random_coa_bench(N):- forall(between(1, N, _),
	   (	random_term(R, []),
			pt_term_coa([R], [I], Coa),
			expand_coa_node(I, Coa, T),
			variant(R,  T))).
%
expand_coa_node(I, V, T):- functor(V, F, N),
	functor(U, F, N),
	coa_to_zip_by_equality(1, V, U, Z),
	maplist(call, Z),
	arg(I, U, T).
%
coa_to_zip_by_equality(I, V, U, [X = Tc|Z]):- arg(I, V, T), !,
	arg(I, U, X),
	map_node_args(U, T, Tc),
	J is I + 1,
	coa_to_zip_by_equality(J, V, U, Z).
coa_to_zip_by_equality(_, _, _, []).
%
map_node_args(U, T, Tc):-
	(	compound(T), compound_name_arity(T, F, K),  K > 0
	->  functor(Tc, F, K),
		map_node_args(1, U, T, Tc)
	;	Tc = T
	).
%
map_node_args(I, U, T, Tc):-arg(I, T, J), !,
	arg(J, U, X),
	arg(I, Tc, X),
	K is I + 1,
	map_node_args(K, U, T, Tc).
map_node_args(_, _, _, _).
17

From what you write, this can provide a better term_factorized/3, If I undertand you correctly, this should be faster and, I guess, also save against the ordering issues with cyclic terms that can make term_factorized/3 misbehave if I recall well. Please just package it up and let us play with it!