factorize result of algdep()

[jdemeyer]

Instead of returning composite polynomials from algdep(), factor the result and return the best factor.

Depends on #22715

CC: @dimpase

Component: basic arithmetic

Author: Jeroen Demeyer

Branch/Commit: 65ab385

Reviewer: Moritz Firsching

Issue created by migration from https://trac.sagemath.org/ticket/22714

[jdemeyer]

Dependencies: #22715

[jdemeyer]

Branch: u/jdemeyer/factorize_result_of_algdep__

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

2f6afb2	Define `__abs__` for p-adic numbers
65ab385	Factorize result of algdep()

[sagetrac-git]

Commit: 65ab385

[mo271]

comment:6

This seems like a very nice enhancement! I have had a very similar setup in my own scripts for identifying algebraic numbers...

I plan to run some tests and review it swiftly.

[mo271]

Reviewer: Moritz Firsching

[mo271]

comment:7

I tried with this polynomial:

sage: x = polygen(AA)
sage: p = x^71 - x^69 - 2*x^68 - x^67 + 2*x^66 + 2*x^65 + x^64 - x^63 - x^62 - x^61 - x^60 - x^59 + 2*x^58 + 5*x^57 + 3*x^56 - 2*x^55 - 10*x^54 - 3*x^53 - 2*x^52 + 6*x^51 + 6*x^50 + x^49 + 9*x^48 - 3*x^47 - 7*x^
....: 46 - 8*x^45 - 8*x^44 + 10*x^43 + 6*x^42 + 8*x^41 - 5*x^40 - 12*x^39 + 7*x^38 - 7*x^37 + 7*x^36 + x^35 - 3*x^34 + 10*x^33 + x^32 - 6*x^31 - 2*x^30 - 10*x^29 - 3*x^28 + 2*x^27 + 9*x^26 - 3*x^25 + 14*x^24 - 8
....: *x^23 - 7*x^21 + 9*x^20 + 3*x^19 - 4*x^18 - 10*x^17 - 7*x^16 + 12*x^15 + 7*x^14 + 2*x^13 - 12*x^12 - 4*x^11 - 2*x^10 + 5*x^9 + x^7 - 7*x^6 + 7*x^5 - 4*x^4 + 12*x^3 - 6*x^2 + 3*x - 6
sage: t = AA.polynomial_root(p, RIF(1.3,1.4)).n(digits=400)
sage: algdep(t, 77)
x^71 - x^69 - 2*x^68 - x^67 + 2*x^66 + 2*x^65 + x^64 - x^63 - x^62 - x^61 - x^60 - x^59 + 2*x^58 + 5*x^57 + 3*x^56 - 2*x^55 - 10*x^54 - 3*x^53 - 2*x^52 + 6*x^51 + 6*x^50 + x^49 + 9*x^48 - 3*x^47 - 7*x^46 - 8*x^45 - 8*x^44 + 10*x^43 + 6*x^42 + 8*x^41 - 5*x^40 - 12*x^39 + 7*x^38 - 7*x^37 + 7*x^36 + x^35 - 3*x^34 + 10*x^33 + x^32 - 6*x^31 - 2*x^30 - 10*x^29 - 3*x^28 + 2*x^27 + 9*x^26 - 3*x^25 + 14*x^24 - 8*x^23 - 7*x^21 + 9*x^20 + 3*x^19 - 4*x^18 - 10*x^17 - 7*x^16 + 12*x^15 + 7*x^14 + 2*x^13 - 12*x^12 - 4*x^11 - 2*x^10 + 5*x^9 + x^7 - 7*x^6 + 7*x^5 - 4*x^4 + 12*x^3 - 6*x^2 + 3*x - 6
sage: f = gp.algdep(t, 77); f
x^74 - x^72 - x^71 - x^70 + x^69 + x^66 + x^65 - 2*x^63 - 2*x^62 + x^61 + 4*x^60 + 2*x^59 - 5*x^57 - 4*x^55 - 4*x^54 + 3*x^53 - x^52 + 15*x^51 + 3*x^50 - 6*x^49 + x^48 - 11*x^47 + 3*x^46 - 2*x^45 + 5*x^43 - 6*x^42 + 15*x^41 - 12*x^40 - 5*x^39 + 8*x^38 - 10*x^37 + 17*x^36 + 2*x^35 - 9*x^34 + 8*x^33 - 9*x^32 - 9*x^31 - x^29 - 6*x^28 + 16*x^27 + x^26 - 3*x^25 + 7*x^24 + x^23 + 3*x^22 - 11*x^21 - x^20 - 4*x^19 + 8*x^18 - 3*x^17 - 5*x^16 + 3*x^14 - 7*x^12 - 4*x^11 - x^10 - 2*x^9 + 7*x^8 - 3*x^7 + 5*x^6 + x^5 - x^4 + 6*x^3 - 6*x^2 + 3*x - 6
sage: factor(f)
[x + 1, 1; x^2 - x + 1, 1; x^71 - x^69 - 2*x^68 - x^67 + 2*x^66 + 2*x^65 + x^64 - x^63 - x^62 - x^61 - x^60 - x^59 + 2*x^58 + 5*x^57 + 3*x^56 - 2*x^55 - 10*x^54 - 3*x^53 - 2*x^52 + 6*x^51 + 6*x^50 + x^49 + 9*x^48 - 3*x^47 - 7*x^46 - 8*x^45 - 8*x^44 + 10*x^43 + 6*x^42 + 8*x^41 - 5*x^40 - 12*x^39 + 7*x^38 - 7*x^37 + 7*x^36 + x^35 - 3*x^34 + 10*x^33 + x^32 - 6*x^31 - 2*x^30 - 10*x^29 - 3*x^28 + 2*x^27 + 9*x^26 - 3*x^25 + 14*x^24 - 8*x^23 - 7*x^21 + 9*x^20 + 3*x^19 - 4*x^18 - 10*x^17 - 7*x^16 + 12*x^15 + 7*x^14 + 2*x^13 - 12*x^12 - 4*x^11 - 2*x^10 + 5*x^9 + x^7 - 7*x^6 + 7*x^5 - 4*x^4 + 12*x^3 - 6*x^2 + 3*x - 6, 1]

