Opened 4 years ago

Closed 3 years ago

# multivariate factorization over QQbar

Reported by: Owned by: gh-BrentBaccala major sage-9.0 algebra gh-tom111 Brent Baccala Travis Scrimshaw N/A 91d467e 91d467ee9ee4008074de452de9ac8a6b0b90359a #25351

We can currently (#8544) factor univariate polynomials over QQbar:

```sage: R.<x>=QQbar[]
sage: (x^2+1).factor()
(x - I) * (x + I)
```

We can factor over multivariate rings, if the polynomial is actually univariate:

```sage: R.<x,y>=QQbar[]
sage: (x^2+1).factor()
(x - I) * (x + I)
```

But we can't do full "absolute factorization" (that's what it's called in the literature):

```sage: R.<x,y>=QQbar[]
sage: (x^2+y^2).factor()
...
NotImplementedError: proof = True factorization not implemented.  Call factor with proof=False.
sage: (x^2+y^2).factor(proof=False)
...
TypeError: no conversion of this ring to a Singular ring defined
```

This ticket implements absolute factorization by calling Singular's `absfact.lib`, which works for rings over Q (or transcendental extensions thereof), and uses a technique described in the attached document to handle number fields.

I flagged it `major` because polynomial factorization is a fundamental feature that blocks the implementation of other things, like primary decomposition of ideals.

### comment:1 Changed 4 years ago by gh-BrentBaccala

• Authors set to Brent Baccala
• Branch set to public/25390
• Commit set to 27f109d5b6b0f9cf699218a5af9ea6ccea22a865
• Status changed from new to needs_review

I added code to call Singular's `absFactorize` if the polynomial's coefficients are in `QQ`.

Here's what's now possible:

```sage: R.<x,y> = QQbar[]

sage: L = (x^2+y^2).factor()
sage: L
((-1*I)*x - y) * (1*I*x - y)
sage: L.value()
x^2 + y^2

sage: p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5*x^2+y^2)^3*(x-y)^4
sage: F = p.factor()
sage: F
...
sage: F.value() == p
True
```

But something like this still doesn't work (because the coefficients aren't in `QQ`):

```sage: (x+QQbar(sqrt(2))*y).factor()
```

It's something, which is better than nothing, so I think we should move it towards `master`, but there's still a lot of work to be done before this ticket is closed.

New commits:

 ​27f109d `Trac #25390: use Singular 'absfact' to factor multivariate polynomials over QQ`

### comment:2 Changed 4 years ago by tscrim

I agree that small improvements are good. We can always do followup tickets to add more functionality. I do have some comments:

Can we use `libsingular` instead of calls to `singular` (which uses the pexpect interface and is much slower)?

Doc tweaks:

```         ALGORITHM:

-        Uses Singular's `absfact` library.
+        Uses Singular's ``absfact`` library.

-        TODO:
+        .. TODO::

-        Implement absolute factorization over number fields
+            Implement absolute factorization over number fields.
```

Also if you could split the long output line so that it is at most ~80 chars per line.

Instead of `elem_map`, you could just use `elem_dict.__getitem__`.

I don't see the need for the internal functions `polynomial_map` and `reverse_polynomial_map`. These little trivial functions make the code harder to follow.

In Python, usually `assert`s are with spaces, e.g., `assert minpoly.degree() == 0`.

### comment:3 Changed 4 years ago by gh-BrentBaccala

• Status changed from needs_review to needs_work

### comment:4 Changed 4 years ago by git

• Commit changed from 27f109d5b6b0f9cf699218a5af9ea6ccea22a865 to 747861ef4df567b68e077ae8a47ccc3a44313eea

Branch pushed to git repo; I updated commit sha1. New commits:

 ​e6124e2 `Trac #25390: convoluted code that I want to save before cleaning` ​b34f2b3 `Trac #25390: correctly handle factorization over number fields` ​2879cc0 `Trac #25390: bug fix (and add a test case)` ​747861e `Trac #25390: code cleanup suggested by tscrim's review`

### comment:5 Changed 4 years ago by gh-BrentBaccala

Over the weekend I remembered a technique that I read years ago for factoring over a number field - factor the norm of the polynomial, which has rational coefficients and is a multiple of the original polynomial.

So, this ticket is moving towards `closed` more quickly that I had thought possible. I think it can now factor any polynomial over `QQbar`.

