[fricas-devel] Noncommutative factorization

[Waldek Hebisch]

I looked ate noncommutative factorization code and AFAICS
'xdpolyf1.spad' has serious problem. One example is:

(58) -> factor((x^2 - 2)*(y - 1)*(x - 1))

2 2 3 2
(58) [- 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))

that is '(x^2 - 2)*(y - 1)*(x - 1)' is treated as irreducible.

More generally, factorization via equation solving directly
gives absolute factorization, that is factorization over algebraic
closure of base field. To get factorization over base field
one needs to recombine factors. IIUC 'xdpolyf1.spad' tries
various tricks to avoid algebraic extentions, but this is
very unlikely to work in general.

[Ray]

Post by Waldek Hebisch
I looked ate noncommutative factorization code and AFAICS
(58) -> factor((x^2 - 2)*(y - 1)*(x - 1))
2 2 3 2
(58) [- 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
that is '(x^2 - 2)*(y - 1)*(x - 1)' is treated as irreducible.
More generally, factorization via equation solving directly
gives absolute factorization, that is factorization over algebraic
closure of base field. To get factorization over base field
one needs to recombine factors. IIUC 'xdpolyf1.spad' tries
various tricks to avoid algebraic extentions, but this is
very unlikely to work in general.

I tried this on my saved version (part of a test -harness) and it works
correctly.
Here is the result
(52) -> aa:=(x^2 - 2)*(y - 1)*(x - 1)

2 2 3 2
(52) - 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x
Type:
XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer)))
(53) -> factor(aa)

2
(53) [- 2 + x , 1 - y, 1 - x]
Type:
List(XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer))))
--
Attached is a test file.
Let me know if you are interested in a full test suite and harness?
I also have Konrad Schrempf's next to last entry solving factoring.
I have personal copies (more than I need) of both with a plethora of
testing :)
I defined "aa" to make sure there was no "cheating"; i.e. inadvertent
peeking by any program

RayR

[Bill Page]

Post by Ray

I looked at noncommutative factorization code and AFAICS
(58) -> factor((x^2 - 2)*(y - 1)*(x - 1))
2 2 3 2
(58) [- 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
that is '(x^2 - 2)*(y - 1)*(x - 1)' is treated as irreducible.

Thank you for this example. As noted in the source code the purpose
of the heuristic for finding variables that can be set to zero is to
provide a shortcut that makes the solution more efficient for larger
systems of equations but this heuristic can sometimes fail to produce
a consistent set of equations. This is solved by considering
successively smaller subsets of variables that might be set to zero.
Your example triggers a case where the existing code fails to try hard
enough. The following patch corrects this problem:

~/ncpoly$ git diff
diff --git a/xdpolyf1.spad b/xdpolyf1.spad
index d91cf4a..dc1d002 100644
--- a/xdpolyf1.spad
+++ b/xdpolyf1.spad
@@ -145,16 +145,24 @@ XDistributedPolynomialFunctions1(ALPHABET:List
Symbol, F:Join(IntegralDomain,Gcd

else
s0: List Equation G := []
- while empty? s0 for fp in v repeat
- --output("try: ", fp::OutputForm)
- s := solve(e1,remove(fp,v))$SystemSolvePackage(F)
- while empty? s0 for s1 in s repeat
- m := map(explicit, s1)
- if #s1>0 and reduce(_and$Boolean, m) then s0:=s1
- --output("s0: ",s0::OutputForm)
-
- if empty? s0 then
- return [p]
+ while empty? s0 repeat
+ while empty? s0 for fp in v repeat
+ --output("try: ", fp::OutputForm)
+ s := solve(e1,remove(fp,v))$SystemSolvePackage(F)
+ while empty? s0 for s1 in s repeat
+ m := map(explicit, s1)
+ if #s1>0 and reduce(_and$Boolean, m) then s0:=s1
+ --output("s0: ",s0::OutputForm)
+
+ if empty? s0 then
+ if empty? lz then
+ return [p]
+ else
+ -- try harder!
+ lz := bisect(lz,e)
+ e1 := eval(e,lz)
+ v := members set(concat map(variables, e1))$Set(Symbol)
+
sz := concat(lz,s0)

-- choose a parameter value to make G retractable to F

with the following result:

(45) -> factor((x^2 - 2)*(y - 1)*(x - 1))

2
(45) [2 - x , 1 - y, - 1 + x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 0.09 (EV) + 0.01 (OT) = 0.10 sec

Post by Ray

More generally, factorization via equation solving directly
gives absolute factorization, that is factorization over algebraic
closure of base field. To get factorization over base field
one needs to recombine factors.

I am not convinced that this is the case.

Post by Ray

IIUC 'xdpolyf1.spad' tries
various tricks to avoid algebraic extentions, but this is
very unlikely to work in general.

The tricks in xdpolyf1 having nothing to do with avoiding algebraic
extensions. radicalSolve is only called if the base ring has
RadicalCategory, otherwise SystemSolvePackage only looks for solutions
in fraction field and xdpolyf1 lifts that solution to the base ring.

Are you claiming that this does not provide the most general
factorization in the base ring? Are you suggesting that if one looked
for factorizations over the algebraic closure and then recombined some
factors it might be possible to general polynomials over the base ring
that are not considered when directly solving the factorization
equations over the fraction field?

Post by Ray
I tried this on my saved version (part of a test -harness) and it works
correctly.
Here is the result
(52) -> aa:=(x^2 - 2)*(y - 1)*(x - 1)
2 2 3 2
(52) - 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x
XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer)))
(53) -> factor(aa)
2
(53) [- 2 + x , 1 - y, 1 - x]
List(XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer))))
--

