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

Question

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

91 c++ floating-point optimization performance g++

我正在开发一些工程模拟.这包括实现一些长方程,例如这个方程,以计算橡胶类材料中的应力:

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;

Run Code Online (Sandbox Code Playgroud)

我使用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上运行得更好,但目前该选项不可用.

Answer 1

Dav*_*men 87

编辑摘要

我的原始答案仅仅指出代码包含大量复制计算,并且许多权力涉及1/3的因素.例如,pow(x, 0.1e1/0.3e1)与...相同cbrt(x).
我的第二次编辑是错误的,而我的第三次编辑推断了这种错误.这使得人们害怕改变以字母'M'开头的符号数学程序的类似oracle的结果.我已经敲除了(即,敲击)那些编辑并将它们推到了当前修订版本的底部.但是,我没有删除它们.我是人类.我们很容易犯错误.
我的第四个编辑开发了一个非常紧凑的表达式正确代表的问题令人费解的表情IF参数l1,l2和l3是正实数,如果a是非零实数.(我们还没有听到OP关于这些系数的具体性质.鉴于问题的性质,这些都是合理的假设.)
此编辑试图回答如何简化这些表达式的一般问题.

首先要做的事情

我使用Maple生成C++代码以避免错误.

Maple和Mathematica有时会错过显而易见的事实.更重要的是,Maple和Mathematica的用户有时会犯错误.代替"经常",或者甚至"几乎总是"代替"有时可能更接近标记".

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

如何从符号数学程序中简化那些大混乱(本编辑)

符号数学程序通常提供提供有关各种参数的信息的能力.使用该能力,特别是如果您的问题涉及分裂或取幂.在手边的例子,你可以帮助枫告诉它该简化的表达l1,l2和l3是正实数,并且a是非零实数.如果是这种情况,请告诉它.这些符号数学程序通常假设手头的数量很复杂.限制域允许程序进行假设,例如^x b ^x =(ab)^x.这只是当a且b是正实数且如果x是真实的.它在复数中无效.

最终,这些符号数学程序遵循算法.帮助它.在生成代码之前尝试使用扩展,收集和简化.在这种情况下,你可以收集那些涉及的因素方面mu以及涉及的因素K.将表达式简化为"最简单的形式"仍然是一门艺术.

当你得到一堆丑陋的生成代码时,不要接受它作为你不能触及的事实.尽量简化它.看一下符号数学程序在生成代码之前的作用.看看我如何将你的表达简化为更简单,更快速的东西,以及沃尔特的答案如何进一步采取了我的几个步骤.没有神奇的食谱.如果有一个神奇的食谱,Maple会应用它并给出Walter给出的答案.

关于具体问题

你在那个计算中做了很多加法和减法.如果您的条款几乎相互抵消,您可能会陷入深深的麻烦.如果你有一个术语比其他术语占主导地位,你就浪费了大量的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;

Run Code Online (Sandbox Code Playgroud)

有了这个,

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.

枫树确实错过了显而易见的事.
例如,有一种更容易编写的方法

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

Run Code Online (Sandbox Code Playgroud)

假设l1,l2和l3是实而非复数,而真正立方根(而不是原则复杂根)将被提取的,上述减少到

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

Run Code Online (Sandbox Code Playgroud)

要么

2.0/(3.0 * l123_pow_1_3)

Run Code Online (Sandbox Code Playgroud)

使用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);

Run Code Online (Sandbox Code Playgroud)

总是仔细检查,但总是简化.

以下是我实现上述步骤的一些步骤:

// 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);

Run Code Online (Sandbox Code Playgroud)

错误的答案,故意保持谦虚

请注意,这很糟糕.这是不对的.

更新

枫树确实错过了显而易见的事.例如,有一种更容易编写的方法

(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);
Run Code Online (Sandbox Code Playgroud)

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

@Deduplicator - 没有浮点数.除非有人启用不安全的数学优化(例如,通过使用gcc或clang指定`-ffast-math`),编译器不能依赖于`pow(x,-1.0/3.0)`等于`x*pow(x, - 4.0/3.0)`.后者可能会下流,而第一个可能不会.为了符合浮点标准,编译器不得将该计算优化为零. (3认同)
好吧,CSE消除*应该*工作,独立于轻松的语义(并且在评论中表示OP).当然,如果它很重要(测量),那应该进行检查(生成装配).关于主导术语,错过公式简化,更好的专业功能和取消的危险,你的观点非常好. (2认同)

Answer 2

mar*_*sky 31

首先要注意的是,这pow是非常昂贵的,所以你应该尽可能地摆脱这个.扫描表达式我看到很多重复的pow(l1 * l2 * l3, -0.1e1 / 0.3e1)和pow(l1 * l2 * l3, -0.4e1 / 0.3e1).所以我希望从预先计算那些:

 const double c1 = pow(l1 * l2 * l3, -0.1e1 / 0.3e1);
