bugs connected to PolynomialFactorizationExplicit

[Ralf Hemmecke]

There is seemingly a problem with the implementation of some functions
from PolynomialFactorizationImplizit in Integer as well as in
FractionInteger.

In fact, I can only guess what the function "factorSquareFreePolynomial"
should actually return.

http://fricas.github.io/api/PolynomialFactorizationExplicit.html#l-polynomial-factorization-explicit-factor-square-free-polynomial

factorSquareFreePolynomial: SparseUnivariatePolynomial % -> Factored
SparseUnivariatePolynomial %
factorSquareFreePolynomial(p) factors the univariate polynomial p
into irreducibles where p is known to be square free and primitive with
respect to its main variable.

Isn't it strange that it must be mentioned that p is square-free when it
is supposed to be irreducible? Are there irreducible polynomials that
are not square-free?

Ralf

=================================================================
(1) -> Z==>Integer;Q==>Fraction Z;SUP==>SparseUnivariatePolynomial
Type:
Void
(2) -> f := (((x-1)*x)^2) :: SUP Z

4 3 2
(2) ? - 2 ? + ?
Type:
SparseUnivariatePolynomial(Integer)
(3) -> squareFree f

2 2
(3) (? - ?)
Type:
Factored(SparseUnivariatePolynomial(Integer))
(4) -> squareFreePolynomial f

2 2
(4) (? - ?)
Type:
Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Integer)))
(5) -> squareFreePolynomial(f)$Z

2 2
(5) (? - ?)
Type:
Factored(SparseUnivariatePolynomial(Integer))
(6) -> factor f

2 2
(6) (? - 1) ?
Type:
Factored(SparseUnivariatePolynomial(Integer))
(7) -> factorSquareFreePolynomial f

2 2
(7) (? - 1) ?
Type:
Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Integer)))
(8) -> factorSquareFreePolynomial(f)$Z
2 2 <menu-bar> <signals> <break>
>> System error:
Interactive interrupt at #x5261AA8C.

(8) -> )clear prop f
(8) -> f := (((x-1)*x)^2) :: SUP Q

4 3 2
(8) ? - 2 ? + ?
Type:
SparseUnivariatePolynomial(Fraction(Integer))
(9) -> squareFree f

2 2
(9) (? - ?)
Type:
Factored(SparseUnivariatePolynomial(Fraction(Integer)))
(10) -> squareFreePolynomial f

2 2
(10) (? - ?)
Type:
Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Fraction(Integer))))
(11) -> squareFreePolynomial(f)$Q
Internal Error
The function squareFreePolynomial with signature hashcode is missing
from domain Fraction(Integer)

(11) -> factor f

2 2
(11) (? - 1) ?
Type:
Factored(SparseUnivariatePolynomial(Fraction(Integer)))
(12) -> factorSquareFreePolynomial f

2 2
(12) (? - 1) ?
Type:
Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Fraction(Integer))))
(13) -> factorSquareFreePolynomial(f)$Q
<menu-bar> <signals> <break>
>> System error:
Interactive interrupt at #x52100D0C.

[Waldek Hebisch]

Ralf Hemmecke wrote:
>
> There is seemingly a problem with the implementation of some functions
> from PolynomialFactorizationImplizit in Integer as well as in
> FractionInteger.

Well, some parts of machinery related to PolynomialFactorizationExplicit
are still unimplemented. ATM my tactic is to avoid puting in a lot
of new code, rather I would like to reduce redundancy and eliminate
some bugs:
- several domains have specialized polynomial factorizers. They
are faster than general code so I do not want to simply drop
them. Rather, I would like to transfer speed improvements to
general code and only after that drop old factorizers
- positive characteristic case has large unimplemented part
(here relation between square-free factorization and full one
is more tricky).
- current code assumes that it can factor polynomials in 3 or more
variables over finite fields without using field extentions
for evaluation points. IIUC this is false, but counter examples
are very rare.
- for speed 'solveLinearPolynomialEquation' and 'gcdPolynomial'
should use modular methods if possible. I believe that using
Zippel method would give us large speed improvements (there
are newer things than Zippel, it is not clear fro me if they
give worthwile improvement over Zippel). Current general
case uses fractions which is slower than almost any alternative.
In particular I do not want to use this code when we have
better specialized version.
- in general, PolynomialFactorizationExplicit is build via
recursion. This is nice from point of view of high level
coding, but means that some expensive computations at
lower level need to be done multiple times. IIUC original
authors assumed that caching can eliminate cost of multiple
evaluation, but caching has its own problems and when we
make random choices it simply does not work.

