New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Root isolation for polynomials with algebraic coefficients #630

Merged
merged 2 commits into from Aug 27, 2018

Conversation

2 participants
@skirpichev
Copy link
Collaborator

skirpichev commented Mar 16, 2018

  • detect cm fields, see sagemath/sage@ba389ae
  • special class for real algebraic fields
  • alter construct_domain to make real algebraic ground domains where needed
  • comparisons without mpmath
  • better representation for gaussian rationals and friends (Complex type?)

see https://groups.google.com/d/msg/sympy/Ujznd13xfgw/LPaqaWnbr40J

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch from 5afcb5f to 5a5dc57 Apr 10, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch from 5a5dc57 to ff62103 Jun 15, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch from c49a8fa to 1b3f4eb Jun 20, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch from 1b3f4eb to e346844 Jun 22, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch 4 times, most recently from 45b4184 to 301b17b Jul 7, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch 7 times, most recently from 841fbf7 to 4e14925 Jul 13, 2018

@skirpichev skirpichev changed the title [wip] Rich comparisons for AlgebraicElement [wip] Root isolation for polynomials with algebraic coefficients Jul 17, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jul 17, 2018

@jksuom, you may want to look on this. If you feel this can be helpful for sympy, I will appreciate review. @cswiercz, this also solves problem in sympy/sympy#9366:

In [1]: aroot = RootOf(x**7 + x + 1, 0, radicals=False)

In [2]: RootOf(y**2 - aroot, 0, radicals=False, extension=True)
Out[2]: 
      ⎛ 2         ⎛ 7           ⎞   ⎞
RootOf⎝y  - RootOf⎝x  + x + 1, 0⎠, 0⎠

In [3]: _.n()
Out[3]: 2.2937127585107e-28 - 0.892493335621313⋅ⅈ

As you can see, there is no changes in rootisolation.py logic. But so far, current code is slow for algebraic domains, take look on the comparison with several rational approximations for coefficients:

In [1]: %time RootOf(x**4 + sqrt(2)*x**3 - I*x + 1, 0, extension=True)
CPU times: user 18 s, sys: 96 ms, total: 18.1 s
Wall time: 18.1 s
Out[1]: 
      ⎛ 4     ___  3             ⎞
RootOf⎝x  + ╲╱ 2 ⋅x  - ⅈ⋅x + 1, 0⎠

In [2]: %time RootOf(x**4 + 46341*x**3/32768 - I*x + 1, 0, extension=True)
CPU times: user 664 ms, sys: 0 ns, total: 664 ms
Wall time: 664 ms
Out[2]: 
      ⎛            3             ⎞
      ⎜ 4   46341⋅x              ⎟
RootOf⎜x  + ──────── - ⅈ⋅x + 1, 0⎟
      ⎝      32768               ⎠

In [3]: %time RootOf(x**4 + 6369051672525773*x**3/4503599627370496 - I*x + 1, 0, extension=True)
CPU times: user 792 ms, sys: 4 ms, total: 796 ms
Wall time: 803 ms
Out[3]: 
      ⎛                       3             ⎞
      ⎜ 4   6369051672525773⋅x              ⎟
RootOf⎜x  + ─────────────────── - ⅈ⋅x + 1, 0⎟
      ⎝       4503599627370496              ⎠

dup_root_upper_bound() wasn't modified yet to return bound for algebraic domains, but it seems to be not the real reason for this slowdown.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Jul 17, 2018

Thanks, this looks interesting. I have only been reading the code, so I cannot see exactly where the time might be spent. There are a couple of points that probably should be investigated. The new _construct_algebraic takes the real and imaginary parts. Will those appear in the unevaluated form (expr + conjugate(expr))/2? Also, what happens when the sign of a real algebraic number is computed?

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jul 17, 2018

so I cannot see exactly where the time might be spent.

I can guess: maybe it's comparisons, or factorization (which is noticeably slower for algebraic domains).

The new _construct_algebraic takes the real and imaginary parts. Will those appear in the unevaluated form (expr + conjugate(expr))/2?

It might be, but I don't consider this as a principal problem. If minimal_polynomial() can handle stuff in this unevaluated form - it's fine:

