RoofOf Modular Exponentiation

[Chris Swierczewski]

Hello,

Has work has been done on modular reduction of `RootOf` objects? (If not, I would like to begin implementing this functionality!) What I mean is this:

import sympy
from sympy.abc import x
f = x**3 + 2*x + 1
a = sympy.RootOf(f,0)
b = a**3

Since $f(a) = 0$ by definition it would be nice to return

print b
-1 - 2*a

instead of

print b
a**3

so that expression swell could be curbed.

If this hasn't been implemented could anyone give suggestions on what tools to keep in mind when overloading `RoofOf.__pow__()` and related operator methods? (I'm proficient at Python and OOP but don't have much experience with CAS design. Even if this sounds like a bad idea for Sympy I'd like to try it out on my own `MyRootOf` class.

--

Chris

[Ondřej Čertík]

On Wed, Nov 26, 2014 at 11:31 AM, Chris Swierczewski <cswi...@gmail.com> wrote:
> Hello,
>
> Has work has been done on modular reduction of `RootOf` objects? (If not, I
> would like to begin implementing this functionality!) What I mean is this:
>
> import sympy
> from sympy.abc import x
> f = x**3 + 2*x + 1
> a = sympy.RootOf(f,0)
> b = a**3
>
> Since $f(a) = 0$ by definition it would be nice to return
>
> print b
> -1 - 2*a
>
> instead of
>
> print b
> a**3
>
> so that expression swell could be curbed.

This is what I get:

In [1]: f = x**3 + 2*x + 1

In [2]: a = RootOf(f, 0)

In [3]: b = a**3

In [4]: a
Out[4]:
⎛ 3 ⎞
RootOf⎝x + 2⋅x + 1, 0⎠

In [5]: b
Out[5]:
3
⎛ 3 ⎞
RootOf⎝x + 2⋅x + 1, 0⎠

So looks like it hasn't been implemented yet.

>
> If this hasn't been implemented could anyone give suggestions on what tools
> to keep in mind when overloading `RoofOf.__pow__()` and related operator
> methods? (I'm proficient at Python and OOP but don't have much experience
> with CAS design. Even if this sounds like a bad idea for Sympy I'd like to
> try it out on my own `MyRootOf` class.

Good question. First, I would implement RootOf.mypow(n) and see if you
can figure out how to get the answer you want. That's the hard part.
What's the algorithm to determine which powers simplify and which not?

Then the next question is how to hook this up into sympy, so that it
simplifies automatically, I guess one would need to add an if
statement into the Pow simplification somewhere.

Ondrej

[Chris Swierczewski]

Hello,

Thanks for the quick response!

> This is what I get:

> ...

> So looks like it hasn't been implemented yet.

Right on. I'll start poking around, then.

> Good question. First, I would implement RootOf.mypow(n) and see if you
> can figure out how to get the answer you want. That's the hard part.
> What's the algorithm to determine which powers simplify and which not?

I was considering starting with a binary method like the one outlined here:

http://en.wikipedia.org/wiki/Modular_exponentiation

(I remember this from a number theory course I took with William Stein back in my undergraduate days!) I guess the simplification convention, if there is one, is to reduce modulo the leading term of the polynomial. So if `alpha` is a root of

f(x) = c x**n + lower_order_terms(x)

then

alpha**n = lower_order_terms(alpha) / c

Any high power of the root can be given in terms of lower powers. The harder part will be making this fast. (To my knowledge the binary method is the most common fast approach.)

> Then the next question is how to hook this up into sympy, so that it
> simplifies automatically, I guess one would need to add an if
> statement into the Pow simplification somewhere.

The more I think about it the more it looks like the bulk of the problem is to implement modular exponentiation of polynomials...

I'll implement a "mypow" first and see what people have to say about making this play nice with everything.

--
Chris Swierczewski
University of Washington
Department of Applied Mathematics
www.cswiercz.info

[Ondřej Čertík]

Hi Chris,

On Wed, Nov 26, 2014 at 11:56 AM, Chris Swierczewski <cswi...@gmail.com> wrote:
> Hello,
>
> Thanks for the quick response!
>
>> This is what I get:
>> ...
>> So looks like it hasn't been implemented yet.
>
> Right on. I'll start poking around, then.
>
>> Good question. First, I would implement RootOf.mypow(n) and see if you
>> can figure out how to get the answer you want. That's the hard part.
>> What's the algorithm to determine which powers simplify and which not?
>
> I was considering starting with a binary method like the one outlined here:
>
> http://en.wikipedia.org/wiki/Modular_exponentiation
>
> (I remember this from a number theory course I took with William Stein back in my undergraduate days!) I guess the simplification convention, if there is one, is to reduce modulo the leading term of the polynomial. So if `alpha` is a root of
>
> f(x) = c x**n + lower_order_terms(x)
>
> then
>
> alpha**n = lower_order_terms(alpha) / c
>
> Any high power of the root can be given in terms of lower powers. The harder part will be making this fast. (To my knowledge the binary method is the most common fast approach.)

I see, yes, that should work. Clever.

>
>> Then the next question is how to hook this up into sympy, so that it
>> simplifies automatically, I guess one would need to add an if
>> statement into the Pow simplification somewhere.
>
> The more I think about it the more it looks like the bulk of the problem is to implement modular exponentiation of polynomials...
>
> I'll implement a "mypow" first and see what people have to say about making this play nice with everything.

Exactly, that's always my approach to first implement the algorithm
and then worry how to best hook it up in sympy (or any other
software).

Let us know if you run into any problems. Thanks for working on this,
that would be a great addition to RootOf.

Ondrej

[Aaron Meurer]

There's a decent chance this algorithm is already implemented and just
not integrated into the RootOf object, so I would search around the
polys code first.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVARTrP6t5KucJ3zjJMhpONdCTirfZmat5-s8uq7CgniFg%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Chris Swierczewski]

