Factored again

[Waldek Hebisch]

In my developement version I switched polynomial factorization
over prime fields to use new factorizer. That gives substantial
speedup. However, profiling showed that in my test case (which
was polynomial with 84 linear factors) abut 10% time was spent
in multiplying factors, of of which 8.4% went into comparing
exponents (POLYCAT-;smaller?)

One trouble is that comparizons of polynomials are expensive.
Second, multiplying factor-by-factor triggers quadratic
behaviour in Factored. For use in factorizers this can
be fixed by passing list of factors to 'makeFr'.
We probably should have specialized comparizon routine
in SMP.

It is not entirely clear to me what to do with Factored.
Ordering can be used to speed up computations only when
we have true order, otherwise we may get wrong results.
Additionally, ordering factors gives speedup only when
all or most factors is prime, for other factors we need
to use quadratic method anyway. My conclusion is that
we can not have really efficient operations with
generality offered by Factored, so probably should
not care too much about speed of Factored.

--
Waldek Hebisch

[Ralf Hemmecke]

On 14.07.23 16:48, Waldek Hebisch wrote:
> In my developement version I switched polynomial factorization
> over prime fields to use new factorizer. That gives substantial
> speedup. However, profiling showed that in my test case (which
> was polynomial with 84 linear factors) abut 10% time was spent
> in multiplying factors, of of which 8.4% went into comparing
> exponents (POLYCAT-;smaller?)

Just a side remark.
Isn't the implementation of POLYCAT-;smaller? a bit weird if it is used
in factorization? It basically, compares the leading coefficients.
So we have...

(1) -> S := SMP(Integer, Symbol)

(31) -> px := x::S

(31) x
Type: SparseMultivariatePolynomial(Integer,Symbol)
(32) -> py := y::S

(32) y
Type: SparseMultivariatePolynomial(Integer,Symbol)
(33) -> smaller?(-py^100,px)

(33) true

Certainly not a problem for the definition of smaller?, as it just wants
a linear order, but for factorization I would have expected that a
polynomial with a smaller degree counts as smaller.

Sorry, I haven't looked into the factorization code to check why such a
linear order makes sense. In fact, I wonder why the factorizer uses the
linear order by smaller? at all. Doesn´t Comparable say that it is just
an arbitrary but total order. I wouldn´t be surprised if I override
smaller? in SMP by some other crazy order that the factorization might
go wrong. Wild speculation without looking at the code, of course, but I
understood that it makes sense to introduce Comparable for a printing
order, but I am not sure about it's usefulness in other algorithms, that
want a few more properties than just totality in order to get some
complexity measure.

> It is not entirely clear to me what to do with Factored.
> Ordering can be used to speed up computations only when
> we have true order, otherwise we may get wrong results.

In what sense does the order influence correctness?

> Additionally, ordering factors gives speedup only when
> all or most factors is prime, for other factors we need
> to use quadratic method anyway. My conclusion is that
> we can not have really efficient operations with
> generality offered by Factored, so probably should
> not care too much about speed of Factored.

Oh, do you mean the sx := sort!(LispLessP, x) in makeFR $ Factored(R) ?
Why should one care about the order in the representation of Factored?

Ralf

[Waldek Hebisch]

On Fri, Jul 14, 2023 at 08:29:37PM +0200, Ralf Hemmecke wrote:
> On 14.07.23 16:48, Waldek Hebisch wrote:
> > In my developement version I switched polynomial factorization
> > over prime fields to use new factorizer. That gives substantial
> > speedup. However, profiling showed that in my test case (which
> > was polynomial with 84 linear factors) abut 10% time was spent
> > in multiplying factors, of of which 8.4% went into comparing
> > exponents (POLYCAT-;smaller?)
>
> Just a side remark.
> Isn't the implementation of POLYCAT-;smaller? a bit weird if it is used in
> factorization? It basically, compares the leading coefficients.
> So we have...
>
> (1) -> S := SMP(Integer, Symbol)
>
> (31) -> px := x::S
>
> (31) x
> Type: SparseMultivariatePolynomial(Integer,Symbol)
> (32) -> py := y::S
>
> (32) y
> Type: SparseMultivariatePolynomial(Integer,Symbol)
> (33) -> smaller?(-py^100,px)
>
> (33) true
>
> Certainly not a problem for the definition of smaller?, as it just wants a
> linear order, but for factorization I would have expected that a polynomial
> with a smaller degree counts as smaller.

