Opened 4 years ago

Closed 4 years ago

#22714 closed enhancement (fixed)

factorize result of algdep()

Reported by: jdemeyer Owned by:
Priority: major Milestone: sage-8.0
Component: basic arithmetic Keywords:
Cc: dimpase Merged in:
Authors: Jeroen Demeyer Reviewers: Moritz Firsching
Report Upstream: N/A Work issues:
Branch: 65ab385 (Commits) Commit: 65ab385bcc4f202920c0f1b23db32f77a4f495d8
Dependencies: #22715 Stopgaps:

Description

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

Change History (11)

comment:1 Changed 4 years ago by jdemeyer

  • Cc dimpase added

comment:2 Changed 4 years ago by jdemeyer

  • Dependencies set to #22715

comment:3 Changed 4 years ago by jdemeyer

  • Branch set to u/jdemeyer/factorize_result_of_algdep__

comment:4 Changed 4 years ago by git

  • Commit set to 65ab385bcc4f202920c0f1b23db32f77a4f495d8

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

2f6afb2Define __abs__ for p-adic numbers
65ab385Factorize result of algdep()

comment:5 Changed 4 years ago by jdemeyer

  • Status changed from new to needs_review

comment:6 Changed 4 years ago by moritz

  • Reviewers set to Moritz Firsching

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.

comment:7 Changed 4 years ago by moritz

  • Status changed from needs_review to positive_review

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!

comment:8 follow-up: Changed 4 years ago by dimpase

  • Status changed from positive_review to needs_work
  • Work issues set to fix complex_double

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)
**********************************************************************

comment:9 in reply to: ↑ 8 Changed 4 years ago by jdemeyer

  • Status changed from needs_work to positive_review

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.

comment:10 Changed 4 years ago by jdemeyer

  • Work issues fix complex_double deleted

comment:11 Changed 4 years ago by vbraun

  • Branch changed from u/jdemeyer/factorize_result_of_algdep__ to 65ab385bcc4f202920c0f1b23db32f77a4f495d8
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.