c++ - 在 C++ 中实现长方程时，如何通过高级方法提高性能

我正在开发一些工程模拟。这涉及实现一些长方程，例如这个方程来计算橡胶类 Material 中的应力:

T = (
    mu * (
            pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
            * (
                pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
                - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
            ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l1
            - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
            - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
        ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l2 * l3
) * N1 / l2 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
        + pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l2
        - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
    ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l1 * l3
) * N2 / l1 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        + pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l3
    ) / a
+ K * (l1 * l2 * l3 - 0.1e1) * l1 * l2
) * N3 / l1 / l2;

我使用 Maple 生成 C++ 代码以避免错误(并节省繁琐的代数时间)。由于此代码执行数千次(如果不是数百万次)，性能是一个问题。不幸的是，到目前为止，数学只是简化了；长方程是不可避免的。

我可以采取什么方法来优化此实现？ 我正在寻找在实现这些方程时应该应用的高级策略，不一定是上面显示的示例的特定优化。

我正在使用 g++ 和 --enable-optimize=-O3 进行编译.

更新:

我知道有很多重复的表达式，我假设编译器会处理这些；到目前为止，我的测试表明确实如此。
l1, l2, l3, mu, a, K都是正实数(不是零)。

我已更换 l1*l2*l3具有等效变量:J .这确实有助于提高性能。

更换 pow(x, 0.1e1/0.3e1)与 cbrt(x)是个好建议。

这将在 CPU 上运行，在不久的将来这可能会在 GPU 上运行得更好，但目前该选项不可用。

最佳答案

编辑摘要

我最初的回答只是指出代码包含大量重复计算，并且许多幂涉及 1/3 的因数。例如，pow(x, 0.1e1/0.3e1)与 cbrt(x) 相同.

我的第二次编辑是错误的，而我的第三次编辑是对这种错误进行推断的。这就是让人们害怕改变以字母“M”开头的符号数学程序的类似预言的结果的原因。我已经删除(即删除)这些编辑并将它们推送到此答案当前修订版的底部。但是，我没有删除它们。我是人类。我们很容易犯错误。

我的第四次编辑开发了一个非常紧凑的表达式，它正确地表示了问题 中复杂的表达式。如果 参数l1 , l2 , 和 l3是正实数，如果 a是一个非零实数。 (我们还没有收到 OP 关于这些系数的具体性质的消息。鉴于问题的性质，这些都是合理的假设。)

此编辑试图回答如何简化这些表达式的一般问题。

第一件事

I use Maple to generate the C++ code to avoid mistakes.

Maple 和 Mathematica 有时会错过显而易见的事情。更重要的是，Maple 和 Mathematica 的用户有时会犯错误。用“经常”，甚至“几乎总是”代替“有时可能更接近标记”。

您可以通过告诉 Maple 相关参数来帮助它简化该表达式。在手头的例子中，我怀疑 l1 , l2 , 和 l3是正实数且 a是一个非零实数。如果是这种情况，请告诉它。这些符号数学程序通常假设手头的数量是复杂的。限制域允许程序做出在复数中无效的假设。

如何简化符号数学程序中的那些大麻烦(此编辑)

符号数学程序通常能够提供有关各种参数的信息。使用这种能力，特别是如果您的问题涉及除法或幂运算。在手头的示例中，您可以通过告诉 Maple l1 来帮助 Maple 简化该表达式。 , l2 , 和 l3是正实数且 a是一个非零实数。如果是这种情况，请告诉它。这些符号数学程序通常假设手头的数量是复杂的。限制域让程序做出诸如 axbx=(ab)x 之类的假设。这仅当 a和 b是正实数，如果 x是真的。它在复数中无效。

最终，那些符号数学程序遵循算法。一起帮它。在生成代码之前尝试扩展、收集和简化。在这种情况下，您可以收集涉及因子 mu 的那些项。以及涉及 K 因子的那些.将表达式简化为“最简单的形式”仍然是一门艺术。