Well, SparseMultivariatePolynomial uses lexicographic order. The
same order is used for Groebner bases. Current answer to your
wish for degree order is "use domain with degree order like
HomogeneousDistributedMultivariatePolynomial". I must say that
I am not fully satisfied with this answer, but lexicographic
order as _default_ order for SparseMultivariatePolynomial
makes a lot of sense. Including coefficients is necessary to
have ordered ring. Note that over fields our univariate factors
are monic. Over integers we normalize factors so that leading
coefficient is positive. So in both cases bigger degree will
win.

> Sorry, I haven't looked into the factorization code to check why such a
> linear order makes sense. In fact, I wonder why the factorizer uses the
> linear order by smaller? at all. Doesn´t Comparable say that it is just an
> arbitrary but total order. I wouldn´t be surprised if I override smaller? in
> SMP by some other crazy order that the factorization might go wrong. Wild
> speculation without looking at the code, of course, but I understood that it
> makes sense to introduce Comparable for a printing order, but I am not sure
> about it's usefulness in other algorithms, that want a few more properties
> than just totality in order to get some complexity measure.

Well, main reason for order is to get repeatable results when
printing. But if you have ordered lists you can find duplicates
using linear scan. Without extra structure like order or hash
function elimination of duplicates is quadratic. When all
factors are prime and unitCanonical holds (so factors are
uniquely represented) gcd is equivalent to finding duplicate
factors and adding exponents, multiplication is eqivalent
to merge of sorted lists with removal of duplicates.

Unfortunatly conditions above are quite restrictive, so in many
cases we need to fall back to quadratic version.

> > It is not entirely clear to me what to do with Factored.
> > Ordering can be used to speed up computations only when
> > we have true order, otherwise we may get wrong results.
>
> In what sense does the order influence correctness?

See above: one may miss fact that to factors are equvalent
and count them as separate factors giving clearly wrong gcd
and slightly wrong result of multiplication.

> > Additionally, ordering factors gives speedup only when
> > all or most factors is prime, for other factors we need
> > to use quadratic method anyway. My conclusion is that
> > we can not have really efficient operations with
> > generality offered by Factored, so probably should
> > not care too much about speed of Factored.
>
> Oh, do you mean the sx := sort!(LispLessP, x) in makeFR $ Factored(R) ?
> Why should one care about the order in the representation of Factored?

Well, we want repeateble results, so evan crazy order is better
than no order. Unfortunately, GGREATERP which is used as
fallback may give different results on different runs.

Currently our factorizers use (pseudo)random numbers, so
without extra ordering effort order of factors will vary.
Over resonable finite fields we have Comparable and can
order. But we also use random evaluation points for
multivariate polynomials, here we need some ordering
effort. Over sufficiently complicated domains our
ordering efforts will fail, so probably we should
concentrate on resonable order for simple domains.

--
Waldek Hebisch

[Ralf Hemmecke]

>> Certainly not a problem for the definition of smaller?, as it just
>> wants a linear order, but for factorization I would have expected
>> that a polynomial with a smaller degree counts as smaller.
>

> Well, SparseMultivariatePolynomial uses lexicographic order.No. That
> is not true.S := SMP(Integer, Symbol);px := x::S;py := y::S;
smaller?(-py^100,px) -- returns true
smaller?( py^100,px) -- returns false

For Gröbner bases lex order just works on the exponent vectors and
doesn't care about the coefficient ring at all. Certainly not a total
order suitable for Comparable. A simpler probably cheaper total order
for SMP would be to just compare degree in the main variable and only if
that is equal apply smaller? to the leading coefficient (a polynomial in
one variable less). That also ends with calling smaller?(r1,r2) on the
base ring R, but most probably recusion is not so deep because the
polynomials differ in degree in the main variable.

> Current answer to your wish for degree order is "use domain with
> degree order like HomogeneousDistributedMultivariatePolynomial".

No, no. I never wished for something. Just my expectations were not met.

> I must say that I am not fully satisfied with this answer, but
> lexicographic order as _default_ order for
> SparseMultivariatePolynomial makes a lot of sense.

Yes, of course. But what order are you talking about, the one given by
smaller? is maybe lexicographic, but not in the sense it is usually used
in Gröbner bases theory. POLYCAT-;smaller? implement an order that works
by virtually introducing 0-coefficients in order to bring both
polynomials to the same degree. Possible, but I have the feeling that it
is costly.

> Including coefficients is necessary to have ordered ring. Note that
> over fields our univariate factors are monic. Over integers we
> normalize factors so that leading coefficient is positive. So in
> both cases bigger degree will win.

I don't check now, but on this subset both implementations might give
the same result, but I still claim that the one that only recurses if
the degrees in the main variable are equal has less work to do.

