优化sympy生成的代码

Question

我在使用SymPy来计算导数的时候，参考了这个问题：https://math.stackexchange.com/questions/726104/apply-chain-rule-to-vector-function-with-chained-dot-and-cross-product，写了以下代码：

from sympy import *
from sympy.physics.mechanics import *
from sympy.printing import print_ccode
from sympy.utilities.codegen import codegen


x1, x2, x3 = symbols('x1 x2 x3')
y1, y2, y3 = symbols('y1 y2 y3')
z1, z2, z3 = symbols('z1 z2 z3')

u = ReferenceFrame('u')

u1=(u.x*x1 + u.y*y1 + u.z*z1)
u2=(u.x*x2 + u.y*y2 + u.z*z2)
u3=(u.x*x3 + u.y*y3 + u.z*z3)

s1=(u1-u2).normalize()
s2=(u2-u3).normalize()
v=cross(s1, s2)
f=dot(v,v)

df_dy2=diff(f, y2)


print_ccode(df_dy2, assign_to='df_dy2')


[(c_name, c_code), (h_name, c_header)] = codegen( ("df_dy2", df_dy2), "C", "test", header=False, empty=False)

print c_code

运行后得到了这个结果：

#include "test.h"
#include <math.h>
double df_dy2(double x1, double x2, double x3, double y1, double y2, double y3, double z1, double z2, double z3) {
   return ((x1 - x2)*(y2 - y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - (x2 - x3)*(y1 - y2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(2*(x1 - x2)*(y1 - y2)*(y2 - y3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(x1 - x2)*(-y2 + y3)*(y2 - y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x1 - x2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x2 - x3)*pow(y1 - y2, 2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x2 - x3)*(y1 - y2)*(-y2 + y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x2 - x3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2)))) + (-(x1 - x2)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + (x2 - x3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(-2*(x1 - x2)*(y1 - y2)*(z2 - z3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x1 - x2)*(-y2 + y3)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x2 - x3)*(y1 - y2)*(z1 - z2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(x2 - x3)*(-y2 + y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L))) + ((y1 - y2)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - (y2 - y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(2*pow(y1 - y2, 2)*(z2 - z3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(y1 - y2)*(-y2 + y3)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) - 2*(y1 - y2)*(y2 - y3)*(z1 - z2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(-y2 + y3)*(y2 - y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) - 2*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))));
}

结果中有很多相同数字的平方根和幂的重复计算，这样会影响代码的可读性和执行速度。不过我不知道该怎么处理这些重复的计算……

问题1：你知道有没有办法让sympy自动处理这些重复计算吗？

问题2：你知道有没有其他工具可以对这段代码进行后处理吗？

问题3：gcc能在编译时优化这个吗？为什么？

Answer 1

这是我自己写的小脚本，灵感来自asmeurers的提示：

def sympyToC( symname, symfunc ):
    tmpsyms = numbered_symbols("tmp")
    symbols, simple = cse(symfunc, symbols=tmpsyms)
    symbolslist = map(lambda x:str(x), list(symfunc.atoms(Symbol)) )
    symbolslist.sort()
    varstring=",".join( " double "+x for x in symbolslist )

    c_code = "double "+str(symname)+"("+varstring+" )\n"
    c_code +=  "{\n"
    for s in symbols:
        #print s
        c_code +=  "  double " +ccode(s[0]) + " = " + ccode(s[1]) + ";\n"
    c_code +=  "  r = " + ccode(simple[0])+";\n"
    c_code +=  "  return r;\n"
    c_code += "}\n"
    return c_code

还有适用于python3.5及以上版本的：

def sympyToC( symname, symfunc ):
    tmpsyms = numbered_symbols("tmp")
    symbols, simple = cse(symfunc, symbols=tmpsyms)
    symbolslist = sorted(map(lambda x:str(x), list(symfunc.atoms(Symbol))))
    varstring=",".join( " double "+x for x in symbolslist )
    c_code = "double "+str(symname)+"("+varstring+" )\n"
    c_code +=  "{\n"
    for s in symbols:
        #print s
        c_code +=  "  double " +ccode(s[0]) + " = " + ccode(s[1]) + ";\n"
    c_code +=  "  r = " + ccode(simple[0])+";\n"
    c_code +=  "  return r;\n"
    c_code += "}\n"
    return c_code

Answer 2

gcc可能会自动优化这个，但如果你想自己动手做，可以看看cse。你可以访问这个链接了解更多信息：http://docs.sympy.org/latest/modules/simplify/simplify.html#module-sympy.simplify.cse_main