Note that you are not using the same polynomial base ring as Waldek,
i.e. Fraction(MultivariatePolynomial([a,b,c,d,e],Integer))) versus
Integer.

[Ray]

Post by Bill Page

Post by Ray
I tried this on my saved version (part of a test -harness) and it works
correctly.
Here is the result
(52) -> aa:=(x^2 - 2)*(y - 1)*(x - 1)
2 2 3 2
(52) - 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x
XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer)))
(53) -> factor(aa)
2
(53) [- 2 + x , 1 - y, 1 - x]
List(XDistributedPolynomial(OrderedVariableList([x,y,z]),Fraction(MultivariatePolynomial([a,b,c,d,e],Integer))))
--

Note that you are not using the same polynomial base ring as Waldek,
i.e. Fraction(MultivariatePolynomial([a,b,c,d,e],Integer))) versus
Integer.

Yes I did.

RayR

[Waldek Hebisch]

Post by Bill Page

Post by Waldek Hebisch
More generally, factorization via equation solving directly
gives absolute factorization, that is factorization over algebraic
closure of base field. To get factorization over base field
one needs to recombine factors.

I am not convinced that this is the case.

Post by Waldek Hebisch
IIUC 'xdpolyf1.spad' tries
various tricks to avoid algebraic extentions, but this is
very unlikely to work in general.

The tricks in xdpolyf1 having nothing to do with avoiding algebraic
extensions. radicalSolve is only called if the base ring has
RadicalCategory, otherwise SystemSolvePackage only looks for solutions
in fraction field and xdpolyf1 lifts that solution to the base ring.
Are you claiming that this does not provide the most general
factorization in the base ring? Are you suggesting that if one looked
for factorizations over the algebraic closure and then recombined some
factors it might be possible to general polynomials over the base ring
that are not considered when directly solving the factorization
equations over the fraction field?

If you could find solution _in the fraction field_ then the
method would be fine. However, in general finding rational
solutions to polynomial system of equations is uncomputable.
Groebner bases decide if there are solutions in algebraic
closure, but you may have algebraic solutions without
rational solutions. If you say that you can find out
if there is rational solution (= factorization) you should
better justify this and explain what special properties
of system you use.

Classical factorization methods (in commutative case)
avoid uncomputability by recombining algebraic factors.
If you do not want recombination you should say how
you avoid uncomputability.

To be clear: your code apparently makes some assumptions.
Both theory (uncomputablity in general case) an practice
suggest that those assumptions may fail unless we can
prove them.

[Bill Page]

Post by Waldek Hebisch
If you could find solution _in the fraction field_ then
the method would be fine. However, in general finding
rational solutions to polynomial system of equations is
uncomputable.

Can you suggest a reference? I could not find this result (well
known?) after a quick search.

Post by Waldek Hebisch
Groebner bases decide if there are solutions in
algebraic closure, but you may have algebraic
solutions without rational solutions. If you say that
you can find out if there is rational solution
(= factorization) you should better justify this and
explain what special properties of system you use.

Yes, that would be nice. Unfortunately SystemSolvePackage in FriCAS
does not make any explicit claim about completeness. But it does refer
to the method of triangular systems.

The only special property that I can think of is that the equations in
the system are at most quadratic. I do not know if that is sufficient.
Another thing is that we are only interested in finding at least one
explicit solution (if it exists). We do not need to know all
solutions.

Post by Waldek Hebisch
Classical factorization methods (in commutative case)
avoid uncomputability by recombining algebraic factors.
If you do not want recombination you should say how
you avoid uncomputability.

I agree.

Post by Waldek Hebisch
To be clear: your code apparently makes some
assumptions. Both theory (uncomputablity in general
case) an practice suggest that those assumptions may
fail unless we can prove them.

OK.

[Ray]

Post by Bill Page

