Yes, the algorithm is very slow and should be improved, but I don't think it's buggy.

A workaround in this case is lreduce(lambda([q,r],zhegalkin_form(q and r)),args(foo)). rreduce is a bit faster, tree_reduce is a bit slower. In general, the order can be critical.

Of course, in general, the size of the Z form of an expression is exponential in the size of the input: Z(or(a[i],i,1,n)) has 2^n-1 terms, so both intermediate and final expression blow-up is a very real possibility.