generic log/exp for univariate polynomials

[videlec]

We implement generic _log_series and _exp_series on univariate polynomials.

Depends on #28618

CC: @mezzarobba

Component: algebra

Author: Vincent Delecroix

Branch/Commit: u/vdelecroix/28592 @ 1475df1

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

[videlec]

Branch: u/vdelecroix/28592

[videlec]

New commits:

6c69313

generic log and exp for univariate polynomials

[videlec]

Commit: 6c69313

[fchapoton]

comment:2

What is the advantage compared to the existing

sage: R.<x> = QQ[[]]
sage: (1 + x).O(4).log()
x - 1/2*x^2 + 1/3*x^3 + O(x^4)

?

[mezzarobba]

comment:3

@Frédéric: The advantage is that it allows us to have specialized implementations for specific coefficient rings using the existing polynomial classes. I think the idea is that, eventually, operations on elements of power series rings should become wrappers for these _XXX_series() methods. It also makes it possible to call truncated operations without the overhead of constructing a power series when you have a polynomial.

@Vincent: From a quick look, the code looks fine to me, but isn't it a bit of a pity to have a quadratic-time implementation while there are simple quasi-linear algorithms? (Integrate f'/f to compute log(f), Newton-invert to get exp(f).) I'm not saying you need to change it, but should you want to do it, here is an implementation done in a small student project I supervised where initial goal was to submit the code to Sage, but the students didn't do it and I never took the time to clean up the code myself.

[videlec]

comment:4

Thanks Mark for the pointer. I just quickly implemented that because I needed it! I should rather backport what is in power series (ODE integration).

[videlec]

comment:5

The solve_linear_ode (from power series) is limited to degree one differential equations and hence does not handle cos and sin. This is kind of unfortunate...

[sagetrac-git]

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

b4407ac

28592: generic log and exp for polynomials

[sagetrac-git]

Changed commit from 6c69313 to b4407ac

[mezzarobba]

comment:8

Thank you for the implementation!

• A small bug:

sage: S.<t> = SR[]
....: S(3)._log_series(1)
0

• A failing test:

File "src/sage/rings/polynomial/polynomial_element.pyx", line 10228, in sage.rings.polynomial.polynomial_element.Polynomial._log_series
Failed example:
    x._log_series()
Expected:
    Traceback (most recent call last):
    ...
    ValueError: can not take logarithm with zero constant coefficient
Got:
    <BLANKLINE>
    Traceback (most recent call last):
      File "/home/marc/co/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 681, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/marc/co/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 1123, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.rings.polynomial.polynomial_element.Polynomial._log_series[10]>", line 1, in <module>
        x._log_series()
    TypeError: _log_series() takes exactly one argument (0 given

• self = self[1:] would be better than self -= a0 over RLF, intervals, ...
• I'm not sure if testing the characteristic instead of just letting the code fail with a division by zero is actually a good idea, because there may be applications where one wants to compute truncated series to a precision < characteristic. But that's a really minor detail, and largely a matter of taste.

Replying to @videlec:

The solve_linear_ode (from power series) is limited to degree one differential equations and hence does not handle cos and sin. This is kind of unfortunate...

FWIW, one way to compute sin/cos without code for solving ODEs of order two is via 1/(1+t²) → arctan → tan → tan²/2.

[videlec]

comment:9

Replying to @mezzarobba:

Thanks for the feedback!

• self = self[1:] would be better than self -= a0 over RLF, intervals, ...

I have a question about this. There is the class function shift(self, n) that is supposed to do the same thing as self[n:] (when n is positive). However, the shift is very often not optimized and goes through lists and calling the parent. Is there a prefered way (for readability/uniformity of the code)?

• I'm not sure if testing the characteristic instead of just letting the code fail with a division by zero is actually a good idea, because there may be applications where one wants to compute truncated series to a precision < characteristic. But that's a really minor detail, and largely a matter of taste.

Fine. It will make the code faster :-)

Replying to @videlec:

The solve_linear_ode (from power series) is limited to degree one differential equations and hence does not handle cos and sin. This is kind of unfortunate...

FWIW, one way to compute sin/cos without code for solving ODEs of order two is via 1/(1+t²) → arctan → tan → tan²/2.

This gives you arctan and tan. But how do you extract cos/sin from there?

[videlec]

comment:10

The implementation fails for non-commutative base rings

sage: M = MatrixSpace(ZZ,2)
sage: x = polygen(M)
sage: p = 1 + M([2,1,0,4]) * x + M([1,0,-1,3]) * x**2 + M([1,1,1,1]) * x**3
sage: p._log_series(10)._exp_series(10) == p
False