Post by Waldek Hebisch
If you could find solution _in the fraction field_ then
the method would be fine. However, in general finding
rational solutions to polynomial system of equations is
uncomputable.

Can you suggest a reference? I could not find this result (well
known?) after a quick search.

I believe that the Godel problem as implemented by Chaitine actually
constructed a Diophantine equation (and program) that provably couldn't
be proven to have a finite or infinite number of integer solutions.
Since any "solution" could be considered a "factor" then this would
present a problem to a factoring program.
Along the same lines:
From:
https://en.wikipedia.org/wiki/Undecidable_problem#Examples_of_undecidable_problems
"
A Diophantine equation is a more general case of Fermat's Last Theorem;
we seek the integer roots of a polynomial in any number of variables
with integer coefficients. Since we have only one equation but n
variables, infinitely many solutions exist (and are easy to find) in the
complex plane; however, the problem becomes impossible if solutions are
constrained to integer values only. Matiyasevich showed this problem to
be unsolvable by mapping a Diophantine equation to a recursively
enumerable set and invoking Gödel's Incompleteness Theorem.
"
And this is over a fully commutative ring; integers!
RayR

[Waldek Hebisch]

Post by Ray
Attached is a test file.
Let me know if you are interested in a full test suite and harness?
I also have Konrad Schrempf's next to last entry solving factoring.
I have personal copies (more than I need) of both with a plethora of
testing :)
I defined "aa" to make sure there was no "cheating"; i.e. inadvertent
peeking by any program

Well, I would like to include noncommutative code in FriCAS.
There were several versions showing up on the list. I need
to find the "correct version" and collect enough tests.
I have files that you posted in the past. ATM the problem
is that I have too many files with several near duplicates
(small differences between files) -- for inclusion I need
a single version with hopefully all fixes.

[Waldek Hebisch]

Post by Bill Page

Post by Waldek Hebisch
If you could find solution _in the fraction field_ then
the method would be fine. However, in general finding
rational solutions to polynomial system of equations is
uncomputable.

Can you suggest a reference? I could not find this result (well
known?) after a quick search.

Matiasevicz theorem (solution to Hibert 10-th problem) says
that existence of integer solutions is uncomputable.
Concerning rationals I probably misrememberd things:
Wikipedia says that it is unknown if existence of
rational solutions is computable or not.

Post by Bill Page

Post by Waldek Hebisch
Groebner bases decide if there are solutions in
algebraic closure, but you may have algebraic
solutions without rational solutions. If you say that
you can find out if there is rational solution
(= factorization) you should better justify this and
explain what special properties of system you use.

Yes, that would be nice. Unfortunately SystemSolvePackage in FriCAS
does not make any explicit claim about completeness. But it does refer
to the method of triangular systems.

IIRC it uses either Groebner bases or triangular systems. In both
cases you can get algebraic solutions (irrational ones).

Post by Bill Page
The only special property that I can think of is that the equations in
the system are at most quadratic. I do not know if that is sufficient.

No. There is old trick (Veronese embedding) which reduces general
systems to systems of degree 2.

Post by Bill Page
Another thing is that we are only interested in finding at least one
explicit solution (if it exists). We do not need to know all
solutions.

One useful thing is Davenport observation: one can get solutions
for highest order terms in combinatorial way. Given high order
terms the rest seem to reduce to sequence of linear systems.

Actually, I would like to know some hard examples: all that we
have now seem to be quite easy by ad-hoc methods.

[Waldek Hebisch]

I looked more at the problem. It seems that Cohn claims
that there are finitely many factorizations, but he jumps
over few subtle points so I need to check his arguments
more carefuly.

AFAICS Caruso (after Davenport) gives correct proof that
factorization of homogeneous polynomials is unique and gives
correct algorithm for factoring homogeneous polynomials.
OTOH what he writes about noncommutative polynomials
is problematic. Namely, he tries to choose monomials
and that looks wrong. AFAICS one gets correct version
of algorithm looking at factorization of leading term,
treating it is sequence of primes. So left factor
coresponds to initial part of sequence, call it h, right
factor to the rest, call it g. If there is no overlap
between h and g, then one get complete factorizations
solving linear equations. Each overlap potentially leads
to non-uniqueness, introducing one parameter family
of solutions. So in general we may get solution depending
polynomially on k parameters when k is number of overlaps.
But linear system obtained are overdetermined and there
is good chance that we can quickly determine some parameters
from other equations. So we really get k parameters only
when there is very special structure (probably quite close
to commutative). Once we obtained candidate solution
from linear equations we can multiply candidate factors
and get system of nonlinear equations in at most k variables.
ATM in final step I do not know better method than
convertion to triangular system.

