SymPy cannot lambdify Product

Question

I'm using SymPy 1.0 and Python 2.7. I want to compute the sum of first 100 integer numbers:

This code runs succesfully

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Sum(x[i], (i, 0, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s, 'numpy')
s_lambda(np.arange(101))

And gives 5050 as expected. But when I try to do the same with a Product instead of a Sum:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 0, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s, 'numpy')
s_lambda(np.arange(101))

I got a NameError: global name 'Product' is not defined What am I doing wrong? Is there a workaround to get what I want?

Edit 1: And what if I don't know in advance the limit of the Product. Let's say something like

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 0, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n) s.doit(), 'numpy')
s_lambda(np.arange(101), 5)

Edit 2: I'm trying to find a workaround. NameError: global name 'Product' is not defined error is given because of this:

lambdastr((sy.DeferredVector('x'), n), p)

That gives:

lambda x,n: (Product(x[i], (i, 0, n)))

While for the Sum we got a working lambda function:

lambda x,n: ((builtins.sum(x[i] for i in range(0, n+1))))

At this point the problem revolves around the definition of the Product function. According to the manual I can inject via a dict my definition of a function

def my_prod(a, b):
    # my implementation
    pass

my_fun = {"Product" : my_prod}
f = sy.lambdify((sy.DeferredVector('x'), n), p, modules=['numpy', my_fun])
f([1,2,3,4,5], 2)

Problem is, list indices must be integers, not Symbol error is raised when I try to use the lambdified function f. I guess this is due to i that is a symbol while it is supposed to be an integer. I can't understand why it's not passed the actual integer before trying to call my_prod as it is for the Sum case.

Answer 1

When the limits of the Product are known in advance

You could work around the problem by calling .doit() to expand the Product into its component parts:

In [104]: s = sy.Product(x[i], (i, 1, 10)); s
Out[104]: Product(x[i], (i, 1, 10))

In [105]: s.doit()
Out[105]: x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[10]

---

For example,

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 1, 10))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(11)))

prints

3628800

---

However, if you use .doit() with sy.Product(x[i], (i, 1, 100)) then you'll get an arithmetic overflow since np.arange(101) has dtype int32 or int64 (depending on your OS) and the product 100!

In [109]: math.factorial(100)
Out[109]: 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

is too large be stored in either an int32 or int64 array value.

In [118]: np.iinfo('int64').max
Out[118]: 9223372036854775807

In [119]: np.iinfo('int64').max < math.factorial(100)
Out[119]: True

Thus,

s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101)))

raises a

RuntimeWarning: overflow encountered in long_scalars

and erroneously prints 0.

---

If you change the input from an array of dtype int64 to a list of Python ints then the product can be computed correctly:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101).tolist()))

prints

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

---

When the limits of the Product are not known in advance

The work-around (AFAICS) becomes more complicated. If you use a debugger to trace the code path followed when Sum is used you'll find that LambdaPrinter._print_Sum is called to convert Sum(x[i], (i, 0, n)) to the expression builtins.sum(x[i] for i in range(0, n+1)).

If we add a _print_Product method to NumPyPrinter (a subclass of LambdaPrinter), then we can get lambdify to successfully convert Product into an expression that NumPy can evaluate:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np
import sympy.printing.lambdarepr as SPL

def _print_Product(self, expr):
    loops = (
        'for {i} in range({a}, {b}+1)'.format(
            i=self._print(i),
            a=self._print(a),
            b=self._print(b))
        for i, a, b in expr.limits)
    return '(prod([{function} {loops}]))'.format(
        function=self._print(expr.function),
        loops=' '.join(loops))
SPL.NumPyPrinter._print_Product = _print_Product

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 1, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n), s, 'numpy')
print(s_lambda(np.arange(101), 5))

prints

120