const double c2 = boost::math::pow<4>(c1);

Run Code Online (Sandbox Code Playgroud)

这里我使用升压战俘功能.

此外,你还有更多pow的指数a.如果a是Integer且在编译器时已知,则还可以替换它们boost::math::pow<a>(...)以获得进一步的性能.我还建议替换术语a / l1 / 0.3e1,a / (l1 * 0.3e1)因为乘法比分割更快.

最后,如果使用g ++,您可以使用-ffast-math允许优化器在转换方程式时更积极的标志.阅读这个标志实际上做了什么,因为它有副作用.

@usr:我认为这就是重点.根据标准,`pow`不是*纯函数,因为它应该在某些情况下设置`errno`.设置`-fno-math-errno`等标志会导致*不*设置`errno`(因此违反标准),但它是一个纯函数,可以这样优化. (5认同)
在我们的代码中,使用`-ffast-math`导致代码变得不稳定或者给出错误的答案.我们在使用英特尔编译器时遇到类似的问题,并且必须使用`-fp-model precise`选项,否则代码会爆炸或给出错误的答案.因此`-ffast-math`可以加快速度,但除了链接问题中列出的副作用外,我还建议您对该选项进行非常谨慎的处理. (4认同)
@tpg2114：[根据我的测试，你只需要`-fno-math-errno`](http://stackoverflow.com/questions/2940367/what-is-more-efficient-using-pow-to-square- or-just-multiply-it-with-itself/2940800#comment53658268_2940800) 使 g++ 能够在循环中提升对 `pow` 的相同调用。对于大多数代码，这是 -ffast-math 中最不“危险”的部分。 (2认同)

Answer 3

cdo*_*nat 20

哇,多么神奇的表达.用Maple创建表达式实际上是次优选择.结果简直难以理解.

选择说话变量名称(不是l1,l2,l3,但是例如高度,宽度,深度,如果这是他们的意思).然后,您更容易理解自己的代码.
计算多次使用的子项,预先将结果存储在带有说话名称的变量中.
你提到,表达式被评估了很多次.我想,在最内层循环中只有少数参数不同.在该循环之前计算所有不变子项.重复第二个内循环,依此类推,直到所有不变量都在循环之外.

从理论上讲,编译器应该能够为您完成所有这些工作,但有时它不能 - 例如,当循环嵌套在不同的编译单元中扩展多个函数时.无论如何,这将为您提供更好的可读性,可理解性和可维护性的代码.

"编译器应该这样做,但有时它不会",这是关键.当然,除了可读性之外. (8认同)
"结果根本不可读" - 为什么这是一个问题？你不会关心词法分析器或解析器生成器的高级语言输出是"不可读的"(人类).这里重要的是代码生成器*(Maple)的*输入是可读和可检查的.如果您想确信它没有错误,那么*不要做的事情就是手动编辑生成的代码. (8认同)
重新*选择说话变量名* - 很多时候,当你在做数学时,那个好的规则不适用.在查看应该在科学期刊中实现算法的代码时,我更愿意看到代码中的符号正是日记中使用的符号.通常,这意味着极短的名称,可能带有下标. (4认同)
如果编译器不需要做某事,那么假设几乎总是错误的. (3认同)
@DavidHammen:嗯,在这种情况下,单字母的*是*"说出名字".例如,当在二维笛卡尔坐标系中进行几何时,"x"和"y"不是*无意义的单字母变量,它们是整个*单词*,具有精确的定义和良好且广泛理解的含义. (3认同)

Answer 4

Wal*_*ter 17

大卫·哈门的答案很好,但仍远未达到最佳状态.让我们继续他的最后一个表达(在撰写本文时)

auto l123 = l1 * l2 * l3;
auto 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);

Run Code Online (Sandbox Code Playgroud)

可以进一步优化.特别是,如果利用一些数学身份,我们可以避免调用cbrt()和调用其中一个pow().让我们一步一步地再做一遍.

// step 1 eliminate cbrt() by taking the exponent into pow()
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a; // avoid division
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)*pow(l1*l1/(l2*l3),athird)
                   + (N2+N2-N3-N1)*pow(l2*l2/(l1*l3),athird)
                   + (N3+N3-N1-N2)*pow(l3*l3/(l1*l2),athird))
  + K*(l123-1.0)*(N1+N2+N3);

Run Code Online (Sandbox Code Playgroud)

请注意,我还优化2.0*N1了N1+N1等.接下来,我们只能进行两次调用pow().

// step 2  eliminate one call to pow
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a;
auto pow_l1l2_athird = pow(l1/l2,athird);
auto pow_l1l3_athird = pow(l1/l3,athird);
auto pow_l2l3_athird = pow_l1l3_athird/pow_l1l2_athird;
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)* pow_l1l2_athird*pow_l1l3_athird
                   + (N2+N2-N3-N1)* pow_l2l3_athird/pow_l1l2_athird
                   + (N3+N3-N1-N2)/(pow_l1l3_athird*pow_l2l3_athird))
  + K*(l123-1.0)*(N1+N2+N3);