In a sense what I write above is equivalent to what xdpolyf1
is doing. However, xdpolyf1 will use _much_ larger number
of variables and hopes that code computing Groebner bases
can cope with them. Above we can use special structure of
systems and have essentially linear time solver. When we
have several parameters and high degree we may get many
linear systems so method above may get expensive. However
using Groebner bases for this looks much more expensive.

Let me say that when highest order term has nontrivial
commutative image, then we can use commutative factorization
and get combinatorial problem of matching commutative
factors to noncommutative ones. AFAICS any matching
uniquely determines parameters of linear systems, so
to find all factorization we need to check all matchings.
In principle there may be exponentially many matchings,
so this may be expensive. But we should be able to
find quickly low degree factors.

From theoretical point of view in this case we
get direct proof that there are finitely many factorizations.
I hope that one can say more, but enough for today.

[Waldek Hebisch]

I looked at noncommutative factorization code and AFAICS
(58) -> factor((x^2 - 2)*(y - 1)*(x - 1))
2 2 3 2
(58) [- 2 + 2 y + 2 x - 2 y x + x - x y - x + x y x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
that is '(x^2 - 2)*(y - 1)*(x - 1)' is treated as irreducible.

<snip>
Before the patch

f101 := (x*z - z*x)^2 - 2

was immediately recognized as irreducible. With the patch I did not
get answer for several minutes (may be I am not patient enough?).
Similarly

f102 := (x*z + z*x)^2 - 2

and

f103 := (x^2 + x + z)^2 - 2

[Bill Page]

On Sun, Nov 4, 2018 at 8:03 AM Waldek Hebisch <***@math.uni.wroc.pl> wrote:
...

Post by Waldek Hebisch
<snip>
Before the patch
f101 := (x*z - z*x)^2 - 2
was immediately recognized as irreducible. With the patch I did not
get answer for several minutes (may be I am not patient enough?).
Similarly
f102 := (x*z + z*x)^2 - 2
and
f103 := (x^2 + x + z)^2 - 2

Thanks again for these great test cases. There was an error in the
'bisect' function that caused an infinite loop triggered by the
earlier patch. Here is a revised patch that corrects this problem.
(Only one additional change at the beginning.)

---

diff --git a/xdpolyf1.spad b/xdpolyf1.spad
index d91cf4a..cc88597 100644
--- a/xdpolyf1.spad
+++ b/xdpolyf1.spad
@@ -75,7 +75,7 @@ XDistributedPolynomialFunctions1(ALPHABET:List
Symbol, F:Join(IntegralDomain,Gcd
remove(0,map((x:G):G+->eval(x,z)$RF,e))

bisect(z:List Equation G, e:List G):List Equation G ==
- if #z<2 then return z
+ if #z<2 then return []
z1 := first(z,#z quo 2)
if groebner(map(numer,eval(e,z)))~=[1] then return z1
z1 := last(z,#z quo 2)
@@ -145,16 +145,25 @@ XDistributedPolynomialFunctions1(ALPHABET:List
Symbol, F:Join(IntegralDomain,Gcd

else
s0: List Equation G := []
- while empty? s0 for fp in v repeat
- --output("try: ", fp::OutputForm)
- s := solve(e1,remove(fp,v))$SystemSolvePackage(F)
- while empty? s0 for s1 in s repeat
- m := map(explicit, s1)
- if #s1>0 and reduce(_and$Boolean, m) then s0:=s1
- --output("s0: ",s0::OutputForm)
-
- if empty? s0 then
- return [p]
+ while empty? s0 repeat
+ while empty? s0 for fp in v repeat
+ --output("try: ", fp::OutputForm)
+ if not groebner(map(numer,e1))=[1] then
+ s := solve(e1,remove(fp,v))$SystemSolvePackage(F)
+ while empty? s0 for s1 in s repeat
+ m := map(explicit, s1)
+ if #s1>0 and reduce(_and$Boolean, m) then s0:=s1
+ --output("s0: ",s0::OutputForm)
+
+ if empty? s0 then
+ if empty? lz then
+ return [p]
+ else
+ -- try harder!
+ lz := bisect(lz,e)
+ e1 := eval(e,lz)
+ v := members set(concat map(variables, e1))$Set(Symbol)
+
sz := concat(lz,s0)

-- choose a parameter value to make G retractable to F

---

[Waldek Hebisch]

Post by Bill Page
earlier patch. Here is a revised patch that corrects this problem.
(Only one additional change at the beginning.)

I have tried:

h3 := x*y*z - x*z*y + z*x*y - z*y*x + y*z*x - y*x*z
factor(h3*h3)

and after about hour I did not get answer.

Since nobody seems to be interested in coding Davenport method
I did that. At

http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04b.input

dcfact.input is actual routine. nc_ini04b.input defines needed
types. With that things like

h2 := x*y - y*x
homo_fact(h2*h3*h2*h3*h2)

work.

[Bill Page]

Post by Waldek Hebisch

Post by Bill Page
earlier patch. Here is a revised patch that corrects this problem.
(Only one additional change at the beginning.)

h3 := x*y*z - x*z*y + z*x*y - z*y*x + y*z*x - y*x*z
factor(h3*h3)
and after about hour I did not get answer.

Yes. This is a good example where the heuristic fails.

Post by Waldek Hebisch
Since nobody seems to be interested in coding Davenport method
I did that. At
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04b.input
dcfact.input is actual routine. nc_ini04b.input defines needed
types. With that things like
h2 := x*y - y*x
homo_fact(h2*h3*h2*h3*h2)
work.

Thanks, that is very nice. And it is impressively fast on all the
examples I tried. However it did find this example where homo_fact
fails to find a factorization:

(110) -> factor((x^2-1)^2)

(110) [1 - x, 1 - x, 1 + x, 1 + x]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.44 (EV) + 0.01 (OT) = 0.45 sec
(111) -> homo_fact((x^2-1)^2)

2 4
(111) [1 - 2 x + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 0.00 (OT) = 0.00 sec

--

It is also interesting to compare the results obtained by Konrad
Schrempf's algorithm. (Start with 'nc_ini14.input' from the ncpoly
repository.)

[Bill Page]

Post by Bill Page

Post by Waldek Hebisch
Since nobody seems to be interested in coding Davenport method
I did that.

...
(111) -> homo_fact((x^2-1)^2)
2 4
(111) [1 - 2 x + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 0.00 (OT) = 0.00 sec

Yes I see that this is not a homogeneous polynomial so Davenport's
algorithm does not handle it correctly.

[Ray]

Post by Bill Page

Post by Bill Page

Post by Waldek Hebisch
Since nobody seems to be interested in coding Davenport method
I did that.

...
(111) -> homo_fact((x^2-1)^2)
2 4
(111) [1 - 2 x + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 0.00 (OT) = 0.00 sec

Yes I see that this is not a homogeneous polynomial so Davenport's
algorithm does not handle it correctly.

So homo_fact((x^2-y^2)^2)
would succeed?

[Bill Page]

Post by Ray
...
So homo_fact((x^2-y^2)^2)
would succeed?

Yes.

(66) -> homo_fact((x^2-y^2)^2)

2 2 2 2
(66) [- y + x , - y + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0 sec

[Bill Page]

Surprisingly a non-homogeneous polynomial of the same degree works OK

(69) -> factor((h3+1)*(h3+1))

(69)
[1 - z y x + z x y + y z x - y x z - x z y + x y z,
1 - z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 15.93 (EV) + 0.03 (OT) = 15.96 sec

but as Waldek noted 'factor(h3*h3)' does not complete in reasonable time.

Post by Bill Page

Post by Ray
...
So homo_fact((x^2-y^2)^2)
would succeed?

Yes.
(66) -> homo_fact((x^2-y^2)^2)
2 2 2 2
(66) [- y + x , - y + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0 sec

[Waldek Hebisch]

Post by Bill Page

Post by Bill Page

Post by Waldek Hebisch
Since nobody seems to be interested in coding Davenport method
I did that.

...
(111) -> homo_fact((x^2-1)^2)
2 4
(111) [1 - 2 x + x ]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.00 (IN) + 0.00 (OT) = 0.00 sec

Yes I see that this is not a homogeneous polynomial so Davenport's
algorithm does not handle it correctly.

Yes, that is only for homogeneous polynomials. In genral case
one needs to implement lift, somewhat like in Caruso preprint
(but Algorithm 2 has problems, one needs to correct it).

[Waldek Hebisch]

One update to what I wrote before. In

J. P. Bell, A. Heinle, and V. Levandovskyy,
On Noncommutative Finite Factorization Domains,
Trans. Amer. Math. Soc. 369 (2017), 2675-2695

there is proof of finite number of factorizations.

I have now implemented the lift part of Davenport-Caruso method.
You fetch code at:

http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact2.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04c.input

As before, dcfact2.input is an actual routine, nc_ini04c.input
set up needed types.

You can try things like:

dc_fact(((x*z - z*x)^2 - 2)*((x*z - z*x)^2 - 3)*((x*z - z*x)^2 - 5))

The idea is like in Caruso preprint, but code differs -- I coded
what follows by working out solutions to lift equations.

The code is unoptimized, for example for homogeneous polynomials
it runs the whole lift, while we know that the factorization must
be homogeneous... Also, code introduces extra parameters
when there is overlap between leading monomials of top factors,
while in many cases this parameter could be immediately eliminated.
Still, it seem to work quite a bit faster than xdpolyf1.spad.

[Bill Page]

On Sat, Nov 10, 2018 at 12:02 PM Waldek Hebisch
<***@math.uni.wroc.pl> wrote:
...

Post by Waldek Hebisch
I have now implemented the lift part of Davenport-Caruso method.
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact2.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04c.input
As before, dcfact2.input is an actual routine, nc_ini04c.input
set up needed types.
dc_fact(((x*z - z*x)^2 - 2)*((x*z - z*x)^2 - 3)*((x*z - z*x)^2 - 5))
The idea is like in Caruso preprint, but code differs -- I coded
what follows by working out solutions to lift equations.
The code is unoptimized, for example for homogeneous polynomials
it runs the whole lift, while we know that the factorization must
be homogeneous... Also, code introduces extra parameters
when there is overlap between leading monomials of top factors,
while in many cases this parameter could be immediately eliminated.
Still, it seem to work quite a bit faster than xdpolyf1.spad.
--

This looks very promising however I did find this apparent failure:

(82) -> h3

(82) - z y x + z x y + y z x - y x z - x z y + x y z
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (OT) = 0.00 sec
(83) -> dc_fact((3+x*z*y)*h3)

(83)
[
- 3 z y x + 3 z x y + 3 y z x - 3 y x z - 3 x z y + 3 x y z - x z y z y x
+
2 2
x z y z x y + x z y z x - x z y x z - x z y x z y + x z y x y z
]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (OT) = 0.00 sec

[Bill Page]

Post by Bill Page
On Sat, Nov 10, 2018 at 12:02 PM Waldek Hebisch
...

Post by Waldek Hebisch
I have now implemented the lift part of Davenport-Caruso method.
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact2.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04c.input
As before, dcfact2.input is an actual routine, nc_ini04c.input
set up needed types.

...

Post by Bill Page
(82) -> h3
(82) - z y x + z x y + y z x - y x z - x z y + x y z
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (OT) = 0.00 sec
(83) -> dc_fact((3+x*z*y)*h3)
(83)
[
- 3 z y x + 3 z x y + 3 y z x - 3 y x z - 3 x z y + 3 x y z - x z y z y x
+
2 2
x z y z x y + x z y z x - x z y x z - x z y x z y + x z y x y z
]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (OT) = 0.00 sec

For comparison here is the output of Konrad Schrempf's program:

(1) -> )r nc_ini14.input
...
(52) -> h3 := x*y*z - x*z*y + z*x*y - z*y*x + y*z*x - y*x*z

(52)
AGGSET
+1 - x - z 0 - y 0 0 0 + + 0 +
| | | |
| 1 0 z 0 0 y 0 | | 0 |
| | | |
| 1 - x 0 y 0 0 | | 0 |
| | | |
| 1 0 0 0 y | | 0 |
| |s = | |
| 1 - z - x 0 | | 0 |
| | | |
| 1 0 x | | 0 |
| | | |
| 1 - z| | 0 |
| | | |
+ 1 + +- 1+
,
MR
,
- z y x + z x y + y z x - y x z - x z y + x y z
Type: FreeDivisionAlgebra(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (IN) + 0.03 (EV) + 0.01 (OT) = 0.04 sec
(53) -> map(x+->x::XDP,factor((3+x*z*y)*h3))

(53) [3 + x z y, - z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (IN) + 0.16 (EV) + 0.01 (OT) = 0.17 sec

[Bill Page]

Post by Bill Page
On Sat, Nov 10, 2018 at 12:02 PM Waldek Hebisch
...

Post by Waldek Hebisch
I have now implemented the lift part of Davenport-Caruso method.
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact2.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04c.input
As before, dcfact2.input is an actual routine, nc_ini04c.input
set up needed types.

...

Here is a more minimal example of the problem:

(86) -> h3 := x*y*z - x*z*y + z*x*y

(86) z x y - x z y + x y z
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (EV) + 0.00 (OT) = 0.00 sec
(87) -> dc_fact(h3*(1+y))

2 2
(87) [z x y - x z y + x y z + z x y - x z y + x y z y]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (EV) + 0.00 (OT) = 0.00 sec

The xdpolyf1 code gives the following result:

(96) -> factor(h3*(1+y))

(96) [z x y - x z y + x y z, 1 + y]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Integer))
Time: 0.14 (EV) + 0.00 (OT) = 0.14 sec

[Bill Page]

Waldek,

Post by Waldek Hebisch
One update to what I wrote before. In
J. P. Bell, A. Heinle, and V. Levandovskyy,
On Noncommutative Finite Factorization Domains,
Trans. Amer. Math. Soc. 369 (2017), 2675-2695
there is proof of finite number of factorizations.
I have now implemented the lift part of Davenport-Caruso method.
...

Since the number of factorizations of a non-commutative polynomial
over a unique factorization domain is finite but not unique there may
be some applications where it maybe interesting to know more than one
or even all possible factorizations. Your current implementation
produces just one factorization. Do you see any opportunity to extend
the Davenport-Caruso method to produce multiple factorizations or a
complete enumeration of factorizations?

I did some experiments with the xdpolyf1 factorizer to produce such
multiple factorizations. This was relatively easy since the solution
algorithm (with the "pruning" heuristic) naturally produces
factorizations in which either the left factor or right factor at each
step has a minimum number of terms. By alternately choosing to
minimize first the right factor and then the left factor it is
possible to explore alternate factorizations. I did not get so far as
to attempt to prove completeness.

[Waldek Hebisch]

Post by Bill Page
On Sat, Nov 10, 2018 at 12:02 PM Waldek Hebisch
...

Post by Waldek Hebisch
I have now implemented the lift part of Davenport-Caruso method.
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact2.input
http://www.math.uni.wroc.pl/~hebisch/fricas/nc_ini04c.input
As before, dcfact2.input is an actual routine, nc_ini04c.input
set up needed types.

(82) -> h3
(82) - z y x + z x y + y z x - y x z - x z y + x y z
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (OT) = 0.00 sec
(83) -> dc_fact((3+x*z*y)*h3)
(83)
[
- 3 z y x + 3 z x y + 3 y z x - 3 y x z - 3 x z y + 3 x y z - x z y z y x
+
2 2
x z y z x y + x z y z x - x z y x z - x z y x z y + x z y x y z
]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (OT) = 0.00 sec

Thanks for testcases (this one and from other post). I overlooked
a case when lift equation has unique solution, but the simple
method used by 'dc_fact1' found wrong one. Now those cases
are hanlded like cases with non-unique solution and passed to
'solve' to sort out. New version at:

http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact3.input

[Bill Page]

That looks great!

As a performance test I tried this:

(79) -> h2323:=((h2*h3*h2*h3)^2);

Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (IN) + 2.71 (EV) + 0.00 (OT) = 2.71 sec
(80) -> dc_fact h2323

(80)
[- y x + x y, - z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (IN) + 315.56 (EV) + 0.13 (OT) = 315.69 sec

Although you have already pointed out that this code is not yet
optimized I think it is interesting to compare this with Konrad's
program that produces the following result on the same problem:

(54) -> h2323:=((h2*h3*h2*h3)^2);

Type: FreeDivisionAlgebra(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (IN) + 0.37 (EV) + 0.00 (OT) = 0.37 sec
(55) -> map(x+->x::XDP,factor h2323)

(55)
[y x - x y, z y x - z x y - y z x + y x z + x z y - x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, y x - x y,
z y x - z x y - y z x + y x z + x z y - x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (IN) + 25.45 (EV) + 0.01 (OT) = 25.46 sec

---

Post by Waldek Hebisch
...
Thanks for testcases (this one and from other post). I overlooked
a case when lift equation has unique solution, but the simple
method used by 'dc_fact1' found wrong one. Now those cases
are hanlded like cases with non-unique solution and passed to
http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact3.input
--

[Waldek Hebisch]

Post by Bill Page
(79) -> h2323:=((h2*h3*h2*h3)^2);
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (IN) + 2.71 (EV) + 0.00 (OT) = 2.71 sec
(80) -> dc_fact h2323
(80)
[- y x + x y, - z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (IN) + 315.56 (EV) + 0.13 (OT) = 315.69 sec
Although you have already pointed out that this code is not yet
optimized I think it is interesting to compare this with Konrad's

Thanks for interestiong observation.

Actually Konrad's program performs quite a different computation.
In particular IIUC Konrad tries to preserve structure of the input.
That is good, but means that routines are hard to compare.

Probably more relevant is comparison to 'homo_fact' which in this
case is doint real work -- for homogeneous polynomial 'dc_fact'
just adds overhead. On my machine

(67) -> homo_fact(h2323)

(67)
[- y x + x y, - z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,^
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.32 (EV) = 0.32 sec

As you can see 'homo_fact' is about 1000 times faster than 'dc_fact'.
You may notice that 'homo_fact' is faster than multiplication that
produced 'h2323'. At first glance it may look reasonable, but
'homo_fact' checks result by multiplication. It turned out
that speed of multiplication depends strongly on order of operations.
Compare:

(70) -> h2323 := a2*(a3*(a2*(a3*(a2*a3*a2*a3))));

Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.12 (EV) = 0.12 sec

(71) -> h2323 := (a2*(a3*(a2*(a3*(a2*a3*a2)))))*a3;

Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 80.09 (EV) = 80.09 sec

'dc_fact' uses about 15-20 sec to multiply factors in unfavourable
order. As a quick hack I added routine to muliplay factors starting
from right, it gave expected speedup. But we probably should work
on noncommutative polynomial multiplication to make it faster.
Current code is in 'xpoly.spad' and main step is:

+/[[[t1.k*t2.k, t1.c*t2.c]$TERM for t2 in p2]
for t1 in p1]

that is we add list of multiples of p2 by single term of p1.
We could probably speed it up by having faster routine to add
list of elements of free abelian monoid. More precisely,
elements of free monoid are represented by ordered list.
For two lists of similar length addition is linear. But when
we add several such lists the "acumulator" gets quite long
and cost become quadratic. We probably could speed it up
by trying to balance length of lists. One idea would be
to add two shortest list as a basic step.

However, the main reason of slowness of this example is
slowness of Groebner bases (used via 'solve'), this takes
about 95% of time. We have system of 34560 equations
in 4 unknowns. This system is highly redundant, there
are only 18 distinct equations. After adding
'removeDuplicates' on list of equations (and workaround
for order of multiplication) on my machine time is down
to 16.25 sec. I did not measure how it splits between
routines, but printouts on the screen suggested that
most time went into 'removeDuplicates'.

Now concerning possible fixes: for users it is rather natural
to create large sets of equations with high redundancy.
So it would be good if 'solve' (or Groebner code) could
handle it. 'removeDuplicates' is not a good solution.
First, it is slow for long lists. And redundancy may
be more subtle than having duplicates: for example we
may have scalar multiples of the same equation. In our
case the system is kind of trivial: two variables can be
determined from linear equations and when values ot those
variables is plugged into other equations one gets enough
linear equations to determine values of two other variables.
Currently 'solve' checks if system is linear, clearly
there is place for smarter handling of linear equations.

Another thing is that 'dc_fact' could spend more effort and
solve linear equations intead of passing them to 'solve'.
In this way we would avoid generationg large number of
useless equations -- the example shows that 'dc_fact'
should do this, but improving 'solve' also would have
benefits.

Of course, for homogeneous polynomials we could directly
use result of 'homo_fact', but AFAICS the same issues
will appear for non-homogeneous polynomials. In fact,
currently 'dc_fact' treats homogeneous polynomials like
all others exactly because I wanted to get more testing
for main part of 'dc_fact' (in particular to get various
edge cases).

--
Waldek Hebisch

[Waldek Hebisch]

Post by Bill Page
That looks great!
(79) -> h2323:=((h2*h3*h2*h3)^2);
Type: XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer))
Time: 0.00 (IN) + 2.71 (EV) + 0.00 (OT) = 2.71 sec
(80) -> dc_fact h2323
(80)
[- y x + x y, - z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 0.00 (IN) + 315.56 (EV) + 0.13 (OT) = 315.69 sec

I have put a new version at:

http://www.math.uni.wroc.pl/~hebisch/fricas/dcfact4.input

Now, on my machine:

(357) -> dc_fact(h2323)

(357)
[- y x + x y, - z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z, - y x + x y,
- z y x + z x y + y z x - y x z - x z y + x y z]
Type: List(XDistributedPolynomial(OrderedVariableList([x,y,z,w,x1,x2,x3,x4,x5]),Fraction(Integer)))
Time: 1.12 (EV) = 1.12 sec

Code tries to solve linear equations if possible, so in this case
does not need to call 'solve'. But there are many special cases
so more space where bugs can hide...

[Waldek Hebisch]

Post by Bill Page
Since the number of factorizations of a non-commutative polynomial
over a unique factorization domain is finite but not unique there may
be some applications where it maybe interesting to know more than one
or even all possible factorizations. Your current implementation
produces just one factorization. Do you see any opportunity to extend
the Davenport-Caruso method to produce multiple factorizations or a
complete enumeration of factorizations?

1) My feeling is that given one factorization it should be easier
to produce other. But that requires some theoretical work.
In particular all examples of non-uniqueness I saw are very
special -- I do not know if there is a general principle or
I just did not saw more tricky examples.

2) ATM 'dc_fact1' produces just one factorization of given degree,
but by using all results from solve that are in base field
we would get all factorizations with given highest degree part.
Currently 'dc_fact' tries low degree left factor first and
do not try higher degree left factors. By calling 'dc_fact1'
with all possible splittings of highest degree part one would
get all possible factorizations into two factors. Recursively
one could get from this all possible factorizations.

[Waldek Hebisch]

I have now put Spad version of factorization code at:

http://www.math.uni.wroc.pl/~hebisch/fricas/xpfact.spad

[Waldek Hebisch]

A new version with some improvements:

http://www.math.uni.wroc.pl/~hebisch/fricas/xpfact2.spad

In particular

factor((x^4 + 5)*(x^4 + x + 7))

is now much faster (previously needed 2397.40 sec on my machine).