The algorithm is tricky enough that I'm planning to write up a short paper describing it, post it on arxiv, and link to it from the documentation. If anybody has a better suggestion for how to incorporate a three page paper explaining the algorithm of a Sage method, please let me know.

I took most of `tscrim`'s suggestions, but I can't figure out how to access a Singular variable from `libsingular`. I posted a question on asksage, but so far no answer.

I also intend to implement factorization over `AA` before flipping this ticket to `needs_review`.

### comment:6 Changed 4 years ago by git

• Commit changed from 747861ef4df567b68e077ae8a47ccc3a44313eea to 247f94497f5e9bb646d10cfbf439c1154e983c43

Branch pushed to git repo; I updated commit sha1. New commits:

 ​b2657e2 `Trac #25390: simplify factorization logic` ​580c11e `Trac #25390: enhance multivariate factorization to handle AA as well as QQbar` ​be078ac `Trac #25390: ensure that multivariate factorization factors are monic` ​7b6c93b `Trac #25390: add a test case` ​247f944 `Trac #25390: slight code rearrangement`

### comment:7 Changed 4 years ago by git

• Commit changed from 247f94497f5e9bb646d10cfbf439c1154e983c43 to efdaac1f0611ed3026d5314471b6fad61ca4e869

Branch pushed to git repo; I updated commit sha1. New commits:

 ​efdaac1 `Merge tag '8.3.beta5' into public/25390`

### comment:8 Changed 4 years ago by gh-BrentBaccala

• Dependencies set to #25351
• Status changed from needs_work to needs_review

Still haven't figured out how to access the Singular variable from `libsingular`; got an answer to that SageMath question that made sense, but `nbruin` (the responder) said that it was "tricky" and he was right.

I spent two weeks trying to figure out to make this code work without doing division in polynomial rings over `QQbar` but it's a lot harder than I thought. I'm attaching a PDF that describes the progress I've made and the problems that I've had.

Right now, the code is fairly straightforward, but does require polynomial division over `QQbar`, and therefore requires Trac #25351 to be applied before it works. I've spent enough time trying to break that dependency, so I'm flipping this ticket to `needs_review` and leaving it the way it is. See the attached PDF for more details.

### comment:9 Changed 4 years ago by tscrim

I didn't realize that my suggestion would have been so complicated. Thank you for looking into it. I have two small additional nitpicks, but otherwise this ticket is a positive review assuming #25351 (which I cannot review as that is much more involved mathematically, you should cc some people who work more in that area, such as jdemeyer).

You do not need to create a list:

```norm_f = prod([numfield_f.map_coefficients(h) for h in numfield.embeddings(QQbar)]).change_ring(QQ)
norm_f = prod(numfield_f.map_coefficients(h) for h in numfield.embeddings(QQbar)).change_ring(QQ)
```
```-        REFERENCE::
+        REFERENCES:

-            Geddes, et. al, "Algorithms for Computer Algebra", Section 8.8
+        - Geddes, et. al, "Algorithms for Computer Algebra", Section 8.8
```

although I think it would be better to reference the book in the master ref file and use `- [Geddes1234]_ Section 8.8` as the reference.

### comment:10 Changed 4 years ago by nbruin

If at all possible, you'll probably want to avoid computing norms by multiplying complex embeddings together. Especially when you have your polynomial over a number field, it's probably much better to avoid floats altogether. For instance, if your number field is given as K=QQ[t]/(h(t)) (with h(t) monic) and you have a polynomial f(x) in K[x], then you can represent your polynomial as F(x,t) in QQ[x,t], via the isomorphism K[x] ~ Q[x,t]/(h(t)). Then you have that

Norm_(K[x]/Q[x]) (f) = Resultant(F,h,t)

Code for computing resultants is usually quite optimized. You can experiment a bit to see which approach works best, but I'd expect the resultant to work quite well.

### comment:11 Changed 4 years ago by git

• Commit changed from efdaac1f0611ed3026d5314471b6fad61ca4e869 to cbc51c4c3ff2ffe3cdd726c9479946e07f39f3ce

Branch pushed to git repo; I updated commit sha1. New commits:

 ​3adbd36 `Trac #25390: use resultant to compute polynomial norm` ​9d6353a `Trac #25390: bug fix (if polynomial names conflict with number field names)` ​cbc51c4 `Trac #25390: update references`

