2013-03-09 2 views
2

Я пишу программу для решения плоской ограниченной проблемы с тремя телами. Его уравнения ниже. Эта функция вычисляет производные по положению и скорости и записывает их в массив.Есть ли быстрый алгоритм, вычисляющий полномочия, умноженные на половину?

valarray<double> force(double t, valarray<double> r) 
{ 
    valarray<double> f(dim); 
    valarray<double>r0(r-rb0); 
    valarray<double>r1(r-rb1);  

    f[0]= 2 * r[1] + r[2] - (1 - mu)*r0[2]/norm3(r0) - mu*r1[2]/norm3(r1); 
    f[1]= - 2 * r[0] + r[3] - mu*r0[3]/norm3(r0) - mu*r1[3]/norm3(r1); 
    f[2] = r[0]; 
    f[3] = r[1]; 
    return f; 
} 

double norm3(valarray<double> x) 
{ 
    return pow(x[2]*x[2]+x[3]*x[3],1.5); 
} 

Поэтому я должен вычислить квадрат вектора положения, а затем поднять его до степени 3/2. Я думаю, что эти операции занимают большую часть времени вычислений.

Теперь я использую функцию pow из math.h. Есть ли еще более быстрый алгоритм для вычисления этой мощности? Я попытался использовать fast inverse square root (и куб его позже), но он дает слишком неточное значение для моих целей и работает дольше (возможно, из-за кубирования).

Спасибо!

+0

Я думаю, что часто используются две конкретные реализации: 1. экспоненциальное тождество ('x^1.5 = exp (ln (x) * 1.5)') и использование [многочлена Тейлора.] (Http: //en.wikipedia .org/wiki/Taylor's_theorem) – 2013-03-09 09:35:08

+1

как вы «кубировали его позже»? используя 'pow' или' x * x * x'? – stefan

+1

Почему вы не можете использовать t = x [2] * x [2] + x [3] * x [3]; return t * sqrt (t)? –

ответ

5

Простым подходом может быть попытка x * sqrt (x), но проверить его, чтобы быть уверенным.

double norm3(valarray<double> x) 
{ 
    double result=x[2]*x[2]+x[3]*x[3]; 
    result=result * sqrt(result); 
    return result; 
} 
1

The FSQRT в семействе процессоров 15h AMD64 занимает 52 циклов. Варианты SSE2 принимают 29 для скалярного значения и 38 для упакованной операции. Версия C sqrt(), вероятно, несколько дополнительных инструкций, но я сомневаюсь, что это намного больше.

Если вы хотите относительно точные результаты, я сомневаюсь, что с какой-либо другой операцией будет намного лучше. Скорее всего, что-либо, производящее хорошую точность с участием pow(), exp() и log() и т. Д., Займет больше времени.

Однако, спрашивая в Интернете, вы не побиваете свои собственные контрольные показатели. Если это критическая часть вашего кода, попробуйте несколько вариантов и посмотрите, что вы получаете.