In [12]: sqrt(2 + root(2 + I, 7))
Out[12]: 
   _______________
  ╱     7 _______ 
╲╱  2 + ╲╱ 2 + ⅈ  

In [13]: minimal_polynomial(_)(x)
Out[13]: 
 28       26        24         22          20          18           16           14           12            10            8           
x   - 28⋅x   + 364⋅x   - 2912⋅x   + 16016⋅x   - 64064⋅x   + 192192⋅x   - 439300⋅x   + 768824⋅x   - 1025360⋅x   + 1026144⋅x  - 747712⋅x

6           4           2        
  + 375424⋅x  - 116480⋅x  + 16901

Also, what happens when the sign of a real algebraic number is computed?

@jksuom, probably, you are about comparisons of real algebraic numbers? You can see how implemented magic methods. They use mpmath capabilites, which may be suboptimal.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Jul 17, 2018

They use mpmath capabilites, which may be suboptimal.

That is a possibility that also occurred to me.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jul 19, 2018

Oh, forgot to provide conjugate()'d example for minimal_polynomial() and expression above:

In [6]: sqrt(2 + root(2 + I, 7))
Out[6]: 
   _______________
  ╱     7 _______ 
╲╱  2 + ╲╱ 2 + ⅈ  

In [7]: minimal_polynomial(_6.conjugate())(x)
---------------------------------------------------------------------------
NotAlgebraic                              Traceback (most recent call last)
<ipython-input-7-58fdbf7062e4> in <module>()
----> 1 minimal_polynomial(_6.conjugate())(x)

~/src/diofant/diofant/polys/numberfields.py in minimal_polynomial(ex, method, **args)
    569                       QQ.frac_field(*ex.free_symbols) if ex.free_symbols else QQ)
    570 
--> 571     result = _minpoly(ex, x, domain)
    572     _, factors = factor_list(result, x, domain=domain)
    573     result = _choose_factor(factors, x, ex, dom=domain)

~/src/diofant/diofant/polys/numberfields.py in _minpoly_compose(ex, x, dom)
    514         res = _minpoly_rootof(ex, x)
    515     else:
--> 516         raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)
    517     return res
    518 

NotAlgebraic: conjugate(sqrt(2 + (2 + I)**(1/7))) doesn't seem to be an algebraic element

Indeed, it currently doesn't work, as you expected. In sympy too. But it should, NotAlgebraic exception - apparently is wrong here. @jksuom, feel free to open an issue in sympy (I can't do this anymore), as it seems there is no such issue opened.

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch 2 times, most recently from 0f2c422 to ac3d72c Aug 11, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 13, 2018

Support for real algebraic domains in rootisolation.py: #673

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch 2 times, most recently from c508d9a to e4549cb Aug 14, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 19, 2018

Change default construction of algebraic domains for polys: #656

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch 5 times, most recently from c236e7e to 2fd9558 Aug 20, 2018

skirpichev added some commits Aug 25, 2018

domains: Add ComplexAlgebraicField subclass
For elements of this domain real and imag properties and
conjugate() method are implemented.
polys: support for complex root isolation of polys with algebraic coe…
…fficients

    In [1]: RootOf(x**4 + sqrt(2)*x**3 - I*x + 1, 0, extension=True)
    Out[1]:
          ⎛ 4     ___  3             ⎞
    RootOf⎝x  + ╲╱ 2 ⋅x  - ⅈ⋅x + 1, 0⎠

    In [2]: _1.interval.as_tuple()
    Out[2]:
    ⎛⎛-283    ⎞  ⎛-283   283⎞⎞
    ⎜⎜─────, 0⎟, ⎜─────, ───⎟⎟
    ⎝⎝ 100    ⎠  ⎝ 200   200⎠⎠

    In [3]: _1.n(3)
    Out[3]: -1.52 + 0.513⋅ⅈ

    In [4]: nroots(_1.poly.as_expr(), 3)[0]
    Out[4]: -1.52 + 0.513⋅ⅈ

Support for isolation of real roots was added in 8e13a2f.

@skirpichev skirpichev force-pushed the skirpichev:new-ANP-2 branch from 2fd9558 to 89a2152 Aug 26, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 26, 2018

