Domino - 4 months ago 12

C Question

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:

`>>>_=var('x')`

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

Answer

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