In fact I was thinking about something which could be called
AlgebraicallyExplicit, that is domains for which we can get
information about generators and relations. In other words,
we can build explicit partial izomorphizm with subset of
algebraic extention of field of rational functions.

> In fact, I can only guess what the function "factorSquareFreePolynomial"
> should actually return.
>
> http://fricas.github.io/api/PolynomialFactorizationExplicit.html#l-polynomial-factorization-explicit-factor-square-free-polynomial
>
> factorSquareFreePolynomial: SparseUnivariatePolynomial % -> Factored
> SparseUnivariatePolynomial %
> factorSquareFreePolynomial(p) factors the univariate polynomial p
> into irreducibles where p is known to be square free and primitive with
> respect to its main variable.
>
> Isn't it strange that it must be mentioned that p is square-free when it
> is supposed to be irreducible? Are there irreducible polynomials that
> are not square-free?

Well, each function has preconditions, that is conditions that
is input must satisfy, and postconditions that is what output
satisfies. In this case:

Precondition:

p is square-free and primitive

Postcondition:

output is factorization of p into irreducibles.

If you know a little about factorization precondition is natural:
core factorization methods only work for square-free and primitive
polynomials. Full factorization routine first computes square-free
factorization and make polynomial primitive. Only after that
calls core factorizer. If we know that polynomial is square-free
and primitive, we can skip first part and directly go to core
factorization. In other words, 'factorSquareFreePolynomial' is
core factorizer and before calling it you may need some proprocessing.

You violated precondition, so can not expect sensible result.
You may get one (this case probably went trough full factorizer),
but there is no warranty.

> (8) -> factorSquareFreePolynomial(f)$Z
> 2 2 <menu-bar> <signals> <break>
> >> System error:
> Interactive interrupt at #x5261AA8C.

Again, you violated precondition, so can not expect sensible result.
In fact, if you look at implementation it is clear that infinite
loop is expected here.

> (8) -> )clear prop f
> (8) -> f := (((x-1)*x)^2) :: SUP Q
>
> 4 3 2
> (8) ? - 2 ? + ?
> Type:
> SparseUnivariatePolynomial(Fraction(Integer))
> (9) -> squareFree f
>
> 2 2
> (9) (? - ?)
> Type:
> Factored(SparseUnivariatePolynomial(Fraction(Integer)))
> (10) -> squareFreePolynomial f
>
> 2 2
> (10) (? - ?)
> Type:
> Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Fraction(Integer))))
> (11) -> squareFreePolynomial(f)$Q
> Internal Error
> The function squareFreePolynomial with signature hashcode is missing
> from domain Fraction(Integer)

Well, Hyperdoc clearly shows that it is unimplemented.

> (11) -> factor f
>
> 2 2
> (11) (? - 1) ?
> Type:
> Factored(SparseUnivariatePolynomial(Fraction(Integer)))
> (12) -> factorSquareFreePolynomial f
>
> 2 2
> (12) (? - 1) ?
> Type:
> Factored(SparseUnivariatePolynomial(SparseUnivariatePolynomial(Fraction(Integer))))
> (13) -> factorSquareFreePolynomial(f)$Q
> <menu-bar> <signals> <break>
> >> System error:
> Interactive interrupt at #x52100D0C.
>

