在C ++中实现长方程时，如何通过高级方法提高性能？
我正在开发一些工程模拟。 这涉及到实施一些长方程，如这个方程来计算橡胶材料的应力：
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 ++代码以避免错误（并且用繁琐的代数节省时间）。 由于此代码执行数千次（如果不是数百万次），所以性能是一个问题。 不幸的是，math只是简化到目前为止; 长方程是不可避免的。
我可以采取什么方法来优化这个实现？ 我正在寻找在实施这样的方程时应该应用的高级策略，而不一定是针对上述示例的具体优化。
我正在编译使用g ++与--enable-optimize=-O3 。
更新：
我知道有很多重复的expression式，我正在使用编译器会处理这些的假设; 我迄今为止的testing表明它的确如此。
l1, l2, l3, mu, a, K都是正实数（不是零）。
我已经用等效variablesreplace了l1*l2*l3 ： J 。 这确实有助于提高性能。
用cbrt(x)代替pow(x, 0.1e1/0.3e1) cbrt(x)是一个很好的build议。
这将在CPU上运行，在不久的将来，这可能会在GPU上运行得更好，但现在该选项不可用。
编辑摘要
- 我原来的回答只是指出，代码中包含了大量的重复计算，而且其中的许多权力都涉及1/3的因素。 例如，
pow(x, 0.1e1/0.3e1)与cbrt(x)相同。 - 我的第二个编辑是错误的，而我的第三个编辑就是这个错误。 这就是人们害怕从以字母“M”开始的符号math程序改变甲骨文式结果的原因。 我击败了（即
罢工）这些编辑，并将其推到当前版本的答案的底部。 但是，我没有删除它们。 我是人类。 我们很容易犯错。 - 我的第四个编辑开发了一个非常紧凑的expression式，正确地表示了问题中的复杂expression式， 如果参数
l1，l2和l3是正实数且a是非零实数。 （关于这些系数的具体性质，我们还没有听说OP，鉴于问题的性质，这些都是合理的假设。） - 此编辑尝试回答如何简化这些expression式的一般问题。
首先是事情
我使用Maple生成C ++代码以避免错误。
枫和Mathematica有时会错过显而易见的事物。 更重要的是，枫和Mathematica的用户有时会犯错误。 代替“时常”或甚至“几乎总是”代替“有时可能更接近标记”。
你可以帮助Maple通过告诉它关于参数的问题来简化expression式。 在这个例子中，我怀疑l1 ， l2和l3是正实数， a是一个非零的实数。 如果是这样，告诉它。 那些象征性的math课程通常假定手头的数量是复杂的。 限制域允许程序进行假设，在复数中无效。
如何从符号math程序中简化这些大混乱（这个编辑）
符号math程序通常提供提供有关各种参数信息的能力。 使用这种能力，特别是如果你的问题涉及分裂或求幂。 在这个例子中，你可以帮助Maple简单地通过告诉它， l1 ， l2和l3是正实数，并且a是一个非零的实数。 如果是这样，告诉它。 那些象征性的math课程通常假定手头的数量是复杂的。 限制域允许程序做出假设，如x b x =（ab） x 。 这只有在a和b是正实数且x是实数的情况下。 这在复数中是无效的。
最终，这些符号math程序遵循algorithm。 帮助它。 在生成代码之前，尝试使用扩展，收集和简化操作。 在这种情况下，你可能已经收集了涉及一个因子mu的因素和涉及因子K的因子。 减lessexpression“最简单的forms”仍然是一个艺术。
当你得到一个丑陋的生成代码，不要接受它是一个你不能碰的事实。 尝试自己简化它。 看看符号math程序在生成代码之前有什么。 看看我如何将你的表情减less到更简单，更快速的状态，以及沃尔特的回答如何进一步把我的几个步骤。 没有魔法配方。 如果有一个神奇的配方，枫树会应用它，并给出了沃尔特的答案。
关于具体的问题
在这个计算中你正在做很多的加法和减法。 如果您有几乎相互抵消的条款，您可能会遇到很大的麻烦。 如果你有一个主宰其他的术语，你正在浪费很多的CPU。
接下来，通过执行重复的计算，你正在浪费大量的CPU。 除非你已经启用了-ffast-math ，这让编译器破坏了IEEE浮点的一些规则，编译器不会（实际上不能）为你简化这个expression式。 它将完全按照你所说的去做。 至less，你应该计算l1 * l2 * l3然后计算这个混乱。
最后，你打了很多电话，这是非常缓慢的。 请注意，其中几个电话的forms（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。
枫叶确实错过了明显的。
例如，写一个更简单的方法
(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 ，问题中令人讨厌的expression式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);
错误的答案，故意保持谦卑
请注意，这是灾难。 这是不对的。
更新
枫叶确实错过了明显的。 例如，写一个更简单的方法
（pow（l1 * l2 * l3，-0.1e1 / 0.3e1）-l1 * l2 * l3 * pow（l1 * l2 * l3，-0.4e1 / 0.3e1）/ 0.3e1）
假设l1 ， l2和l3是实数而不是复数，并且要提取实数立方根（而不是复数根的原理），则上述数值减less到零。 零的计算重复多次。
第二次更新
如果我已经完成了math计算（ 不能保证我已经完成了math计算），那么问题中令人讨厌的expression就会降低
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是正实数。
首先要注意的是， pow是非常昂贵的，所以你应该尽可能地去掉这个。 通过expression式扫描，我看到很多重复的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);
我正在使用boost pow函数。
而且，你有更多的指数a 。 如果a是Integer并且在编译器时间已知，那么也可以用boost::math::pow<a>(...)来replace它们以获得进一步的性能。 我还build议用a / (l1 * 0.3e1)replacea / l1 / 0.3e1这样a / l1 / 0.3e1 ，因为乘法运算速度快于除法运算。
最后，如果使用g ++，则可以使用-ffast-math标志，以使优化程序在转换方程时更加积极。 阅读这个标志实际上做了什么 ，因为它有副作用。
哇，这真是一种expression。 用Maple创buildexpression式实际上是一个次优的select。 结果是根本不可读的。
- select了说variables名称（不是l1，l2，l3，但是例如高度，宽度，深度，如果这就是他们的意思）。 那么你就更容易理解你自己的代码。
- 计算多次使用的子项，并将结果存储在带有说出名称的variables中。
- 你提到，expression式被评估了很多次。 我想，在最内层循环中只有很less的参数有所不同。 在该循环之前计算所有不变的子项。 重复第二个内部循环，依此类推，直到所有不variables都在循环之外。
理论上编译器应该能够为你做所有这些工作，但是有时却不能 – 例如当循环嵌套遍布不同编译单元的多个函数时。 无论如何，这将给你更好的可读性，可理解性和可维护性的代码。
大卫·哈姆曼的答案是好的，但仍然远非最佳。 让我们继续他的最后一个expression（写这篇文章的时候）
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);
这可以进一步优化。 特别是，如果利用一些math标识，我们可以避免调用cbrt()和调用pow()其中一个cbrt() 。 我们一步一步来做吧。
// 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);
请注意，我还优化了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);
由于对pow()的调用是迄今为止成本最高的操作，所以尽可能地减less它（下一个代价高昂的操作是调用cbrt() ，这是我们消除的）。
如果是整数，那么对pow的调用可以优化为调用cbrt （加上整数幂），或者如果athird是半整数，我们可以使用sqrt （加上整数幂）。 此外，如果任何机会l1==l2或l1==l3或l2==l3一个或两个调用pow可以被消除。 所以，如果这样的机会真实存在的话，值得把这些视为特例。
- “多less”有多less？
- 多久时间？
- 重新计算这个公式之间的所有参数是否改变？ 或者你可以caching一些预先计算的值？
-
我试图手动简化该公式，想知道它是否可以保存任何东西？
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;
[ADDED]我已经在最后三行公式上做了更多的工作，并且把它归结为这个美：
T = X9 / X0 * ( (X4 * XD - XF - XG) * N1 + (X5 * XD - XE - XG) * N2 + (X5 * XD - XE - XF) * N3) + XB * (N1 + N2 + N3)
让我一步一步展示我的工作：
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)
这可能有点简单，但是我实际上已经发现使用Hornerforms对多项式（能量函数的插值）进行了很好的加速，其基本上将ax^3 + bx^2 + cx + d重写为d + x(c + x(b + x(a))) 。 这将避免大量重复调用pow()并阻止你像做单独调用pow(x,6)和pow(x,7)那样pow(x,7)而不是仅仅做x*pow(x,6) 。
这不是直接适用于你当前的问题，但如果你有整数次幂的高阶多项式，它可以帮助。 你可能不得不小心数值稳定性和溢出问题，因为操作的顺序对于这个很重要（尽pipe一般来说我实际上认为Hornerforms对此有所帮助，因为x^20和x通常相隔很多个数量级）。
另外作为一个实用的提示，如果你还没有这样做，试着简化枫的expression式首先。 你可能可以让它去做大部分常见的子expression式消除。 我不知道在程序中会影响代码生成器的多less，但是我知道在生成代码之前做一个FullSimplify的Mathematica会导致巨大的差异。
看起来你有很多重复的操作正在进行。
pow(l1 * l2 * l3, -0.1e1 / 0.3e1) pow(l1 * l2 * l3, -0.4e1 / 0.3e1)
你可以预先计算这些，所以你没有反复调用pow函数，这可能是昂贵的。
你也可以预先称呼
l1 * l2 * l3
因为你反复使用这个术语。
如果你有一个Nvidia CUDA显卡，你可以考虑将计算任务卸载到显卡上 – 这本身就更适合计算复杂的计算。
https://developer.nvidia.com/how-to-cuda-c-cpp
如果不是，您可能需要考虑多个线程进行计算。
无论如何，你能否象征性地提供计算。 如果有vector操作，你可能真的想要使用blas或lapack进行调查，在某些情况下可以并行运行操作。
It is conceivable (at the risk of being off-topic?) that you might be able to use python with numpy and/or scipy. To the extent that it was possible, your calculations might be more readable.
As you explicitly asked about high level optimizations, it might be worth trying different C++ compilers. Nowadays, compilers are very complex optimization beasts and CPU vendors might implement very powerful and specific optimizations. But please note, some of them are not free (but there might be a free academic program).
- GNU compiler collection is free, flexible, and available on many architectures
- Intel compilers are very fast, very expensive, and may also produce good results for AMD architectures (I believe there is an academic program)
- Clang compilers are fast, free, and might produce similar results to GCC (some people say they are faster, better, but this might differ for each application case, I suggest to make your own experiences)
- PGI (Portland Group) is not free as the Intel compilers.
- PathScale compilers might perform good results on AMD architectures
I've seen code snippets differ in execution speed by the factor of 2, only by changing the compiler (with full optimizations of course). But be aware of checking the identity of the output. To aggressive optimization might lead to different output, which is something you definitely want to avoid.
祝你好运！