### comment:12 Changed 4 years ago by gh-BrentBaccala

I made the changes suggested by tscrim and nbruin; thanks for your feedback.

Using the resultant was more of a pain that I thought it should be. The difficulty lay mainly in needing to introduce a new variable, different from the others. If anybody can suggest a better way that what I came up with, please let me know.

### comment:13 Changed 4 years ago by nbruin

Variable names might be haunting us here:

```sage: k=NumberField(x^2+1,"x")
sage: R=PolynomialRing(k,"x,y")
sage: P=PolynomialRing(QQ,k.gens()+R.gens())
ValueError: variable name 'x' appears more than once
```

Unless you can guarantee the variable names of the number field and the polynomial ring aren't clashing, it seems you cannot go about the problem this way. While it would be reasonable to assume that a user wouldn't intentionally throw this kind of awful naming at you, it's quite conceivable that automatic code would end up doing something like this. Example:

```sage: matrix(k,2,2,[k.0,0,0,1]).charpoly()
x^2 + (-x - 1)*x + x
```

So it would seem the names of your variables for the trivariate ring should be unrelated to the names that occur for the others. Assuming you have a polynomial g in a ring P over a number field k over QQ, I think you can do your conversion with something along the lines of

```R=g.parent()
k=R.base_ring()
P=PolynomialRing(QQ,len(k.gens()+R.gens()),'T')
G=sum(P({(0,)+tuple(k):1})*v.polynomial()(P.0) for k,v in g.dict().iteritems())
NG=G.resultant(k.polynomial()(P.0),P.0)
RQ=PolynomialRing(QQ,R.gens())
NGQ=NG((0,)+RQ.gens())
```

Doing your conversions like this has the advantage that you don't rely on names of generators at all.

### comment:14 Changed 4 years ago by git

• Commit changed from cbc51c4c3ff2ffe3cdd726c9479946e07f39f3ce to 5b23d3f2fdc111d45757db62532ab69f096b124b

Branch pushed to git repo; I updated commit sha1. New commits:

 ​5b23d3f `Trac #25390: improve norm calculation code`

### comment:15 Changed 4 years ago by gh-BrentBaccala

OK, I re-coded it the way nbruin suggested.

### comment:16 Changed 4 years ago by nbruin

```+            norm_f = sum(norm_ring({tuple(k)[1:]:v})
+                         for k,v in norm_flat.dict().iteritems())
```

Please don't think that unpacking a polynomial into a dictionary is the best way of handling it! To "lift" number field elements, I didn't see another way, but in this case

```norm_f = norm_flat((0,)+norm_ring.gens())
```

would work just fine, as far as I can see. If you have a good reason to use the dict upacking (is it faster?) then it's fine, of course; but otherwise using the simpler code is probably preferable.

### comment:17 Changed 4 years ago by git

• Commit changed from 5b23d3f2fdc111d45757db62532ab69f096b124b to a3eb8efc74656da1dbdd05a1762c91e266defa82

Branch pushed to git repo; I updated commit sha1. New commits:

 ​a3eb8ef `Trac #25390: simplify norm calculation code`

### comment:18 Changed 4 years ago by gh-BrentBaccala

I thought there had to a better way to do it!

### comment:20 Changed 4 years ago by git

• Commit changed from a3eb8efc74656da1dbdd05a1762c91e266defa82 to ac30e3a0f45a35bf50c98ee2ea079c5fedd2937e

Branch pushed to git repo; I updated commit sha1. New commits:

 ​ac30e3a `Merge tag '8.4' into public/25390`

### comment:21 Changed 3 years ago by git

• Commit changed from ac30e3a0f45a35bf50c98ee2ea079c5fedd2937e to 94a6c6c80ea91f54417ae983f4e664e58e52e2ee

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

 ​730a40d `Trac #25351: remove @handle_AA_and_QQbar decorator from reduce method` ​a9fde71 `Trac #25351: remove @handle_AA_and_QQbar decorator from _normal_basis_libsingular` ​4f547af `Trac #25351: Put @handle_AA_and_QQbar decorator on groebner_basis(), and remove` ​d187260 `Trac #25351: improve @handle_AA_and_QQbar to handle keyword arguments` ​f83f2bd `Trac #25351: add test cases to methods decorated with @handle_AA_and_QQbar` ​fb145a9 `Trac #25351: modify a test case to check keywords arguments to @handle_AA_and_QQbar` ​758e4ea `Merge tag '8.7.rc0' into public/25351` ​7ba975e `Merge tag '8.7' into public/25351` ​bd8ab1b `Merge tag '8.8.beta2' into public/25351` ​94a6c6c `Merge branch 'public/25351' into public/25390`