Is there a way to fix it or should I forbid non-commutative base rings?

[videlec]

comment:11

Replying to @videlec:

Replying to @mezzarobba:

• self = self[1:] would be better than self -= a0 over RLF, intervals, ...

This is not the same operation. There is a shift in the degree with self[1:].

[videlec]

comment:12

Replying to @videlec:

The implementation fails for non-commutative base rings

sage: M = MatrixSpace(ZZ,2)
sage: x = polygen(M)
sage: p = 1 + M([2,1,0,4]) * x + M([1,0,-1,3]) * x**2 + M([1,1,1,1]) * x**3
sage: p._log_series(10)._exp_series(10) == p
False

Is there a way to fix it or should I forbid non-commutative base rings?

Note that the naive series expansion that I initially wrote work for non-commutative rings!

[videlec]

comment:13

Replying to @videlec:

Replying to @videlec:

The implementation fails for non-commutative base rings

sage: M = MatrixSpace(ZZ,2)
sage: x = polygen(M)
sage: p = 1 + M([2,1,0,4]) * x + M([1,0,-1,3]) * x**2 + M([1,1,1,1]) * x**3
sage: p._log_series(10)._exp_series(10) == p
False

Is there a way to fix it or should I forbid non-commutative base rings?

Note that the naive series expansion that I initially wrote work for non-commutative rings!

... modulo #28610

[videlec]

comment:14

Last question: the following is kind of wrong

sage: x = polygen(ZZ)
sage: (1+x)._log_series(5)
- 1/4*x^4 + 1/3*x^3 - 1/2*x^2 + x

as the answer belongs to QQ[x] and not ZZ[x]. The reason is that the integral method silently makes a coercion to the fraction field. Should integral be more careful with promotion of base rings?

[sagetrac-git]

Changed commit from b4407ac to 1475df1

[sagetrac-git]

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

1b4bf57	28592: polynomial ring promotion for integral division
bd177ed	allow non-positive precisions in _mul_trunc_
1475df1	29592: improved log/exp

[videlec]

comment:16

Needs review again! Now one can take logarithm over MatrixSpace(Zmod(77), 2)['x'].

[videlec]

comment:19

Boils down to

sage: cm.get_action(QQ['x'], ZZ, operator.truediv)
Right inverse action by Rational Field on Univariate Polynomial Ring in x over Rational Field
with precomposition on right by Natural morphism:
  From: Integer Ring
  To:   Rational Field
sage: cm.get_action(Qp(7)['x'], ZZ, operator.truediv)

[videlec]

comment:20

which boils down to

sage: cm.get_action(QQ['x'], ZZ, operator.mul)
Right scalar multiplication by Integer Ring on Univariate Polynomial Ring in x over Rational Field
sage: cm.get_action(Qp(7)['x'], ZZ, operator.mul)

[mezzarobba]

comment:21

Replying to @videlec:

FWIW, one way to compute sin/cos without code for solving ODEs of order two is via 1/(1+t²) → arctan → tan → tan²/2.

This gives you arctan and tan. But how do you extract cos/sin from there?

I meant tan(z/2) in the last step, sorry; and then sin and cos using the rational parametrisation of the circle.

[mezzarobba]

comment:22

Replying to @videlec:

Replying to @videlec:

Replying to @mezzarobba:

• self = self[1:] would be better than self -= a0 over RLF, intervals, ...

This is not the same operation. There is a shift in the degree with self[1:].

Is there?

sage: P.<x> = QQ[]
sage: (1 + x)[1:]
x

[mezzarobba]

comment:23

Replying to @videlec:

The implementation fails for non-commutative base rings

Hmm, right, I hadn't thought of that.

Is there a way to fix it or should I forbid non-commutative base rings?

I don't know, sorry.

[videlec]

Dependencies: #28618

[videlec]

comment:24

Replying to @videlec:

grrr, there is a bug with p-adic polynomial rings and the coercion model

sage: cm = get_coercion_model()
sage: cm.bin_op(QQ['x'].one(), ZZ.one(), operator.truediv).parent()
Univariate Polynomial Ring in x ...
sage: cm.bin_op(Qp(7)['x'].one(), ZZ.one(), operator.truediv).parent()
Fraction Field of Univariate Polynomial Ring in x ...

Opened #28618...

[embray]

comment:25

Ticket retargeted after milestone closed

[mkoeppe]

comment:26

Batch modifying tickets that will likely not be ready for 9.1, based on a review of the ticket title, branch/review status, and last modification date.

[mkoeppe]

comment:28

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[mkoeppe]

comment:29

Setting a new milestone for this ticket based on a cursory review.