Ok, here now is a final version of this PR. ComplexAlgebraicField introduced. I'm not sure this is a long-term variant (I would like to unify somehow RationalField and algebraic numbers), but with this subclass things are much more compact and, I hope, much more clear. I'll appreciate review.

I don't know why root isolation over algebraic fields is so slow (c.f. rational case), but the root algorithm itself seems working as expected, e.g. the dependence of running time on the polynomial degree is cubic.

@skirpichev skirpichev changed the title [wip] Root isolation for polynomials with algebraic coefficients Root isolation for polynomials with algebraic coefficients Aug 26, 2018

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Aug 26, 2018

why root isolation over algebraic fields is so slow

Probably comparisons are slow in real algebraic fields. In principle, root isolation could always be done over rationals, only with polynomials of higher degree. But that might still take less time.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 26, 2018

Probably comparisons are slow in real algebraic fields.

Maybe. I get some speed bonus after reimplementing comparisons to not use mpmath. But that's not only possibility. E.g. factorization is slow over algebraic fields.

In principle, root isolation could always be done over rationals, only with polynomials of higher degree. But that might still take less time.

Mathematica always convert Root's with algebraic coefficients to Root's over integers. Maybe I also could do such expansions much faster, but currently this way seems to be slower in the Diofant:

In [1]: %time RootOf(x**5 + sqrt(2)*x**3 - I*x + 1, 0, extension=True)
CPU times: user 17.3 s, sys: 104 ms, total: 17.4 s
Wall time: 17.5 s
Out[1]: 
      ⎛ 5     ___  3             ⎞
RootOf⎝x  + ╲╱ 2 ⋅x  - ⅈ⋅x + 1, 0⎠

In [2]: %time expand_func(_1)
CPU times: user 3min 16s, sys: 684 ms, total: 3min 17s
Wall time: 3min 27s
Out[2]: 
      ⎛ 20      16      15      12      11      10      8      7      6      5    4      2       ⎞
RootOf⎝x   - 4⋅x   + 4⋅x   + 6⋅x   - 8⋅x   + 6⋅x   + 4⋅x  + 4⋅x  - 4⋅x  + 4⋅x  + x  + 2⋅x  + 1, 2⎠
@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Aug 26, 2018

In [20]: p = Poly(x**5 + sqrt(2)*x**3 - I*x + 1, x, extension=True)
In [23]: %time p.norm()
CPU times: user 2.89 ms, sys: 6.61 ms, total: 9.5 ms
Wall time: 6.91 ms
Out[23]: 
Poly(x**20 - 4*x**16 + 4*x**15 + 6*x**12 - 8*x**11 + 6*x**10 + 4*x**8 + 4*x**7
 - 4*x**6 + 4*x**5 + x**4 + 2*x**2 + 1, x, domain='QQ')
@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 26, 2018

As I already said - resultant computations are fast. Selection of an appropriate root is slow.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Aug 26, 2018

Yes. I was thinking of root isolation only.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 26, 2018

Currently, I use mpmath evaluation in _eval_expand_func helper, but I doubt that I can improve speed of this by order of magnitude, if I replace mpmath with something else. Anyway, we must isolate all roots of resultant like Out[23] (BTW, SymPy sort roots wrongly in this example) and then check roots one-by-one until we found some, which solves original polynomial equation.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Aug 26, 2018

It seems that all roots have to be isolated and checked one-by-one. I was hoping that it would not take such a long time.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Aug 26, 2018

Maybe the process could be accelerated by starting from a floating point approximation of the root in question.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Aug 26, 2018

I was hoping that it would not take such a long time.

Maybe I did something wrongly.

Maybe the process could be accelerated by starting from a floating point approximation of the root in question.

I was thinking about this (using nroots), but I believe it should be the last step in algorithm optimization.

@skirpichev skirpichev merged commit 49cfebb into diofant:master Aug 27, 2018

3 checks passed

codecov/patch 100% of diff hit (target 97%)
Details
codecov/project Absolute coverage decreased by -<1% but relative coverage increased by +2% compared to edf8bda
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@skirpichev skirpichev deleted the skirpichev:new-ANP-2 branch Aug 27, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment