Skip to content

Conversation

GiacomoPope
Copy link
Contributor

@GiacomoPope GiacomoPope commented Feb 16, 2024

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.

┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.2, Release Date: 2023-12-03                    │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
sage: F = GF(2**127 - 1)
sage: E = EllipticCurve(F, [0,6,0,1,0])
sage: P = E.random_point()
sage: k = randint(0, 2**127)
sage: 
sage: %timeit int(k) * P
3.77 ms ± 98.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage: %timeit P * int(k)
3.74 ms ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage: %timeit ZZ(k) * P
235 µs ± 2.62 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit P * ZZ(k)
234 µs ± 4.81 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

This is particularly annoying if you write Sagemath code in a .py file and write
2**100 * P, which will be parsed such that 2**100 is an int 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 Python int values. This results in the following times for the same example:

┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.3.beta8, Release Date: 2024-02-13              │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable.     ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: F = GF(2**127 - 1)
sage: E = EllipticCurve(F, [0,6,0,1,0])
sage: P = E.random_point()
sage: k = randint(0, 2**127)
sage: 
sage: %timeit int(k) * P
222 µs ± 1.6 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit P * int(k)
225 µs ± 4.54 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit ZZ(k) * P
219 µs ± 1.88 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit P * ZZ(k)
230 µs ± 15.6 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Copy link
Member

@JohnCremona JohnCremona left a 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?

@GiacomoPope
Copy link
Contributor Author

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 int(k) * P versus ZZ(k) * P. I am happy to do additional work to include something if anyone has an idea though.

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

@GiacomoPope GiacomoPope requested a review from nbruin February 16, 2024 16:22
@JohnCremona
Copy link
Member

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 int(k) * P versus ZZ(k) * P. I am happy to do additional work to include something if anyone has an idea though.

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".

@nbruin
Copy link
Contributor

nbruin commented Feb 19, 2024

Hm, this is a bit strange. The points of an elliptic curve has a secret life as scheme morphism, because a k-point amounts to a map Spec(k) -> E. So the multiplication routine __mul__ that you end up in is the one for scheme morphisms. It looks like it tries to preempt the coercion system if the other operand is also a scheme morphism and punts to coercion_model.binop if not. The two lines below seem to lead to the same difference in timing:

sage: %timeit coercion_model.bin_op(int(10^8),P,operator.mul)
sage: %timeit coercion_model.bin_op(10^8,P,operator.mul)

Can you figure out which code is actually used when int is used? Obviously there is an elliptic addition routine through which this is routed and how that differs from what is discovered with Integer? I'm thinking that the code for int might just benefit from punting to gp when appropriate, just as for Integer.

I'm not so sure that overriding __mul__ is the right thing to do here: elliptic curve points are supposed to be scheme points. I don't think anyone has actually used that these are scheme morphisms, but perhaps we shouldn't break that (and your proposed code would break that).

Another one would be to somehow improve action discovery and allow int there as well. We could take a cue from how it's done for multiplication, where discovery does work with python ints as well as Integers. Perhaps actions deserve that too, given that integers act on any abelian group ...

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 int in there somewhere to speed it up.

EDIT took a bit of a look. I think we get the difference via:

sage: y=P
sage: yp=parent(y)
sage: xint=int(3)
sage: xInt=3
sage: xintp=parent(x)
sage: xIntp=parent(xInt)
sage: xintp=parent(xint)
sage: A1=coercion_model.get_action(xintp,yp,op,xint,y)
sage: A2=coercion_model.get_action(xIntp,yp,op,xInt,y)

you can now compare A1._act_ and A2._act_. The code for A2._act_ indeed calls pari. The code for A1._act_ does not. It is very generic code. It actually has a doctest that checks elliptic curve point multiplication in a way that wouldn't trigger this code path!

I guess we'd be happy if A1 would amount to whatever A2 is doing. Clearly, some action discovery with int is possible. So perhaps we should find where and how A1 is discovered and let it discover something more specific for elliptic curve points there?

I think what it comes down to is this line (although I haven't check that that actually gets triggered)

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)

which constructs an IntegerMulAction object. It looks like this is a last-resort object for an integer mul action. In our case there is an action of ZZ on our object and it can be found. It just cannot be found from parent(int(3)). So I think that code there should actually first retry detect_element_action but now with ZZ as parent and possibly the element coerced in there. We have at this point already checked that parent_is_integers so the conversion should be safe.

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 int gets identified in a rudimentary way.