> adding exponents, multiplication is eqivalent to merge of sorted
> lists with removal of duplicates.

Oh, now I see. You just perform multiplication in Factored, i.e., adding
exponents. No real multiplication of polynomials. I did not understand
this from your previous mails. OK, then ordering makes, of course, sense.

Ralf

[Waldek Hebisch]

On Fri, Jul 14, 2023 at 10:43:51PM +0200, Ralf Hemmecke wrote:
> > > Certainly not a problem for the definition of smaller?, as it just
> > > wants a linear order, but for factorization I would have expected
> > > that a polynomial with a smaller degree counts as smaller.
> >
> > Well, SparseMultivariatePolynomial uses lexicographic order.No. That is
> > not true.S := SMP(Integer, Symbol);px := x::S;py := y::S;
> smaller?(-py^100,px) -- returns true
> smaller?( py^100,px) -- returns false
>
> For Gröbner bases lex order just works on the exponent vectors and
> doesn't care about the coefficient ring at all.

Well, I would say that for Groebner bases order _is_ on exponent
vectors. But given order in base ring and order on exponents,
we have unique extention to polynomials satisfying axioms of
ordered ring such that monomial with coefficient 1 are ordered
via exponents.

> Certainly not a total
> order suitable for Comparable. A simpler probably cheaper total order
> for SMP would be to just compare degree in the main variable and only if
> that is equal apply smaller? to the leading coefficient (a polynomial in
> one variable less). That also ends with calling smaller?(r1,r2) on the
> base ring R, but most probably recusion is not so deep because the
> polynomials differ in degree in the main variable.

Well, for ordered sets currently Comparable gives the same order
as order from OrderedSet. Not stricly necessary, but simplifies
things because we can reuse implementations.

Concerning order in SparseMultivariatePolynomial, both 'degree' and
'leadingCoefficient' are expensive, but 'degree' is more expensice.
Namely, 'degree' produces list of pairs (symbol, exponent) (element of
IndexedExponents). To do this 'degree' must recurse. Similar
recursion is used for 'leadingCoefficient', but 'leadingCoefficient'
can just return its result. 'degree' must create appropriate
data structure. AFAICS there are several small and moderate
inefficiences in 'smaller?' that together lead to slowness.

> > Current answer to your wish for degree order is "use domain with degree
> > order like HomogeneousDistributedMultivariatePolynomial".
>
> No, no. I never wished for something. Just my expectations were not met.
>
> > I must say that I am not fully satisfied with this answer, but
> > lexicographic order as _default_ order for SparseMultivariatePolynomial
> > makes a lot of sense.
>
> Yes, of course. But what order are you talking about, the one given by
> smaller? is maybe lexicographic, but not in the sense it is usually used
> in Gröbner bases theory. POLYCAT-;smaller? implement an order that works
> by virtually introducing 0-coefficients in order to bring both
> polynomials to the same degree. Possible, but I have the feeling that it
> is costly.

That is just _one_ of several small things slowing 'smaller?'.

> > Including coefficients is necessary to have ordered ring. Note that
> > over fields our univariate factors are monic. Over integers we
> > normalize factors so that leading coefficient is positive. So in
> > both cases bigger degree will win.
>
> I don't check now, but on this subset both implementations might give the
> same result, but I still claim that the one that only recurses if the
> degrees in the main variable are equal has less work to do.

You missed significant part: 'smaller?' is in PolynomialCategory,
it only sees degree as element of IndexedExponents. In particular
'smaller?' has no idea that there is main variable. It works
essentially the same if you give it
HomogeneousDistributedMultivariatePolynomial. 'smaller?' simply
calls 'degree' to get representation of exponent vector, than
calls comparison functions to compare exponents. Then calls
'leadingCoefficient' to extract leading coefficient. Main cost
is due to genrality: all those things have to be SPADCALL-s.
Due to SPADCALL-s real work is probably less than overhead of
dispatch.

Specialized version for SparseMultivariatePolynomial could be
much faster, it still would have do work to extract leading
coefficient, but could do this in a loop using inlinable code
instead of recursion via SPADCALL-s.

--
Waldek Hebisch

[Ralf Hemmecke]

On 15.07.23 00:31, Waldek Hebisch wrote:
> Concerning order in SparseMultivariatePolynomial, both 'degree' and
> 'leadingCoefficient' are expensive, but 'degree' is more expensice.
> Namely, 'degree' produces list of pairs (symbol, exponent) (element of
> IndexedExponents). To do this 'degree' must recurse.