当你得到一堆丑陋的生成代码时，不要把它当作你不能触及的真理。尝试自己简化。看看符号数学程序在生成代码之前有什么。看看我是如何将你的表达简化为更简单、更快的，以及 Walter's answer让我更进一步。没有神奇的配方。如果有一个神奇的配方，枫就会应用它并给出沃尔特给出的答案。

关于具体问题

你在那个计算中做了很多加法和减法。如果您的条款几乎相互取消，您可能会遇到很大的麻烦。如果您有一个术语支配其他术语，那么您将浪费大量 CPU。

接下来，您通过执行重复计算浪费了大量 CPU。除非您已启用 -ffast-math ，这让编译器打破了 IEEE 浮点数的一些规则，编译器不会(实际上，不能)为您简化该表达式。相反，它会完全按照您的吩咐去做。至少，您应该计算 l1 * l2 * l3在计算那个烂摊子之前。

最后，您经常拨打 pow 的电话。，这是非常缓慢的。请注意，其中一些调用的形式为 (l1*l2*l3)(1/3)。其中许多调用 pow只需一次调用 std::cbrt 即可执行:

l123 = l1 * l2 * l3;
l123_pow_1_3 = std::cbrt(l123);
l123_pow_4_3 = l123 * l123_pow_1_3;

有了这个，

X * pow(l1 * l2 * l3, 0.1e1 / 0.3e1)变成 X * l123_pow_1_3 .

X * pow(l1 * l2 * l3, -0.1e1 / 0.3e1)变成 X / l123_pow_1_3 .

X * pow(l1 * l2 * l3, 0.4e1 / 0.3e1)变成 X * l123_pow_4_3 .

X * pow(l1 * l2 * l3, -0.4e1 / 0.3e1)变成 X / l123_pow_4_3 .

Maple 确实错过了显而易见的事情。
例如，有一种更简单的写法

(pow(l1 * l2 * l3, -0.1e1 / 0.3e1) - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)

假设 l1 , l2 , 和 l3是实数而不是复数，并且要提取实数立方根(而不是主复数根)，以上简化为

2.0/(3.0 * pow(l1 * l2 * l3, 1.0/3.0))

或

2.0/(3.0 * l123_pow_1_3)

使用 cbrt_l123而不是 l123_pow_1_3 ，问题中令人讨厌的表达简化为

l123 = l1 * l2 * l3; 
cbrt_l123 = cbrt(l123);
T = 
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);

始终仔细检查，但也始终简化。

以下是我达到上述目标的一些步骤:

// Step 0: Trim all whitespace.
T=(mu*(pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1+pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l2-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1+pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l3)/a+K*(l1*l2*l3-0.1e1)*l1*l2)*N3/l1/l2;

