咨询Sympy生成任意尺寸矩阵C代码时Sum函数的正确转换方法

问题背景 #

尝试用Sympy生成包含任意大小矩阵逻辑的C代码时，发现Sum函数未被正确转换为可执行的C代码，生成的代码直接保留了Sum符号，导致代码无法编译运行。

复现代码（Python） #

from sympy import symbols, Sum
from sympy.utilities.codegen import codegen

n = symbols('n', integer=True)
i = symbols('i')
my_sum = Sum(i, (i,0,n))
pprint(codegen(("f", my_sum), "C", "test", header=False, empty=False)[0][1])

生成的无效C代码 #

#include "test.h"
#include <math.h>
double f(int n) {
   double f_result;
   f_result = Sum(i, (i, 0, n));
   return f_result;
}

预期的正确逻辑示例（C++） #

用户提供的符合预期的C++实现（已修正语法错误）：

#include <math.h>
#include <functional>
#include <ranges>
#include <algorithm>

double f(int n) {
   double f_result;

   auto sympy_sum_expr = [](int i) {return static_cast<double>(i);};
   auto sympy_sum = [](std::function<double(int)> expr, int from, int to) -> double 
   {
      double sum = 0.0;
      std::for_each(std::ranges::iota_view(from, to + 1).begin(), 
                    std::ranges::iota_view(from, to + 1).end(),
                    [&](int i) { sum += expr(i); });
      return sum;
   };

   f_result = sympy_sum(sympy_sum_expr, 0, n);
   return f_result;
}

解决方法 #

Sympy的codegen模块默认不会将Sum展开为循环代码，它更偏向于直接映射数学表达式而非生成迭代逻辑。针对不同场景，有两种可行方案：

方案1：化简为闭合数学公式（优先选择） #

对于存在闭合解的简单求和，先通过Sympy将Sum化简为代数表达式，再生成代码：

from sympy import symbols, Sum, simplify
from sympy.utilities.codegen import codegen

n = symbols('n', integer=True)
i = symbols('i')
my_sum = Sum(i, (i,0,n))
# 化简求和公式为 n*(n+1)/2
simplified_sum = simplify(my_sum)
pprint(codegen(("f", simplified_sum), "C", "test", header=False, empty=False)[0][1])

生成的有效C代码：

#include "test.h"
#include <math.h>
double f(int n) {
   double f_result;
   f_result = n*(n + 1)/2.0;
   return f_result;
}

这种方式效率最高，适合有明确代数解的求和场景。

方案2：手动生成循环代码（适合无闭合解的复杂求和） #

如果求和无法化简为闭合公式，可通过解析Sum对象的参数，手动拼接循环形式的C代码：

from sympy import symbols, Sum

n = symbols('n', integer=True)
i = symbols('i')
my_sum = Sum(i, (i,0,n))

# 提取求和变量、上下界和表达式
sum_var = my_sum.variables[0]
lower_bound = my_sum.limits[0][1]
upper_bound = my_sum.limits[0][2]
expr = my_sum.function

# 生成循环代码
c_code = f"""#include <math.h>
double f(int n) {{
   double f_result = 0.0;
   for (int {sum_var} = {lower_bound}; {sum_var} <= {upper_bound}; {sum_var}++) {{
      f_result += {expr};
   }}
   return f_result;
}}"""

print(c_code)

输出的有效C代码：

#include <math.h>
double f(int n) {
   double f_result = 0.0;
   for (int i = 0; i <= n; i++) {
      f_result += i;
   }
   return f_result;
}

这种方式可以灵活处理复杂求和表达式，根据Sum的参数动态生成迭代逻辑。

内容的提问来源于stack exchange，提问作者Ingo Thomas