### comment:22 Changed 3 years ago by git

• Commit changed from 94a6c6c80ea91f54417ae983f4e664e58e52e2ee to ba45dd6ed636ee5ae903b375cc55362d20cb93c8

Branch pushed to git repo; I updated commit sha1. New commits:

 ​1a3b28c `Merge tag '8.8.beta7' into public/25390` ​7830854 `Trac #25351: ensure Python 3 compatibility` ​ebfbf2d `Merge tag '8.8.beta5' into public/25351` ​b602cb1 `Merge tag '8.8.beta7' into public/25351` ​ba45dd6 `Merge branch 'public/25351' into public/25390`

### comment:23 Changed 3 years ago by git

• Commit changed from ba45dd6ed636ee5ae903b375cc55362d20cb93c8 to f36b1167d0843a62dc17409a7a9f29e410209277

Branch pushed to git repo; I updated commit sha1. New commits:

 ​10554b0 `Trac #25390: fix some doctests` ​f36b116 `Merge tag '8.8.rc0' into public/25390`

### comment:24 Changed 3 years ago by gh-BrentBaccala

• Description modified (diff)

### comment:25 Changed 3 years ago by git

• Commit changed from f36b1167d0843a62dc17409a7a9f29e410209277 to dd4840c3f09656c7bbb64c8b80097f16e46104bd

Branch pushed to git repo; I updated commit sha1. New commits:

 ​64dd01e `Merge tag '8.9.rc1' into public/25390` ​4a955fb `Trac #25390: fix failing doctest` ​3fd38a8 `Trac #25390: remove a duplicate check for _factor_multivariate_polynomial` ​dd4840c `Trac #25390: add a TODO item - maybe use univariate factorization code`

### comment:26 Changed 3 years ago by git

• Commit changed from dd4840c3f09656c7bbb64c8b80097f16e46104bd to a9c1e8b3c0ec2558de2a02ee3930562f08b91ac0

Branch pushed to git repo; I updated commit sha1. New commits:

 ​a9c1e8b `Trac #25390: replace iteritems() with items() for Python 3 compatibility`

### comment:27 Changed 3 years ago by git

• Commit changed from a9c1e8b3c0ec2558de2a02ee3930562f08b91ac0 to 6305cc109f76bcfe3921f8a555973cbd1901a558

Branch pushed to git repo; I updated commit sha1. New commits:

 ​6305cc1 `Merge tag '9.0.beta3' into public/25390`

### comment:28 follow-up: ↓ 30 Changed 3 years ago by tscrim

Two minor things: The first is the `.. TODO::` block should be indented. The second is more stylistic, but I would avoid the blankline after the start of the `for` loops.

### comment:29 Changed 3 years ago by git

• Commit changed from 6305cc109f76bcfe3921f8a555973cbd1901a558 to 91d467ee9ee4008074de452de9ac8a6b0b90359a

Branch pushed to git repo; I updated commit sha1. New commits:

 ​91d467e `Trac #25390: indent TODO block and remove some blank lines`

### comment:30 in reply to: ↑ 28 Changed 3 years ago by gh-BrentBaccala

Two minor things: The first is the `.. TODO::` block should be indented. The second is more stylistic, but I would avoid the blankline after the start of the `for` loops.

OK, done.

I didn't remove the blank lines on the second for loop because it's so long, I feel the blank lines improve it readability.

### comment:31 Changed 3 years ago by tscrim

• Reviewers set to Travis Scrimshaw
• Status changed from needs_review to positive_review

Thanks. I disagree about the space because I have never really seen that elsewhere in Python code and it looks strange to be because of that. Yet, I don't think it is important enough to hold up the ticket.

### comment:32 Changed 3 years ago by chapoton

• Milestone set to sage-9.0

### comment:33 Changed 3 years ago by vbraun

• Branch changed from public/25390 to 91d467ee9ee4008074de452de9ac8a6b0b90359a
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.