Another possibility might have been to introduce an option irreducible=True or something simimlar, but I cannot imagine a situation, where one is interested in the smaller factors. Also the lines

    factors = [p for p, e in R(f).factor()]
    return min(factors, key=lambda f: abs(f(z)))

don't slow down the computation in a any significant way.

Great enhancement!

[dimpase]

Work Issues: fix complex_double

[dimpase]

comment:8

I do not understand why this does not touch src/sage/rings/complex_double.pyx, which has its own (unfixed) algdep (calling PARI directly), and needless to say it still fails on FreeBSD 11.0 with clang:

sage -t --warn-long 41.0 src/sage/rings/complex_double.pyx
**********************************************************************
File "src/sage/rings/complex_double.pyx", line 2363, in sage.rings.complex_double.ComplexDoubleElement.algdep
Failed example:
    p = z.algdep(5); p
Expected:
    x^3 + 1
Got:
    x^5 + x^2
**********************************************************************
File "src/sage/rings/complex_double.pyx", line 2365, in sage.rings.complex_double.ComplexDoubleElement.algdep
Failed example:
    p.factor()
Expected:
    (x + 1) * (x^2 - x + 1)
Got:
    (x + 1) * x^2 * (x^2 - x + 1)
**********************************************************************

[jdemeyer]

comment:9

Replying to @dimpase:

I do not understand why this does not touch src/sage/rings/complex_double.pyx, which has its own (unfixed) algdep (calling PARI directly), and needless to say it still fails on FreeBSD 11.0 with clang:

I created #22759 for that.

[jdemeyer]

Changed work issues from fix complex_double to none

[vbraun]

Changed branch from u/jdemeyer/factorize_result_of_algdep__ to 65ab385