Ответ 1
Редактировать резюме
- В моем первоначальном ответе мы отметили, что код содержит много реплицированных вычислений и что многие из факторов задействованы в коэффициентах 1/3. Например,
pow(x, 0.1e1/0.3e1)совпадает сcbrt(x). - Мое второе редактирование было просто неправильным, и мой третий экстраполирован на эту ошибку. Именно это заставляет людей бояться изменять результаты, подобные оракулу, из символических математических программ, начинающихся с буквы "М". Я удалил (т.е.
strike) эти изменения и подтолкнул их к нижней части текущей ревизии этого ответа. Однако я их не удалял. Я человек. Нам легко ошибиться. - Мое четвертое редактирование разработало очень компактное выражение, которое корректно представляет свернутое выражение в вопросе IF параметры
l1,l2иl3- положительные действительные числа, а еслиa- ненулевое вещественное число. (Нам еще предстоит услышать от ОП относительно специфики этих коэффициентов. Учитывая характер проблемы, это разумные предположения.) - Это редактирование пытается ответить на общую проблему упрощения этих выражений.
Сначала первые вещи
Я использую Maple для генерации кода на С++, чтобы избежать ошибок.
Клен и Математика иногда пропускают очевидное. Еще важнее то, что пользователи Maple и Mathematica иногда ошибаются. Подставляя "часто" или, возможно, даже "почти всегда", вместо "иногда, вероятно, ближе к знаку".
Вы могли бы помочь Maple упростить это выражение, рассказав ему о соответствующих параметрах. В приведенном примере я подозреваю, что l1, l2 и l3 являются положительными действительными числами и что a является ненулевым вещественным числом. Если это так, скажите это. Эти символические математические программы обычно предполагают, что имеющиеся величины являются сложными. Ограничение домена позволяет программе делать допущения, которые недопустимы в комплексных числах.
Как упростить эти большие беспорядки из символических математических программ (это редактирование)
Символьные математические программы обычно предоставляют возможность предоставлять информацию о различных параметрах. Используйте эту способность, особенно если ваша проблема связана с делением или возведением в степень. В приведенном примере вы могли бы помочь Maple упростить это выражение, сообщив ему, что l1, l2 и l3 являются положительными действительными числами и что a является ненулевым вещественным числом. Если это так, скажите это. Эти символические математические программы обычно предполагают, что имеющиеся величины являются сложными. Ограничение домена позволяет программе делать предположения, такие как x b x= (ab) x. Это только если a и b - положительные действительные числа, а если x вещественное. Он недействителен в комплексных числах.
В конечном счете, эти символические математические программы следуют алгоритмам. Помогите ему. Попробуйте играть с расширением, сбором и упрощением перед созданием кода. В этом случае вы могли бы собрать эти термины с коэффициентом mu и с коэффициентом K. Сокращение выражения до его "простейшей формы" остается немного искусством.
Когда вы получаете уродливый беспорядок сгенерированного кода, не принимайте его как истину, которую вы не должны касаться. Попытайтесь упростить это самостоятельно. Посмотрите, что символическая математическая программа имела перед собой сгенерированный код. Посмотрите, как я уменьшил ваше выражение на что-то намного более простое и намного более быстрое, и как ответ Уолтера занял несколько шагов дальше. Нет волшебного рецепта. Если бы был волшебный рецепт, Клен применил бы его и дал бы ответ, который дал Уолтер.
Об определенном вопросе
В этом расчете вы делаете много сложения и вычитания. У вас могут быть серьезные проблемы, если у вас есть условия, которые почти отменяют друг друга. Вы тратите много 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.
Клен пропустил очевидное.
Например, гораздо проще писать
(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);
Неверный ответ, намеренно сдержанный
Обратите внимание, что это поражено. Это неправильно.
Обновление
Клен пропустил очевидное. Например, гораздо проще писать
(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 - положительные действительные числа.
Забастовкa >