Opened 7 years ago

Closed 6 years ago

Last modified 6 years ago

#18282 closed enhancement (fixed)

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

Reported by: mmezzarobba Owned by:
Priority: major Milestone: sage-6.9
Component: commutative algebra Keywords:
Cc: rws, jpflori Merged in:
Authors: Marc Mezzarobba Reviewers: Ralf Stephan
Report Upstream: N/A Work issues:
Branch: 8ed1967 (Commits, GitHub, GitLab) Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by mmezzarobba)

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.

Change History (21)

comment:1 Changed 7 years ago by git

  • Commit changed from 64339bcbe4c00d54447e76c20862f93ea92f77c1 to 4fdccb24c57c1886aa8849bb1f8e6fe8c4cae1ea

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

012631eHorner form of symbolic expressions
4fdccb2Improve pol(expr) for expr in SR

comment:2 Changed 7 years ago by mmezzarobba

Last edited 7 years ago by mmezzarobba (previous) (diff)

comment:3 Changed 7 years ago by mmezzarobba

  • Description modified (diff)
  • Summary changed from Improve the evaluation of polynomials on symbolic expressions to Improve the default evaluation method for univariate polynomials

comment:4 Changed 7 years ago by rws

  • Component changed from symbolics to commutative algebra

comment:5 Changed 7 years ago by git

  • Commit changed from 4fdccb24c57c1886aa8849bb1f8e6fe8c4cae1ea to 1f7feb6ad9527fb1a237294d68a0a1fcb08484da

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

e62e761Horner form of symbolic expressions
4d59c7dImprove pol(expr) for expr in SR
1f7feb6Simplify Polynomial.__call__

comment:6 Changed 7 years ago by git

  • Commit changed from 1f7feb6ad9527fb1a237294d68a0a1fcb08484da to 89671731fdf80e827f11b120a6655659c28c0225

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

8967173Simplify Polynomial.__call__

comment:7 Changed 7 years ago by mmezzarobba

  • Status changed from new to needs_review

comment:8 follow-up: Changed 7 years ago by 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? In which ring? You've evaluated x, so the answer should lie in QQy?, 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.

comment:9 in reply to: ↑ 8 ; follow-up: Changed 7 years ago by mmezzarobba

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 QQy?, 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.

comment:10 in reply to: ↑ 9 Changed 7 years ago by nbruin

Replying to mmezzarobba:

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.

comment:11 Changed 6 years ago by jpflori

  • Cc jpflori added

comment:12 Changed 6 years ago by mmezzarobba

  • Branch changed from u/mmezzarobba/wip/pol_of_symb to u/mmezzarobba/pol_call
  • Commit changed from 89671731fdf80e827f11b120a6655659c28c0225 to 0cbf0e5395205c323f5a06bc0f9502bf52a88b8c
  • Milestone changed from sage-6.7 to sage-6.8
Last edited 6 years ago by mmezzarobba (previous) (diff)

comment:13 Changed 6 years ago by git

  • Commit changed from 0cbf0e5395205c323f5a06bc0f9502bf52a88b8c to 84d6a1b874c241e8eed8196cc42fd9075be40206

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__

comment:14 Changed 6 years ago by mmezzarobba

  • Summary changed from Improve the default evaluation method for univariate polynomials to Fixes, cleanup and improvements to the default evaluation method for univariate polynomials

Rebased.

comment:15 follow-up: Changed 6 years ago by bruno

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^2^500, at least on "easy" values such as 1 or -1 or on finite fields.

What do you think of my proposal?

comment:16 in reply to: ↑ 15 ; follow-up: Changed 6 years ago by mmezzarobba

Replying to bruno:

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!

comment:17 in reply to: ↑ 16 Changed 6 years ago by bruno

Replying to mmezzarobba:

Replying to bruno:

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.

comment:18 Changed 6 years ago by rws

  • Branch changed from u/mmezzarobba/pol_call to public/18282

comment:19 Changed 6 years ago by rws

  • Commit changed from 84d6a1b874c241e8eed8196cc42fd9075be40206 to 8ed19674f4788e97c9f3697fdb67a06c2ef2c781
  • Milestone changed from sage-6.8 to sage-6.9
  • Reviewers set to Ralf Stephan
  • Status changed from needs_review to positive_review

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:

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

comment:20 Changed 6 years ago by vbraun

  • Branch changed from public/18282 to 8ed19674f4788e97c9f3697fdb67a06c2ef2c781
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:21 Changed 6 years ago by mmezzarobba

  • Commit 8ed19674f4788e97c9f3697fdb67a06c2ef2c781 deleted

Thanks!

Note: See TracTickets for help on using tickets.