Fixes, cleanup and improvements to the default evaluation method for univariate polynomials

[mezzarobba]

Simplify Polynomial.__call__(), fix the following issues (note that in calls of the form p(x,y,...), x is the outermost variable), and add a few tests.

sage: Pol_x.<x> = QQ[]
sage: Pol_xy.<y> = Pol_x[]
sage: pol = 1000*x^2*y^2 + 100*y + 10*x + 1
sage: pol(y, 0)
1000*x^2*y^2 + 100*y + 10*x + 1

sage: pol(y, 0)
1000*x^2*y^2 + 100*y + 10*x + 1

sage: pol(~y, 0) # not the same bug as above
((10*x + 1)*y^2 + 100*y + 1000*x^2)/y^2

sage: pol(x, y, x=1)
1000*y^2 + 10*y + 101

sage: zero = Pol_xy(0)
sage: zero(1).parent()
Integer Ring

sage: zero = QQ['x'](0)
sage: a = matrix(ZZ, [[1]])
sage: zero(a).parent() # should be over QQ
Full MatrixSpace of 1 by 1 dense matrices over Integer Ring

sage: zero = GF(2)['x'](0)
sage: zero(1.).parent() # should raise an error
Real Field with 53 bits of precision

sage: pol(y=x, x=1)
1111

sage: pol = QQ['x'](range(10))
sage: pol(x) # technically not a bug, but should be expanded
((((((((9*x + 8)*x + 7)*x + 6)*x + 5)*x + 4)*x + 3)*x + 2)*x + 1)*x

Also implement a method to compute the Horner form of a polynomial expression, in order not to lose the “feature” illustrated in the last example.

CC: @rwst @jpflori

Component: commutative algebra

Author: Marc Mezzarobba

Branch: 8ed1967

Reviewer: Ralf Stephan

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

[sagetrac-git]

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

012631e	Horner form of symbolic expressions
4fdccb2	Improve pol(expr) for expr in SR

[sagetrac-git]

Changed commit from 64339bc to 4fdccb2

[sagetrac-git]

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

e62e761	Horner form of symbolic expressions
4d59c7d	Improve pol(expr) for expr in SR
1f7feb6	Simplify Polynomial.__call__

[sagetrac-git]

Changed commit from 4fdccb2 to 1f7feb6

[sagetrac-git]

Changed commit from 1f7feb6 to 8967173

[sagetrac-git]

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

8967173

Simplify Polynomial.__call__

[nbruin]

comment:8

sage: pol(y, 0)
1000*x^2*y^2 + 100*y + 10*x + 1

If I understand correctly, you now make this evaluate to 100*y+1? In which ring? You've evaluated x, so the answer should lie in QQ['y'], which is not a ring that exists at this moment. What about pol(y,I)? should that return an answer in QuadraticField(-1,name="I")['y']?

-1 on this: The method you're calling is evaluating a univariate polynomial over an arbitrary ring. You don't know what "evaluation" of the coefficient would mean and, more importantly, whether it's supported at all. Furthermore, the rings in which the answer is supposed to lie likely don't exist yet. Letting sage choose which rings should be constructed likely leads to difficult to predict behaviour (you'd basically be relying on the common parents the coercion framework cooks up, and then I think it's better to let the user rely on it him/herself).

I think that pol(y,0) should be an error because there's an unhandled coefficient 0 present. It's not clear at all that the second argument should be used for the variable that gets mentioned first in QQ['x']['y']. One might think that pol(y0,x0)==pol(y0)(x0), but that's not the case either. Error really is safer.

[mezzarobba]

comment:9

Replying to @nbruin:

sage: pol(y, 0)
1000*x^2*y^2 + 100*y + 10*x + 1

If I understand correctly, you now make this evaluate to 100*y+1?

Yes:

sage: pol(y, 0) # with patch
100*y + 1

Note that the previous implementation wouldn't work in the case of pol(y, 0), but would happily compute

sage: sage: pol(y+1, 0) # without patch
100*y + 101

In which ring? You've evaluated x, so the answer should lie in QQ['y'], which is not a ring that exists at this moment.

Not exactly:

sage: pol(y, 0).parent() # with patch
Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Rational Field

The logic here is that yes, I have evaluated x, so the evaluated coefficients lie in ℚ, but then I'm evaluating the resulting element of ℚ[Y] on y ∈ ℚ[x][y], not y ∈ ℚ[y]. So the answer should lie in ℚ[x][y].

In contrast,

sage: pol(0, x).parent() # with patch
Univariate Polynomial Ring in x over Rational Field

since x.parent() is QQ[x], not QQ[x][y].

What about pol(y,I)? should that return an answer in QuadraticField(-1,name="I")['y']?

sage: pol(y, I).parent() # with patch - TBI, see #18036
Symbolic Ring
sage: pol(y, I.pyobject()).parent() # with patch
Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Number Field in I with defining polynomial x^2 + 1

Furthermore, the rings in which the answer is supposed to lie likely don't exist yet. Letting sage choose which rings should be constructed likely leads to difficult to predict behaviour (you'd basically be relying on the common parents the coercion framework cooks up,