Run Code Online (Sandbox Code Playgroud)

由于呼叫pow()是迄今为止最昂贵的操作,因此尽可能减少它们是值得的(下一个昂贵的操作是cbrt()我们消除的呼叫).

如果任何机会a是整数,则调用pow可以优化为调用cbrt(加上整数幂),或者如果athird是半整数,我们可以使用sqrt(加上整数幂).此外,如果任何机会l1==l2或l1==l3或l2==l3一个或两个调用pow可以消除.因此,如果现实存在这样的机会,将这些视为特殊情况是值得的.

Answer 5

Vla*_*ein 12

有多少是"很多"？
多久时间？
难道所有的参数这个公式重新计算之间的变化呢？或者你可以缓存一些预先计算的值吗？

我试图手动简化该公式,想知道它是否可以保存任何东西？

C1 = -0.1e1 / 0.3e1;
C2 =  0.1e1 / 0.3e1;
C3 = -0.4e1 / 0.3e1;

X0 = l1 * l2 * l3;
X1 = pow(X0, C1);
X2 = pow(X0, C2);
X3 = pow(X0, C3);
X4 = pow(l1 * X1, a);
X5 = pow(l2 * X1, a);
X6 = pow(l3 * X1, a);
X7 = a / 0.3e1;
X8 = X3 / 0.3e1;
X9 = mu / a;
XA = X0 - 0.1e1;
XB = K * XA;
XC = X1 - X0 * X8;
XD = a * XC * X2;

XE = X4 * X7;
XF = X5 * X7;
XG = X6 * X7;

T = (X9 * ( X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
  + (X9 * (-XE + X5 * XD - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
  + (X9 * (-XE - XF + X6 * XD) / l3 + XB * l1 * l2) * N3 / l1 / l2;

Run Code Online (Sandbox Code Playgroud)

[增加]我已经在最后三行公式上做了更多的工作,并把它归结为这个美丽:

T = X9 / X0 * (
      (X4 * XD - XF - XG) * N1 + 
      (X5 * XD - XE - XG) * N2 + 
      (X5 * XD - XE - XF) * N3)
  + XB * (N1 + N2 + N3)

Run Code Online (Sandbox Code Playgroud)

让我一步一步地展示我的工作:

T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / l1 / l2;


T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / (l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / (l1 * l3) 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / (l1 * l2);

T = (X9 * (X4 * XD - XF - XG) + XB * l1 * l2 * l3) * N1 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) + XB * l1 * l2 * l3) * N2 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XF) + XB * l1 * l2 * l3) * N3 / (l1 * l2 * l3);

T = (X9 * (X4 * XD - XF - XG) + XB * X0) * N1 / X0 
  + (X9 * (X5 * XD - XE - XG) + XB * X0) * N2 / X0 
  + (X9 * (X5 * XD - XE - XF) + XB * X0) * N3 / X0;

T = X9 * (X4 * XD - XF - XG) * N1 / X0 + XB * N1 
  + X9 * (X5 * XD - XE - XG) * N2 / X0 + XB * N2
  + X9 * (X5 * XD - XE - XF) * N3 / X0 + XB * N3;


T = X9 * (X4 * XD - XF - XG) * N1 / X0 
  + X9 * (X5 * XD - XE - XG) * N2 / X0
  + X9 * (X5 * XD - XE - XF) * N3 / X0
  + XB * (N1 + N2 + N3)

Run Code Online (Sandbox Code Playgroud)

那明显,是吧？:) FORTRAN,IIRC,专为高效的公式计算而设计("FOR"用于公式). (2认同)

Answer 6

neo*_*cpp 7

这可能有点简洁,但我实际上通过使用Horner Form找到了多项式(能量函数的插值)的良好加速,其基本上重写ax^3 + bx^2 + cx + d为d + x(c + x(b + x(a))).这将避免大量的重复调用,pow()并阻止你做一些愚蠢的事情,比如单独调用pow(x,6)而pow(x,7)不是仅仅做x*pow(x,6).

这不是直接适用于您当前的问题,但如果您有具有整数幂的高阶多项式,它可以提供帮助.您可能必须注意数值稳定性和溢出问题,因为操作顺序对此非常重要(尽管通常我认为Horner Form对此有帮助,因为x^20并且x通常相隔很多个数量级).

另外作为一个实用的提示,如果你还没有这样做,请先尝试简化maple中的表达式.您可以让它为您执行大多数常见的子表达式消除.我不知道它对该程序中的代码生成器有多大影响,但我知道Mathematica在生成代码之前执行FullSimplify可能会产生巨大的差异.

归档时间：	10 年，2 月前
查看次数：	4436 次
最近记录：	10 年，2 月前