Hello Aaron,

> There's a decent chance this algorithm is already implemented and just
> not integrated into the RootOf object, so I would search around the
> polys code first.

Thanks for the tip. I've been getting to know the `rem()` and `reduce()` functions a little more. From what I've seen so far I still think it'll require some work since the point of efficient modular exponentiation is to avoid computing `p(x)**n` first before reducing it modulo `x**p`. Nonetheless, it'll be a good first step.

--
Chris

[Chris Swierczewski]

I have a preliminary implementation below that I can write up as a pull request if people would prefer to take a look that way. I have not performed any of the sympy tests yet but will do so later this afternoon.

class MyRootOf(RootOf):

    __doc__ = RootOf.__doc__

    

    def __init__(self, *args, **kwds):

        RootOf.__init__(self, *args, **kwds)  

        

    def _eval_power(self, exponent):

        if exponent >= self.poly.degree():

            exponent = int(exponent)

            result = 1

            base = self.poly.gen

            modulus = self.poly

            while exponent > 0:

                if (exponent % 2):

                    result = reduced(result*base, [modulus])[1]

                exponent >>= 1

                base = reduced(base*base, [modulus])[1]        

            return result.subs(self.poly.gen,self).as_expr()

        return RootOf._eval_power(self, exponent)

This is a basic implementation of a modular exponentiation algorithm. The basic idea is that since a root $b$ satisfies $p(b) = c*b**n = q(b)$ then $b**n = - (1/c) * q(b)$. Here is a simple test:

f = 3*y**5 + 2*y - 1

root = MyRootOf(f,0)

print root**5

-2*RootOf(3*y**5 + 2*y - 1, 0)/3 + 1/3

From what I gathered the method _eval_power() is called every time a Sympy Expr subtype is exponentiated. This makes it so that p.eval(b) or p.subs(b) uses the above modular exponentiation algorithm. Just overloading __pow__() doesn't cut it.

I have some preliminary timings as well demonstrating that expansions of polynomial expressions in this new RootOf is about x1.5 - x2 faster than the original behavior which keeps track of all the higher powers.

I tried finding a way to make expressions like 

root**4 * (1 + root)

symplify appropriately. But after diving into how Sympy handle's Mul's I see that the default behavior of expressions like x*(1+x) is to not automatically expand. I doubt that this behavior should be changed since it'll impact so many other computations. However, if somebody know of how I can default expand expressions involving RootOf / MyRootOf objects I would greatly appreciate the feedback and suggestions.

[Aaron Meurer]

I wouldn't worry about the fact that things don't expand
automatically. This sort of thing is intentional. If users want things
to expand, they will call expand(), or simplify() (which will probably
end up calling expand()), and things will simplify then thanks to the
_eval_power routine.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/bf63ac31-569a-4bcb-b9fd-ae6662d8c6e2%40googlegroups.com.

[Chris Swierczewski]

I wouldn't worry about the fact that things don't expand
automatically. This sort of thing is intentional. If users want things
to expand, they will call expand(), or simplify() (which will probably
end up calling expand()), and things will simplify then thanks to the
_eval_power routine.

Thanks for the confirmation.

Is there a collection of methods which are called when certain actions are performed on an expression? That is, I eventually discovered that _eval_power() is called whenever the object is exponentiated. Are there similar methods for other operations?

Also, for anyone reading this who is interested, I created an issue tracking this new feature here: https://github.com/sympy/sympy/issues/8543

--

Chris

[Aaron Meurer]

On Mon, Dec 1, 2014 at 3:43 PM, Chris Swierczewski <cswi...@gmail.com> wrote:
>
>> I wouldn't worry about the fact that things don't expand
>> automatically. This sort of thing is intentional. If users want things
>> to expand, they will call expand(), or simplify() (which will probably
>> end up calling expand()), and things will simplify then thanks to the
>> _eval_power routine.
>
>
> Thanks for the confirmation.
>
> Is there a collection of methods which are called when certain actions are
> performed on an expression? That is, I eventually discovered that
> _eval_power() is called whenever the object is exponentiated. Are there
> similar methods for other operations?

Not all operations have it, but those that do are usually named _eval_something.

Aaron Meurer

>
> Also, for anyone reading this who is interested, I created an issue tracking
> this new feature here: https://github.com/sympy/sympy/issues/8543
>
> --
> Chris
>

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/CACx6PThHLwA5uj3qOf-xhvrW1bbiP6QWvmJXpjhROp-5dHDojA%40mail.gmail.com.