// Step 1:
//   l1*l2*l3 -> l123
//   0.1e1 -> 1.0
//   0.4e1 -> 4.0
//   0.3e1 -> 3
l123 = l1 * l2 * l3;
T=(mu*(pow(l1*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l1-pow(l2*pow(l123,-1.0/3),a)*a/l1/3-pow(l3*pow(l123,-1.0/3),a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l2/3+pow(l2*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l2-pow(l3*pow(l123,-1.0/3),a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l3/3-pow(l2*pow(l123,-1.0/3),a)*a/l3/3+pow(l3*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 2:
//   pow(l123,1.0/3) -> cbrt_l123
//   l123*pow(l123,-4.0/3) -> pow(l123,-1.0/3)
//   (pow(l123,-1.0/3)-pow(l123,-1.0/3)/3) -> 2.0/(3.0*cbrt_l123)
//   *pow(l123,-1.0/3) -> /cbrt_l123
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T=(mu*(pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1-pow(l2/cbrt_l123,a)*a/l1/3-pow(l3/cbrt_l123,a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l2/3+pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2-pow(l3/cbrt_l123,a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l3/3-pow(l2/cbrt_l123,a)*a/l3/3+pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 3:
//   Whitespace is nice.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)*a/l1/3
       -pow(l3/cbrt_l123,a)*a/l1/3)/a
   +K*(l123-1.0)*l2*l3)*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l2/3
       +pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)*a/l2/3)/a
   +K*(l123-1.0)*l1*l3)*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l3/3
       -pow(l2/cbrt_l123,a)*a/l3/3
       +pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a
   +K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 4:
//   Eliminate the 'a' in (term1*a + term2*a + term3*a)/a
//   Expand (mu_term + K_term)*something to mu_term*something + K_term*something
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +K*(l123-1.0)*l2*l3*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +K*(l123-1.0)*l1*l3*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l3))*N3/l1/l2
 +K*(l123-1.0)*l1*l2*N3/l1/l2;

// Step 5:
//   Rearrange
//   Reduce l2*l3*N1/l2/l3 to N1 (and similar)
//   Reduce 2.0/(3.0*cbrt_l123)*cbrt_l123/l1 to 2.0/3.0/l1 (and similar)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/3.0/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/3.0/l3))*N3/l1/l2
 +K*(l123-1.0)*N1
 +K*(l123-1.0)*N2
 +K*(l123-1.0)*N3;

// Step 6:
//   Factor out mu and K*(l123-1.0)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*(  ( pow(l1/cbrt_l123,a)*2.0/3.0/l1
         -pow(l2/cbrt_l123,a)/l1/3
         -pow(l3/cbrt_l123,a)/l1/3)*N1/l2/l3
      + (-pow(l1/cbrt_l123,a)/l2/3
         +pow(l2/cbrt_l123,a)*2.0/3.0/l2
         -pow(l3/cbrt_l123,a)/l2/3)*N2/l1/l3
      + (-pow(l1/cbrt_l123,a)/l3/3
         -pow(l2/cbrt_l123,a)/l3/3
         +pow(l3/cbrt_l123,a)*2.0/3.0/l3)*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 7:
//   Expand
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1*N1/l2/l3
      -pow(l2/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l3/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l1/cbrt_l123,a)/l2/3*N2/l1/l3
      +pow(l2/cbrt_l123,a)*2.0/3.0/l2*N2/l1/l3
      -pow(l3/cbrt_l123,a)/l2/3*N2/l1/l3
      -pow(l1/cbrt_l123,a)/l3/3*N3/l1/l2
      -pow(l2/cbrt_l123,a)/l3/3*N3/l1/l2
      +pow(l3/cbrt_l123,a)*2.0/3.0/l3*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 8:
//   Simplify.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);

错误答案，故意保留以示谦虚

请注意，这是受到打击。这是错误的。

更新

Maple 确实错过了显而易见的事情。例如，有一种更简单的写法

(pow(l1 * l2 * l3, -0.1e1 / 0.3e1) - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)

假设 l1 , l2 , 和 l3是实数而不是复数，并且要提取实数立方根(而不是主复数根)，则上述减少为零。这种零的计算重复多次。

第二次更新

如果我做对了数学(有没有保证我做对了数学)，问题中令人讨厌的表达会简化为

l123 = l1 * l2 * l3; 
cbrt_l123_inv = 1.0 / cbrt(l123);
nasty_expression =
    K * (l123 - 1.0) * (N1 + N2 + N3) 
    - (  pow(l1 * cbrt_l123_inv, a) * (N2 + N3) 
       + pow(l2 * cbrt_l123_inv, a) * (N1 + N3) 
       + pow(l3 * cbrt_l123_inv, a) * (N1 + N2)) * mu / (3.0*l123);

以上假设l1 , l2 , 和 l3是正实数。

关于c++ - 在 C++ 中实现长方程时，如何通过高级方法提高性能，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/32907258/

c++ - 在 C++ 中实现长方程时，如何通过高级方法提高性能

上一篇：c++ - 什么是 SDL 渲染器？

下一篇：c++ - 如何为 C++ 中的函数名分配别名？