Hmmmm.... maybe I misinterpret something...

https://github.com/fricas/fricas/blob/master/src/algebra/polycat.spad#L647

This smaller? calls degree and degree is implemented in SMP.

The representation is recursive.

https://github.com/fricas/fricas/blob/master/src/algebra/multpoly.spad#L108

Now I would have expected the following. (Sorry this is without having
yet studied the implementation of SMP.)

Suppose I have two real polynomial not just elements of R, then they are
univariate polynomials in the main variable.
So the polynomial with the smaller main variable counts as smaller.
If they have the same main variable then check the degree in this
variable. And if the degrees are equal, recurse to the leading coefficient.

But hey, I see

https://github.com/fricas/fricas/blob/master/src/algebra/multpoly.spad#L647

degree p ==
p case R => 0
degree(leadingCoefficient(p.ts)) + monomial(degree(p.ts), p.v)

Oh, no! That computes the degree over all variables. It basically treat
the polynomial like in a distributed fashion. Sorry, I have not
recognized this before.

For an implementation of the smaller? from Compare in SMP, I would
definitely use the univariate recursive structure, i.e. consider the
degree function from D and eventually break ties with smaller? from R.

https://github.com/fricas/fricas/blob/master/src/algebra/multpoly.spad#L108C1-L108C1

D := SparseUnivariatePolynomial(%)
VPoly := Record(v : VarSet, ts : D)
Rep := Union(R, VPoly)

Certainly my previous mail must have been confusing. I was thinking of a
exactly this univariate degree function for smaller?.

And certainly that smaller? function must be implemented in SMP itself,
because it makes use of the particular representation.

Ralf

[Waldek Hebisch]

On Sat, Jul 15, 2023 at 01:12:08AM +0200, Ralf Hemmecke wrote:
>
> But hey, I see
>
> https://github.com/fricas/fricas/blob/master/src/algebra/multpoly.spad#L647
>
> degree p ==
> p case R => 0
> degree(leadingCoefficient(p.ts)) + monomial(degree(p.ts), p.v)
>
> Oh, no! That computes the degree over all variables. It basically treat the
> polynomial like in a distributed fashion. Sorry, I have not recognized this
> before.
>
> For an implementation of the smaller? from Compare in SMP, I would
> definitely use the univariate recursive structure, i.e. consider the degree
> function from D and eventually break ties with smaller? from R.

I have now version of 'smaller?' for SMP. I tried to get the same
results as older version. Good news is that when it works it
seem to be about 4 times faster than old version. Bad news is
that at least using sbcl it breaks in strange ways: printouts
show correct values, but it from time to time it dies with
Lisp errors. To complicate matter, it seem that after
succesfull call with given argument some time later it can
fail with the same arguments. It seems to work with ecl. ATM
it looks to me like a miscomplilation. But it is not clear
whom to blame, if our declaration/initialization is wrong
or it the error is on Lisp side.

If anybody wants to try this code and investigate problems
the patch is in attachement.

BTW: There is a chunk for ZMOD, this one seem to be working
fine.

--
Waldek Hebisch

[Grégory Vanuxem]

Hello Waldek,

I just cloned fricas.git
(configure.ac : 2023-07-17 20:19:21.130766889 +0200 on my drive).

And I guess your patch can not be used with recent Git version, am I right?

With SBCL-2.3.6 no problem without your patch. But if I (successfully)
apply your patch here is the output of 'make check':

=====================================
26 failing files:

bugs2007.output: 18
bugs2008.output: 8
bugs2009.output: 9
bugs2010.output: 7
bugs2012.output: 5
bugs2013.output: 4
bugs2014.output: 5
bugs2015.output: 4
bugs2016.output: 2
bugs2017.output: 24
bugs2018.output: 3
bugs2019.output: 1
bugs2020.output: 1
bugs2021.output: 4
bugs2022.output: 1
bugs2023.output: 4
derham.output: 21
distro.output: 1
ellip.output: 2
integ.output: 1072
limit.output: 63
lode.output: 13
mantepse.output: 24
polylift.output: 1
psgenfcn.output: 2
series3.output: 14
=====================================

That should be reproducible I think, so I did not attach the .output files.
If that's necessary I can put a tarball on my Google Drive and give a
link to it.
__
Greg

Le dim. 16 juil. 2023 à 22:40, Waldek Hebisch
<de...@fricas.math.uni.wroc.pl> a écrit :

> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/ZLRVr78748xllTd3%40fricas.math.uni.wroc.pl.

[Waldek Hebisch]

On Mon, Jul 17, 2023 at 09:12:41PM +0200, Grégory Vanuxem wrote:
> Hello Waldek,
>
> I just cloned fricas.git
> (configure.ac : 2023-07-17 20:19:21.130766889 +0200 on my drive).
>
> And I guess your patch can not be used with recent Git version, am I right?

Thank for looking at this. I am not sure what you mean. Patch is
with respect to recent git version, it can be applied to current
trunk. Concerning "can be used": it can be used to reproduce the
problem...

> With SBCL-2.3.6 no problem without your patch. But if I (successfully)
> apply your patch here is the output of 'make check':
>
> =====================================
> 26 failing files:

<snip>

>
> That should be reproducible I think,

I do not know if all sbcl versions produce exactly the same
set of failures. But there is a lot of failures, enough
to reproduce the problem. I do not think that knowing
exact set of failurs gives useful info.

--
Waldek Hebisch

[Grégory Vanuxem]

Le mar. 18 juil. 2023 à 00:44, Waldek Hebisch
<de...@fricas.math.uni.wroc.pl> a écrit :
>

> On Mon, Jul 17, 2023 at 09:12:41PM +0200, Grégory Vanuxem wrote:
> > Hello Waldek,
> >
> > I just cloned fricas.git
> > (configure.ac : 2023-07-17 20:19:21.130766889 +0200 on my drive).
> >
> > And I guess your patch can not be used with recent Git version, am I right?
>
> Thank for looking at this. I am not sure what you mean. Patch is
> with respect to recent git version, it can be applied to current
> trunk.

Ok.

In fact, I hesitated to ask for a simple test case where
FriCAS on sbcl with modified multpoly.spad dies sometimes with
Lisp errors. But since I had not upgraded FriCAS recently,
I decided to rebuild the latest git version and
was not expecting such a big _temporary_ regression.

For sure, I see this patch as "for testing purposes".

> Concerning "can be used": it can be used to reproduce the
> problem...

I understand now :)

> > With SBCL-2.3.6 no problem without your patch. But if I (successfully)
> > apply your patch here is the output of 'make check':
> >
> > =====================================
> > 26 failing files:
> <snip>
> >
> > That should be reproducible I think,
>
> I do not know if all sbcl versions produce exactly the same
> set of failures. But there is a lot of failures, enough
> to reproduce the problem. I do not think that knowing
> exact set of failurs gives useful info.

BTW, FriCAS built on Clozure CL does not exhibit this set of failures.
As ECL as you said.

> --
> Waldek Hebisch
>
> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/ZLXEJd99sX2bj4zI%40fricas.math.uni.wroc.pl.

[Grégory Vanuxem]

Waldek,

I think in your patch, from line 49
+ not(pv.v = qv.v) =>
+ smaller?(leadingCoefficient(p), 0$R)$R => 1
+ -- pv.v = qv.v

Should read
+ not(pv.v = qv.v) =>
+ smaller?(leadingCoefficient(p), 0$R)$R => 1
+ - 1
+ -- pv.v = qv.v

Or something like that. All tests pass with this modification.

In CCL interpreter:
(eql (the fixnum 1) (the fixnum nil))
issues a warning but returns nil whereas sbcl
aborts computation.
It could be that because from my point of view the
triage routine can return "nothing" in this piece of code
and the compiler should complain about that.

From my tests, distro.input (with )set break br)
the error comes from smaller? and complains in triage(p,q) = 1
Apparently triage(p,q) returns nil sometimes and the eql_SI macro
is used (src/lisp/primitives.lisp).

[Waldek Hebisch]

On Tue, Jul 18, 2023 at 04:36:38AM +0200, Grégory Vanuxem wrote:
> Waldek,
>
> I think in your patch, from line 49
> + not(pv.v = qv.v) =>
> + smaller?(leadingCoefficient(p), 0$R)$R => 1
> + -- pv.v = qv.v
>
> Should read
> + not(pv.v = qv.v) =>
> + smaller?(leadingCoefficient(p), 0$R)$R => 1
> + - 1
> + -- pv.v = qv.v
>
> Or something like that. All tests pass with this modification.

Yes, thanks. This should have triggered error in Spad compiler
(wrong return type), but has not...

Now one can test speed. Simple test could be:

pF := PF(nextPrime(2000000))
lp := [x - random()$pF for i in 1..84];
)set messages time on
sort(smaller?, lp);
[sort(smaller?, lp) for i in 1..10000];

Current version needs 2.55 sec for the last loop. With the patch
this drops to 0.69 sec.

--
Waldek Hebisch