-
-
Notifications
You must be signed in to change notification settings - Fork 656
Faster scalar multiplication for python types #37369
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
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is it possible to have a doctest which shows anything useful?
This is a great catch, thanks for finding and apparently solving the problem. I don't think that I am competent to judge whether the fix is correct and has no unfortunate side effects. Try Nils Bruin?
I thought about the doctest, but the only way to even see what happens is either timing it or doing a trace through the calls made for For the side-effects I am hoping that the CI running the full tests will catch something if it's an issue but I wouldnt be surprised if it doesnt change anything at all as I believe the doctests are parsed by sage before running so no python types will be used across the tests? I'll ping Nils and see what he thinks. I know another type of fix for this would have been to look into the coersion system but as Nils pointed out in devel this is not trivial: https://groups.google.com/g/sage-devel/c/RvmljIs7wBI |
Yes, one would have to do a test using something explicit like comparing 100*P with int(100)*P, but doctests do not seem to be a good way to illustrate a change of this kind, "look how much faster this is". |
Hm, this is a bit strange. The points of an elliptic curve has a secret life as scheme morphism, because a
Can you figure out which code is actually used when I'm not so sure that overriding Another one would be to somehow improve action discovery and allow So ... more concretely for the more limited approach: can you find which two different routines are responsible for computing the multiple of the point in these two cases -- hopefully as much of the call chain as possible. We can then take a look if we should put some special case for EDIT took a bit of a look. I think we get the difference via:
you can now compare I guess we'd be happy if A1 would amount to whatever A2 is doing. Clearly, some action discovery with I think what it comes down to is this line (although I haven't check that that actually gets triggered) sage/src/sage/structure/parent.pyx Lines 2680 to 2683 in 30b3d78
which constructs an It would need some testing and further opinions from other experts on the coercion system, but it seems to me that this would clearly be a spot where the action of |
Looking at the output of the snippet you share above: sage: A1
Left Integer Multiplication by Set of Python objects of class 'int' on Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727
sage: A2
Left action by Integer Ring on Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727 Rather than recasting the This makes me think that: action = detect_element_action(self, S, self_on_left, self_el, S_el) Should not be none for
To coerce all |
well ... the coercion system is asked to figure out an action of an So it's happy to find an "IntegerMultiplicationAction" which it then wraps as by the class I think it's OK to find the ZZ action (if it's there) if you're going to construct an "IntegerMultiplicationAction" anyway. So that was what I was going to construct. This will mean a change in (internal) behaviour in a lot of places, so you'd have to test the change very carefully across the whole library. Be careful with changing the order in which tests are done. Some of the coercion code is carefully constructed to give preference to a fast code path that is often applicable. You don't want to mess that up by doing rare tests earlier on.
I don't think it's Whether
I'd be worried about the cost you might be placing in an otherwise possibly optimized code path. |
My understanding of the codeblock which you shared above: from sage.structure.coerce_actions import detect_element_action
action = detect_element_action(self, S, self_on_left, self_el, S_el)
if action is not None:
return action
if parent_is_integers(S) and not self.has_coerce_map_from(S):
from sage.structure.coerce_actions import IntegerMulAction
try:
return IntegerMulAction(S, self, not self_on_left, self_el)
except TypeError:
_record_exception() Is that the I guess the thing here is that What I don't understand is why Sage would ever want to do something with The coercion system is very complex though and I'm loathed to start ripping it apart without spending lots of time testing and benchmarking the changes. I'll think about it a little but I'm not sure I'll get the best solution quickly... My handwavey thought is I don't know where a check should be made but my gut feeling is that as soon something has a parent of |
Had a little time to test something quick. I added the following to if op is operator.mul: # elements define special action methods.
try:
_register_pair(self, S, "action") # avoid possible infinite loops
# detect actions defined by _rmul_, _lmul_, _act_on_, and _acted_upon_ methods
from sage.structure.coerce_actions import detect_element_action
action = detect_element_action(self, S, self_on_left, self_el, S_el)
if action is not None:
return action
if parent_is_integers(S) and not self.has_coerce_map_from(S):
# New code starts...
# Try the above again, but first coerce integer-like type to Integer
S_el = _Integer(S_el)
S = parent(S_el)
action = detect_element_action(self, S, self_on_left, self_el, S_el)
if action is not None:
return action
# New code ends...
from sage.structure.coerce_actions import IntegerMulAction
try:
return IntegerMulAction(S, self, not self_on_left, self_el)
except TypeError:
_record_exception()
finally:
_unregister_pair(self, S, "action") But this is no good, as the verification of the action fails: ---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[1], line 12
10 xIntp = parent(xInt)
11 xintp = parent(xint)
---> 12 A1 = coercion_model.get_action(xintp,yp,op,xint,y)
File ~/sage/sage/src/sage/structure/coerce.pyx:1723, in sage.structure.coerce.CoercionModel.get_action()
1721 return None
1722
-> 1723 cpdef get_action(self, R, S, op=mul, r=None, s=None) noexcept:
1724 """
1725 Get the action of R on S or S on R associated to the operation op.
File ~/sage/sage/src/sage/structure/coerce.pyx:1760, in sage.structure.coerce.CoercionModel.get_action()
1758 pass
1759 action = self.discover_action(R, S, op, r, s)
-> 1760 action = self.verify_action(action, R, S, op)
1761 self._action_maps.set(R, S, op, action)
1762 return action
File ~/sage/sage/src/sage/structure/coerce.pyx:1813, in sage.structure.coerce.CoercionModel.verify_action()
1811
1812 if action.left_domain() is not R or action.right_domain() is not S:
-> 1813 raise RuntimeError("""There is a BUG in the coercion model:
1814 Action found for R %s S does not have the correct domains
1815 R = %s
RuntimeError: There is a BUG in the coercion model:
Action found for R <built-in function mul> S does not have the correct domains
R = Set of Python objects of class 'int'
S = Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727
(should be Integer Ring, Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727)
action = Left action by Integer Ring on Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727 (<class 'sage.structure.coerce_actions.ActedUponAction'>) So something needs to happen higher up the chain, or the verification needs to be looked at... More to think about. The snippet causing the error is in cpdef get_action(self, R, S, op=mul, r=None, s=None) noexcept:
"""
Get the action of R on S or S on R associated to the operation op.
EXAMPLES::
sage: cm = sage.structure.element.get_coercion_model()
sage: ZZx = ZZ['x']
sage: cm.get_action(ZZx, ZZ, operator.mul)
Right scalar multiplication by Integer Ring
on Univariate Polynomial Ring in x over Integer Ring
sage: cm.get_action(ZZx, QQ, operator.mul)
Right scalar multiplication by Rational Field
on Univariate Polynomial Ring in x over Integer Ring
sage: QQx = QQ['x']
sage: cm.get_action(QQx, int, operator.mul)
Right scalar multiplication by Integer Ring
on Univariate Polynomial Ring in x over Rational Field
with precomposition on right by Native morphism:
From: Set of Python objects of class 'int'
To: Integer Ring
sage: A = cm.get_action(QQx, ZZ, operator.truediv); A
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: x = QQx.gen()
sage: A(x+10, 5)
1/5*x + 2
"""
try:
return self._action_maps.get(R, S, op)
except KeyError:
pass
action = self.discover_action(R, S, op, r, s)
action = self.verify_action(action, R, S, op)
self._action_maps.set(R, S, op, action)
return action EDIT From the examples above, there's this interesting one: sage: QQx = QQ['x']
sage: cm.get_action(QQx, int, operator.mul)
Right scalar multiplication by Integer Ring
on Univariate Polynomial Ring in x over Rational Field
with precomposition on right by Native morphism:
From: Set of Python objects of class 'int'
To: Integer Ring This looks maybe like the kind of thing we want for the elliptic curve case? |
Yes, your trial was what I had in mind too. What we'd like to return is what in our case would be
but that doesn't work because ActedUpAction is strict enough that it wants actual Parent type parents. The |
OK, I started working on an alternative fix to the one proposed above following some other code snippets. I now get the following timings: sage: %timeit ZZ(k) * P
232 µs ± 3.88 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit int(k) * P
229 µs ± 6.15 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each) Looking at the following actions: sage: F = GF(2**127 - 1)
....: E = EllipticCurve(F, [0,6,0,1,0])
....: P = E.random_point()
....: op = operator.mul
....: y = P
....: yp = parent(y)
....: xint = int(3)
....: xInt = 3
....: xintp = parent(x)
....: xIntp = parent(xInt)
....: xintp = parent(xint)
....: A1 = coercion_model.get_action(xintp,yp,op,xint,y)
....: A2 = coercion_model.get_action(xIntp,yp,op,xInt,y)
sage: A1
Left action by Integer Ring on Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727
with precomposition on left by Native morphism:
From: Set of Python objects of class 'int'
To: Integer Ring
sage: A2
Left action by Integer Ring on Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 170141183460469231731687303715884105727 To see how I have done this, see the latest commit. I have no idea if this is best practice but it "looks like" the surrounding code... I will let the CI run to see how it does. |
Ah, nice find of the precomposed action. That should do the trick. My concerns would be for performance (does the composed action add significant overhead?) and the premature definition of the coercion: you can just define that inside the if. Another concern could be that this action might cause memory leaks if it references something strongly that wasn't referenced strongly before. The only data structure at risk for this would be the elliptic curve, so I guess it would be good to test this code with loads of elliptic curves and check they don't pile up in memory. It looks OK to me, but memory management in the face of global dictionaries (as involved for action lookups) is tricky. Finally, but this has more to do with the coercion system discovery process itself, I suspect that if you ask the system to discover the ZZ-action on an elliptic curve by giving it a multiple and a point, the cost might actually depend on n, because some paths in the system end up performing the operation to see if it works. If it does that at all, it should do it with "cheap" elements, but I'm not sure it does: in particular, I think it will do it with the actual arguments provided. |
Hmm, all of the points you make above are valid. For the memory management, I suppose the same statement could be made about the original ZZ action with For the speed, I hope my understanding is correct that this PR should only kick in for pretty rare edge cases so almost all operations will simply be unaffected by my change.
Which definition do you mean? I am happy to rewrite the function? I also do not know if the following snippet is best practices: # Try the above again, but first coerce integer-like type to Integer
# with a connecting coersion
R_el = _Integer(S_el) # Map integer-like to Integer
R = parent(R_el) # Is this the best way to get the Integer parent?
# Compute the coercsion from whatever S is to the Integer class
# This should always work because parent_is_integers(S) is True
connecting = R._internal_coerce_map_from(S) Maybe I could do something smarter here...
Hmm this is an interesting point, but feels orthogonal to what we're trying to fix here? I'm happy to try and think about it, but I'm not sure the best method. There definitely is a measurable change for the sizes though: sage: a = 1
sage: %timeit a * a
36.5 ns ± 0.37 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
sage: a = 2**256
sage: %timeit a * a
49.8 ns ± 0.598 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
sage: a = 2**1024
sage: %timeit a * a
373 ns ± 9.74 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
sage:
sage: a = 1
sage: %timeit a + a
34.9 ns ± 0.899 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
sage: a = 2**256
sage: %timeit a + a
39.1 ns ± 1.11 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
sage: a = 2**1024
sage: %timeit a + a
198 ns ± 3.77 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each) |
For your latest code, what's the top and bottom for? |
Also, the Conda 3.11 CI gives plenty of failing examples. |
I don't know what you're talking about.
Yep! Looks like the CI is doing a good job of showing me how my change here is not good. I guess this coercion is even more delicate than I realised. Now why it's breaking is a bit of a mystery... I'll try and find out tomorrow what I could have done to cause such huge breaks EDIT I tried wrapping the new addition in |
I cannot think of situations in actual python code where an explicit
integer multiplier would be used except for some very small ones, e.g. 2 .
It's surely more likely to arise when some other sage function has returned
an int where it should always return an int: and that is also likely to
only be the case for very small ints.
…On Tue, 20 Feb 2024, 07:30 nbruin, ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/sage/structure/parent.pyx
<#37369 (comment)>:
> @@ -2678,6 +2678,28 @@ cdef class Parent(sage.structure.category_object.CategoryObject):
return action
if parent_is_integers(S) and not self.has_coerce_map_from(S):
+ # Try the above again, but first coerce integer-like type to Integer
+ # with a connecting coersion
+ R_el = _Integer(S_el) # Map integer-like to Integer
+ R = parent(R_el) # Is this the best way to get the Integer parent?
+
+ # Compute the coercsion from whatever S is to the Integer class
+ # This should always work because parent_is_integers(S) is True
+ connecting = R._internal_coerce_map_from(S)
move this line
—
Reply to this email directly, view it on GitHub
<#37369 (review)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAXU32K4PTQVSWEF2LDBVL3YURGIPAVCNFSM6AAAAABDMFZG22VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTQOBZGY3DCNBTGU>
.
You are receiving this because your review was requested.Message ID:
***@***.***>
|
This change seems to have totally broken Sage but I don't understand why at all. I think I need help from someone more experienced with the coercion system... I can't explain why this new action would cause so many issues... |
I found at least one of the bugs, which I have patched with a simple import statement for |
I do the same! But the integers are unlikely to be large if entered literally in your code, surely? |
For cryptography, at least for some of the isogeny-based stuff I write, we often need to clear large scalar factors. E.g. for a point of order 2^n we may want the point of order two. We then write 2**(n-1) * P, this kind of thing is exactly how I discovered this bug when writing some code and benchmarking :) |
OK -- you have convinced me that this could really happen in real life. And has probably happened to me too... |
Hopefully your latest change helps. As far as I can see, sage/src/sage/structure/parent.pyx Lines 134 to 139 in 87cc761
it looks like it's a runtime-delayed import. of So rather than importing
You could get ZZ directly via a similar guarded There are some typos of the word "coercion" in your comments at present; that should be easy to fix. |
Thanks for the feedback and the better way of handling I have also identified the failing doctest which comes from: Ensure (:trac:`27893`) is fixed::
sage: K.<f> = QQ[]
sage: gmpy2.mpz(2) * f
2*f It seems that The current commit simply checks if Does anyone have opinions about this? |
The coercion system is a delicate beast and requires cooperation of the parents involved to work properly (particularly, some coercion info is cached on the parents). I don't think we can generally expect the coercion system to work reliably with non-sage types and I wouldn't mind if we don't raise expectations for that. There are some convenience exceptions, like My concern would be that more common cases would get slowed down if |
Ok, looks like I squashed all the bugs. Now the question is how best to handle this edge case for the |
Oh it looks like we wrote comments at the same time as each other! I think I agree that leaving As such, I feel this PR is ready for a review. After the work we've done together the CI passes, the The timings on my machine support that accidental large sage: F = GF(2**127 - 1)
....: E = EllipticCurve(F, [0,6,0,1,0])
....: P = E.random_point()
....: k = randint(0, 2**127)
sage:
sage: %timeit int(k) * P
228 µs ± 3.5 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit P * int(k)
227 µs ± 5.31 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit ZZ(k) * P
233 µs ± 3.96 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit P * ZZ(k)
235 µs ± 5.79 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each) |
I have been following the thread and think you have done a good job, but I don't really known enough to judge. I think it is OK for @nbruin to do the review. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Excellent job. I think we need a doctest for this, though. Something like
With #37369, registered multiplication actions by `ZZ` are also discovered and used when an `int` is multiplied. Previously it was only discovering the generic Integer Multiplication Action that all additive groups have.
sage: E = EllipticCurve(GF(17),[1,1])
sage: coercion_model.discover action(ZZ, E, operator.mul)
Left action by Integer Ring on Elliptic Curve defined by y^2 = x^3 + x + 1 over Finite Field of size 17
sage: coercion_model.discover action(int, E, operator.mul)
...
Previously this detected
Left Integer Multiplication by Set of Python objects of class 'int' on Elliptic Curve defined by y^2 = x^3 + x + 1 over Finite Field of size 17
and now we should get some precomposed action out of it. As you've seen, the coercion system is a delicate construction where changes can have highly non-local effects. So any test to confirm its behaviour is extremely helpful. If people want to change the system that is fine, but at least they'll have to think about the effects their change has. Unfortunately, the printing of the action hardly tests its actual behaviour, but here it does confirm a different code path taken, at least.
EDIT Actually, it looks like there is already a test that should be adapted:
sage/src/sage/structure/parent.pyx
Lines 2609 to 2614 in 391a721
sage: coercion_model.get_action(E, int, operator.mul) | |
Right Integer Multiplication by Set of Python objects of class 'int' | |
on Elliptic Curve defined by y^2 = x^3 + x over Rational Field | |
sage: coercion_model.get_action(int, E, operator.mul) | |
Left Integer Multiplication by Set of Python objects of class 'int' | |
on Elliptic Curve defined by y^2 = x^3 + x over Rational Field |
Why did that not trigger? With the change here, that test should have failed, shouldn't it?
A bonus would be to have a test where we still do discover a Left IntegerMultiplication action. The following would do:
class P(Parent):
def an_element(self):
return 0
def __repr__(self):
return "P object"
PP = P()
coercion_model.discover_action(int,PP,operator.mul)
(with ZZ
it returns a similar action. I think it would be good to have a test in there that documents when the "add/double" action is actually discovered. I don't know of a ZZ-module implementation somewhere that doesn't implement an explicit ZZ-action or scalar multiplication, so perhaps this pathological shell of a parent is the cheapest option for testing this behaviour)
I believe this is because over the rationals there's no additional action -- only over finite fields is the fast pari action even called? As for the old case, we still have this for sage: type(coercion_model.discover_action(ZZ, E, operator.mul))
<class 'sage.structure.coerce_actions.ActedUponAction'>
sage: type(coercion_model.discover_action(int, E, operator.mul))
<class 'sage.categories.action.PrecomposedAction'>
sage: type(coercion_model.discover_action(gmpy2.mpz, E, operator.mul))
<class 'sage.structure.coerce_actions.IntegerMulAction'> Or alternatively: sage: type(coercion_model.discover_action(ZZ, EllipticCurve(QQ,[1,1]), operator.mul))
<class 'sage.structure.coerce_actions.IntegerMulAction'> |
Ah, thanks! So rational points on elliptic curves over QQ are examples of ZZ-modules without an explicitly declared action by ZZ. That example does just fine in documenting the behaviour then. In that case I think we just want to add the example that should that for elliptic curves over finite fields another action for |
OK -- the latest commit has a doctest explaining the change of this PR and has an example of the difference in actions for ZZ and int. If you think this is OK then I think we are done! As a stupid aside left for someone else, or me when i feel like it, sage: type(coercion_model.discover_action(numpy.int64, E, operator.mul))
<class 'sage.categories.action.PrecomposedAction'>
sage: type(coercion_model.discover_action(gmpy2.mpz, E, operator.mul))
<class 'sage.structure.coerce_actions.IntegerMulAction'> |
Thanks so much for your help with this PR -- coercion is scary! |
I think it's this code: sage/src/sage/rings/integer_ring.pyx Lines 595 to 599 in 30fecca
numpy types are specifically recognized and given a |
Documentation preview for this PR (built with commit 137cb27; changes) is ready! 🎉 |
There is an annoying problem at the moment. If you compute the scalar multiplication of a point on an elliptic curve over a finite field, a fast call to Pari is made when the scalar is a Sage type, but when the scalar is a Python
int
then the multiplication for_acted_upon_
is skipped by the coercion system and__mul__
is called instead which is a much slower Sage implementation using double and add.As an example of this see the following run times.
This is particularly annoying if you write Sagemath code in a
.py
file and write2**100 * P
, which will be parsed such that2**100
is anint
and the very slow method is used. For this 127-bit characteristic example above, the slowdown is more than a factor of 10.This PR writes small methods for
__mul__
and__rmul__
which wrap_acted_upon_
and allow the speed up for Pythonint
values. This results in the following times for the same example: