Domino Domino - 1 month ago 5x
C Question

Converting squared and cube terms into multiplication

I am trying to convert a big expression from sage into valid C code using ccode() from sympy. However, my expression has many squared and cube terms. As pow(x,2) is far slower than x*x, I'm trying to expand those terms in my expression before the conversion.
Based on this conversation, I wrote the following code :

from sympy import Symbol, Mul, Pow, pprint, Matrix, symbols
from sympy.core import numbers

def pow_to_mul(expr):
Convert integer powers in an expression to Muls, like a**2 => a*a.
pows = list(expr.atoms(Pow))
pows = [p for p in pows if p.as_base_exp()[1]>=0]
if any(not e.is_Integer for b, e in (i.as_base_exp() for i in pows)):

raise ValueError("A power contains a non-integer exponent")
repl = zip(pows, (Mul(*[b]*e,evaluate=False) for b,e in (i.as_base_exp() for i in pows)))
return expr.subs(repl)

It partially works, but fails as long as the power is argument of a multiplication:

>>>print pow_to_mul((x^3+2*x^2)._sympy_())
2*x**2 + x*x*x
>>>print pow_to_mul((x^2/(1+x^2)+(1-x^2)/(1+x^2))._sympy_())
x**2/(x*x + 1) - (x*x - 1)/(x*x + 1)

Why? And how can I change that ?
Thank you very much,


If you compile with -ffast-math the compiler will do this optimization for you. If you are using an ancient compiler or cannot affect the level of optimization used in the build process you may pass a user defined function to ccode (using SymPy master branch):

>>> ccode(x**97 + 4*x**7 + 5*x**3 + 3**pi, user_functions={'Pow': [
...   (lambda b, e: e.is_Integer and e < 42, lambda b, e: '*'.join([b]*int(e))),
...   (lambda b, e: not e.is_Integer, 'pow')]})
'pow(x, 97) + 4*x*x*x*x*x*x*x + 5*x*x*x + pow(3, M_PI)'