@GiacomoPope
Copy link
Contributor Author

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 int to and Integer and requesting the action again, I feel like Left Integer Multiplication by Set of Python objects of class 'int' on... is something which a user would never want? Surely we would always want Left action by Integer Ring on Abelian group of points... even if one accidentally uses int rather than Int (here the fact python int exists at all is usually because coercion has been missed?

This makes me think that:

action = detect_element_action(self, S, self_on_left, self_el, S_el)

Should not be none for Integer or int in this case, which means patching the function from

from sage.structure.coerce_actions import detect_element_action

To coerce all int to Integer?

@nbruin
Copy link
Contributor

nbruin commented Feb 19, 2024

Rather than recasting the int to and Integer and requesting the action again, I feel like Left Integer Multiplication by Set of Python objects of class 'int' on... is something which a user would never want? Surely we would always want Left action by Integer Ring on Abelian group of points... even if one accidentally uses int rather than Int (here the fact python int exists at all is usually because coercion has been missed?

well ... the coercion system is asked to figure out an action of an int on a sage object. Such actions are not discovered for int (I think because only actions between Parents are discovered). The system does detect that the int object is sufficiently integer-like to try and see if there is an IntegerMultiplicationAction, which exists as soon as addition is supported (which it actually detects by shudder trying to see if it can add an element. There is a disturbing amount of detection in the coercion framework based on an_element, where it just tries if it can do something with an element. I seriously disagree with that design decision, but I don't think I have the resources to fix that)

So it's happy to find an "IntegerMultiplicationAction" which it then wraps as by the class int (that's just some wrapping -- whether the parent on that is registered as int or Integers() will not make much difference), but it misses the fact that there is actually an explicit action of ZZ declared, which does something more efficient that basic add/double.

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.

This makes me think that:

action = detect_element_action(self, S, self_on_left, self_el, S_el)

Should not be none for Integer or int in this case, which means patching the function from

I don't think it's none for ZZ. I think that's how the right action is found. I think action detection for non-parent types automatically fails. Note that only later does the code test if the element is integer-like. I think that's considered a rare case, or at least one for which the code path doesn't need to be optimized.

Whether detect_element_action can find stuff for python-types ... you'd have to dig into why that was originally explicitly excluded an reevaluate.

To coerce all int to Integer?

I'd be worried about the cost you might be placing in an otherwise possibly optimized code path.

@GiacomoPope
Copy link
Contributor Author

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 detect_element_action() is called and returns something not None when we multiply by an Integer but this is None for an int which causes the fall through to the attempt to do IntegerMulAction which succeeds in the case for int.

I guess the thing here is that detect_element_action() is searching for whether there's a sensible _acted_upon_ method which will work. For the case of scalar multiplication for elliptic curves, the _acted_upon_ method then does ZZ(input) anyway but the issue with the int is that this is never reached.

What I don't understand is why Sage would ever want to do something with Integer which it couldn't do with an int (or at least coerce to Integer from int).

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 int should just always be considered Integer as soon as possible. I know when working in .sage files, or in the shell this is essentially automatic, but having weird slowdowns like this one when int comes through either via a call to randint() rather than ZZ.random_element() or when just writing a number in a .py file is not ideal...

I don't know where a check should be made but my gut feeling is that as soon something has a parent of <class 'int'>, we should be coercing this to ZZ.

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Feb 19, 2024

Had a little time to test something quick. I added the following to discover_action()

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 coerce.pyx:

    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?

@nbruin
Copy link
Contributor

nbruin commented Feb 19, 2024

Yes, your trial was what I had in mind too. What we'd like to return is what in our case would be

sage.structure.coerce_actions.ActedUponAction(int,E,True)

but that doesn't work because ActedUpAction is strict enough that it wants actual Parent type parents. IntegerMulAction shows things apparently don't break badly if there is an action object that does have int as an "acting" parent.

The IntegerAction objects do the wrapping already. It's just that the IntegerMulAction doesn't examine the parent sufficiently to see if another action is more appropriate. Hence, I think that the init of IntegerMulAction should basically do a detect_action(ZZ,...) and if that is successful, configure the object to use S._acted_upon rather than its own add/double implementation.

@GiacomoPope
Copy link
Contributor Author

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.

@nbruin
Copy link
Contributor

nbruin commented Feb 19, 2024

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.

@GiacomoPope
Copy link
Contributor Author

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 Integer on many curves? I can do some basic tests though. Not sure if what I add here will make problems which dont already exist.

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.

and the premature definition of the coercion: you can just define that inside the if.

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...

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 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)

@grhkm21
Copy link
Contributor

grhkm21 commented Feb 19, 2024

For your latest code, what's the top and bottom for?

@grhkm21
Copy link
Contributor

grhkm21 commented Feb 19, 2024

Also, the Conda 3.11 CI gives plenty of failing examples.

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Feb 19, 2024

what's the top and bottom for?

I don't know what you're talking about.

Also, the Conda 3.11 CI gives plenty of failing examples.

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 if connecting is not None: just in case there was some weird case here, either this will fix a mystery issue or the CI will fail again!

@JohnCremona
Copy link
Member

JohnCremona commented Feb 20, 2024 via email

@GiacomoPope
Copy link
Contributor Author

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...

@GiacomoPope
Copy link
Contributor Author

I found at least one of the bugs, which I have patched with a simple import statement for ZZ but I think this is not what I should do generally. If the CI tests pass with this change then I will start to work at optimising what I have written.

@JohnCremona
Copy link
Member

@JohnCremona I write my Sage scripts in .py files and often forget to wrap everything in ZZ. I don’t know if I am in the very minority

I do the same! But the integers are unlikely to be large if entered literally in your code, surely?

@GiacomoPope
Copy link
Contributor Author

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 :)

@JohnCremona
Copy link
Member

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...