> --
> 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/21f9e895-d148-ef07-f8f8-335d95f07304%40hemmecke.org.
>
> --------------65B11D43D9146833BA47C104
> Content-Type: text/plain; charset=UTF-8;
> name="sfp-bug.input"
> Content-Transfer-Encoding: base64
> Content-Disposition: attachment;
> filename="sfp-bug.input"
>
> KWNsZWFyIGNvbXBsZXRlClo9PT5JbnRlZ2VyO1E9PT5GcmFjdGlvbiBaO1NVUD09PlNwYXJz
> ZVVuaXZhcmlhdGVQb2x5bm9taWFsCmYgOj0gKCgoeC0xKSp4KV4yKSA6OiBTVVAgWgpzcXVh
> cmVGcmVlIGYKc3F1YXJlRnJlZVBvbHlub21pYWwgZgpzcXVhcmVGcmVlUG9seW5vbWlhbChm
> KSRaCmZhY3RvciBmCmZhY3RvclNxdWFyZUZyZWVQb2x5bm9taWFsIGYKZmFjdG9yU3F1YXJl
> RnJlZVBvbHlub21pYWwoZikkWgoKKWNsZWFyIHByb3AgZgpmIDo9ICgoKHgtMSkqeCleMikg
> OjogU1VQIFEKc3F1YXJlRnJlZSBmCnNxdWFyZUZyZWVQb2x5bm9taWFsIGYKc3F1YXJlRnJl
> ZVBvbHlub21pYWwoZikkUQpmYWN0b3IgZgpmYWN0b3JTcXVhcmVGcmVlUG9seW5vbWlhbCBm
> CmZhY3RvclNxdWFyZUZyZWVQb2x5bm9taWFsKGYpJFEK
> --------------65B11D43D9146833BA47C104--
>


--
Waldek Hebisch

[Ralf Hemmecke]

>> (8) -> factorSquareFreePolynomial(f)$Z
>> 2 2 <menu-bar> <signals> <break>
>> >> System error:
>> Interactive interrupt at #x5261AA8C.
>
> Again, you violated precondition, so can not expect sensible result.
> In fact, if you look at implementation it is clear that infinite
> loop is expected here.

Thank you, Waldek. Actually, when I prepared my mail I got the
impression that factorSquareFreePolynomial is actually meant in the way
you describe it, i.e., the input must be already square-free. Arrrhhh, I
should have paid more attention to the exact formulation of the
++-docstring. I must have been reading with closed eyes.

I was looking for a square-free factorizer and had chosen the wrong
function. :-(

Anyway, unimplemented code is bad (although I sometimes also appreciate
that I don't have to implement every function before I can play with the
functionality of my new domain that I am actually interested in).

Yes, it looks that factorization is implemented through quite some
domains/packages. I would be happy if the algebra library were clearer.

I also see SparseUnivariatePolynomial and NewSparseUnivariatePolynomial.
The latter is, in fact, an extension of SUP. I have the impression that
it would make sense to merge these domain constructors into SUP.
Unfortunately, I've not yet looked into this too deeply, but I guess,
just moving over the functionality of NewSUP into SUP would be the way
to go. In fact, here it would be nice to have Aldor's "extend" option,
then one could simply "call" NewSUP by the name SUP and still have
separate compilation. Dream, dream. ;-)

Ralf

[Ralf Hemmecke]

> Well, some parts of machinery related to PolynomialFactorizationExplicit
> are still unimplemented.

I guess it is not clear to everyone that PolynomialFactorizationExplicit
is a general idea to create efficient code for certain domains.

That is a domain R implements a function for a domain where R is a kind
of coefficient domain. So Integer implements a function

squareFreePolynomial: SUP(Integer) -> Factored SUP Integer

for SUP = SparseUnivariatePolynomial.

Then the implementation of squareFree in SparseUnivariatePolynomial(R)
where R has PolynomialFactorizationExplicit can simply be

squareFree(x: %): Factored % == squareFreePolynomial(x)$F

and that should be efficient, because specific properties of the
coefficient domain are dealt with in the coefficient domain itself.

Unfortunately, it is not done that way in FriCAS. I see from HyperDoc,
then squareFree of SUP(Integer) is implemented by the default package of
UnivariatePolynomialCategory and there it is:

if R has GcdDomain then
squareFree(p : %) ==
squareFree(p)$UnivariatePolynomialSquareFree(R, %)

Ralf

[Waldek Hebisch]

Ralf Hemmecke wrote:
>
> > Well, some parts of machinery related to PolynomialFactorizationExplicit
> > are still unimplemented.
>
> I guess it is not clear to everyone that PolynomialFactorizationExplicit
> is a general idea to create efficient code for certain domains.

Well, general idea in FriCAS is that domain and conditional sections
can provide more efficient code in special cases. However,
design of PolynomialFactorizationExplicit has serious efficiency
problems. Virtue of PolynomialFactorizationExplicit is
simplicity.



--
Waldek Hebisch

[Ralf Hemmecke]

> Well, general idea in FriCAS is that domain and conditional sections
> can provide more efficient code in special cases.

Yes, I know, but then one has to test for every possible "coefficient"
domain. That doesn't really scale well.

Take for example, SUP. Once it is implemented, one has to add special
code for any new coefficient domain (in case one wants that).

Well we are luckily open source so that's not a big problem, but still,
it's somewhat ugly.

> However, design of PolynomialFactorizationExplicit has serious
> efficiency problems. Virtue of PolynomialFactorizationExplicit is
> simplicity.

Yes simplicity is what I like in this idea. Now you say that there is an
efficience panelty. Can you elaborate on that? I somehow don't quite see
it. Is this because the compiler is not able to simple take the function
from the coefficient domain (for example squareFreePolynomial$R) and
inline it as squareFree$SUP(R)?

Since such a design is also in Aldor's algebra library, maybe Peter
Broadbery can say a few words how this is handled in Aldor.

Ralf

[Waldek Hebisch]

Ralf Hemmecke wrote:
>
> > Well, general idea in FriCAS is that domain and conditional sections
> > can provide more efficient code in special cases.
>
> Yes, I know, but then one has to test for every possible "coefficient"
> domain. That doesn't really scale well.
>
> Take for example, SUP. Once it is implemented, one has to add special
> code for any new coefficient domain (in case one wants that).
>
> Well we are luckily open source so that's not a big problem, but still,
> it's somewhat ugly.

My point is that PolynomialFactorizationExplicit really does not
give you more possibilities than that.

> > However, design of PolynomialFactorizationExplicit has serious
> > efficiency problems. Virtue of PolynomialFactorizationExplicit is
> > simplicity.
>
> Yes simplicity is what I like in this idea. Now you say that there is an
> efficience panelty. Can you elaborate on that? I somehow don't quite see
> it. Is this because the compiler is not able to simple take the function
> from the coefficient domain (for example squareFreePolynomial$R) and
> inline it as squareFree$SUP(R)?

It would take too much time to explain it well. But basicaly you
need global view for efficient factorization. Inability to
inline is just tip of iceberg. One example problem (not the
most important, but easiest to explain) is clearing denominators.
With global view this can be done once and _all_ other steps
work without fractions. With recursive approach there may be
fractions hidden in inner domains and at upper levels there
is no way to avoid them. Another issue is variable ordering:
simple heuristic of using variable of lowest degree as main
variable helps a lot. But you can not do this when variables
are hidden in inner domains.

As an excercise you can think how to organize routines in
src/algebra/amodgcd.spad as nice recursion in style of
PolynomialFactorizationExplicit. You will see that
routines depend crucialy on information that is normaly
hidden by domains. And formulating explicitely properties
of domains leads to so called principal ideal rings (PIR)
which are not integral domains, but have some useful
properties of PID-s.

--
Waldek Hebisch