Well, yes, but that's already the case when you just do a + b! Would you really want the evaluation of elements of ℤ[x][y] on y = y0 ∈ ℚ to raise an error?

I think that pol(y,0) should be an error because there's an unhandled coefficient 0 present. It's not clear at all that the second argument should be used for the variable that gets mentioned first in QQ['x']['y']. One might think that pol(y0,x0)==pol(y0)(x0), but that's not the case either. Error really is safer.

Perhaps, yes, but I didn't invent this feature. It has been present for years, and people use it! There are even examples in the sage library that rely on pol(y,x) working when pol ∈ R[y] where R is not a polynomial ring (but another ring with callable elements). So really this ticket is only about making the implementation understandable, and fixing lots of corner cases such as those mentioned above.

[nbruin]

comment:10

Replying to @mezzarobba:

Perhaps, yes, but I didn't invent this feature. It has been present for years, and people use it! There are even examples in the sage library that rely on pol(y,x) working when pol ∈ R[y] where R is not a polynomial ring (but another ring with callable elements). So really this ticket is only about making the implementation understandable, and fixing lots of corner cases such as those mentioned above.

OK, that seems to be the case indeed. Thanks for cleaning things up a bit.

[mezzarobba]

Changed commit from 8967173 to 0cbf0e5

[mezzarobba]

Changed branch from u/mmezzarobba/wip/pol_of_symb to u/mmezzarobba/pol_call

[sagetrac-git]

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

5ffd321	#18282 Horner form of symbolic expressions
4a519fd	#18282 Improve pol(expr) for expr in SR
84d6a1b	#18282 Simplify Polynomial.__call__

[sagetrac-git]

Changed commit from 0cbf0e5 to 84d6a1b

[mezzarobba]

comment:14

Rebased.

[bgrenet]

comment:15

In view of #18600 (and older sisters #18518 and #18585), I think __call__ may be at the same time improved and adapted to high degree sparse polynomials: Basically, when a univariate polynomial p is evaluated on a value v (other than a symbolic variable), it uses a Horner scheme that ranges over all the coefficients of p, including the zeroes. One could easily modify lines 755-756 in src/sage/rings/polynomial/polynomial_element.pyx to make the loop ranges over the nonzero coefficients. I've not checked already whether it is slower in a sensible way for really dense polynomials or not, but it would be faster for sparse polynomials. And as a side benefit, one would be able to evaluate polynomials such as x<sup>2</sup>500, at least on "easy" values such as 1 or -1 or on finite fields.

What do you think of my proposal?

[mezzarobba]

comment:16

Replying to @bgrenet:

One could easily modify lines 755-756 in src/sage/rings/polynomial/polynomial_element.pyx to make the loop ranges over the nonzero coefficients. I've not checked already whether it is slower in a sensible way for really dense polynomials or not, but it would be faster for sparse polynomials.

Sounds reasonable—or perhaps one should override __call__ in Polynomial_generic_sparse, I don't know. In any case I have no time to spend on this ticket now, but please feel free to add improvements if you want!

[bgrenet]

comment:17

Replying to @mezzarobba:

Replying to @bgrenet:

One could easily modify lines 755-756 in src/sage/rings/polynomial/polynomial_element.pyx to make the loop ranges over the nonzero coefficients. I've not checked already whether it is slower in a sensible way for really dense polynomials or not, but it would be faster for sparse polynomials.

Sounds reasonable—or perhaps one should override __call__ in Polynomial_generic_sparse, I don't know. In any case I have no time to spend on this ticket now, but please feel free to add improvements if you want!

Actually, I did some tests. It appears that it should be better to implement a __call__ method in Polynomial_generic_sparse to avoid hindering performances. The other solution would be to have tests inside the current __call__ method to check whether the parent is sparse or not. One advantage is to avoid code duplication (for *args and *kwds), though I guess it is not the right solution.

I'd better let this ticket as it is (and actually try to review it...) and open a new ticket for sparse polynomials.

[rwst]

Changed branch from u/mmezzarobba/pol_call to public/18282

[rwst]

Changed commit from 84d6a1b to 8ed1967

[rwst]

Reviewer: Ralf Stephan

[rwst]

comment:19

In Pynac-0.4.3 (maybe 0.3.9.3) sin(pi/5) expands immediately to 1/4*sqrt(-2*sqrt(5) + 10) (actually every sin/cos/tan value expressible with sqrts of depth 3), so I'll change the doctest to sin(pi/7).

As I trust Nils on the general purpose, and I can see nothing missing in the code, also the patchbot is happy and my patch is only a doctest change I'll take the liberty to set positive.

---

New commits:

8ed1967

Merge branch 'u/mmezzarobba/pol_call' of trac.sagemath.org:sage into tmp4

[vbraun]

Changed branch from public/18282 to 8ed1967

[mezzarobba]

comment:21

Thanks!

[mezzarobba]

Changed commit from 8ed1967 to none