@nbruin
Copy link
Contributor

nbruin commented Feb 20, 2024

Hopefully your latest change helps. As far as I can see, _Integer is defined here:

cdef object _Integer
cdef bint is_Integer(x) noexcept:
global _Integer
if _Integer is None:
from sage.rings.integer import Integer as _Integer
return type(x) is _Integer or type(x) is int

it looks like it's a runtime-delayed import. of sage.rings.integer.Integer, which is the type of sage integers; not its parent (which is ZZ). There is probably a good reason for doing the delayed import and there may be a good reason to not just do the from ... import all the time -- it is an operation that does have quite a bit of overhead, even if the module is already imported.

So rather than importing ZZ and using that, I'd probably stick closer to the example here and do:

    global _Integer
    if _Integer is None:
        from sage.rings.integer import Integer as _Integer
    ZZ_el = _Integer(S_el)
    ZZ = ZZ_el.parent()

You could get ZZ directly via a similar guarded from .. import into a cache instead. I don't know if the construction through the type is more efficient/sager than the conversion call on ZZ. I would expect that the conversion could trigger more coercion shenanigans, but perhaps it doesn't.

There are some typos of the word "coercion" in your comments at present; that should be easy to fix.

@GiacomoPope
Copy link
Contributor Author

Thanks for the feedback and the better way of handling _Integer. I have fixed this, as well as the typos in the latest commit.

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 ZZ._internal_coerce_map_from(S) is None for the case where S is type gmpy2.mpz. As the connecting map is None things fail.

The current commit simply checks if connecting is None and I hope this fixes the last of the failure, but I dont know if I like this hack. Another option would be to look at ZZ._internal_coerce_map_from and make this work for gmpy2.mpz in the same way it works for int.

Does anyone have opinions about this?

@nbruin
Copy link
Contributor

nbruin commented Feb 20, 2024

The current commit simply checks if connecting is None and I hope this fixes the last of the failure, but I dont know if I like this hack. Another option would be to look at ZZ._internal_coerce_map_from and make this work for gmpy2.mpz in the same way it works for int.

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 int, because they can accidentally end up as input with interactive use. Beyond that, however, I don't think we need to go out of our way, so if gmpy2.mpz doesn't work properly, I think that's fine. We have gmp wrapped in Integer already -- and it's a pretty lightweight wrapper.

My concern would be that more common cases would get slowed down if ZZ._internal_coerce_map_from gets adjusted. It can probably be done without cost for the common code path but you'd probably have to check that's the case.

@GiacomoPope
Copy link
Contributor Author

Ok, looks like I squashed all the bugs. Now the question is how best to handle this edge case for the gmpy2.int class, whether we do as above and simply skip it, or try including this as an option for ZZ._internal_coerce_map_from() I'll do some experiments.

@GiacomoPope
Copy link
Contributor Author

Oh it looks like we wrote comments at the same time as each other!

I think I agree that leaving gmpy2.mpz using the slower fall back to avoid a slowdown in _internal_coerce_map_from is worth it.

As such, I feel this PR is ready for a review. After the work we've done together the CI passes, the ell_point.py is unchanged and the changes to parent.pyx are minimal.

The timings on my machine support that accidental large int will not cause slowdown:

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)

@JohnCremona
Copy link
Member

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.

Copy link
Contributor

@nbruin nbruin left a 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: 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)

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Feb 21, 2024

Why did that not trigger? With the change here, that test should have failed, shouldn't it?

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 gmpy2.mpz. We could for example have something along the lines of:

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'>

@nbruin
Copy link
Contributor

nbruin commented Feb 21, 2024

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 int is discovered, that is in accordance with the action for ZZ (although we can't really check the latter fact -- only that int gives a precomposed action).

@GiacomoPope
Copy link
Contributor Author

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, gmpy2.mpz really is a sore thumb of the lot. Even numpy ints have a nice coercion:

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'>

@GiacomoPope
Copy link
Contributor Author

Thanks so much for your help with this PR -- coercion is scary!

@nbruin
Copy link
Contributor

nbruin commented Feb 22, 2024

As a stupid aside left for someone else, or me when i feel like it, gmpy2.mpz really is a sore thumb of the lot. Even numpy ints have a nice coercion:

I think it's this code:

elif is_numpy_type(S):
import numpy
if issubclass(S, numpy.integer):
return True
return None

numpy types are specifically recognized and given a _coerce_map_from result. The problem with an mpz is of course that it can be very efficiently be coerced into an Integer (just bitstring copy or possibly even by reference -- although that might mess up memory management). Whatever the coercion system comes up with will be quite inefficient. But it's pretty clear you could get it to work; probablty just by testing for the type here. Given that ZZ(gmpy2.mpz(100)) just works, I think there's not an obstacle. It will mainly just invite stupid code, though, so not supporting it might actually be better for that purpose.

Copy link

Documentation preview for this PR (built with commit 137cb27; changes) is ready! 🎉

@vbraun vbraun merged commit 17a8ebd into sagemath:develop Mar 31, 2024
@GiacomoPope GiacomoPope deleted the faster_mul branch April 1, 2024 01:41
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants