Эффективное заполнение итерации Ньютона для корня куба

Question

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

Answer 1

Это обманчиво сложный вопрос. Здесь - хороший обзор некоторых возможных подходов.

Answer 2

Ввиду "гнили ссылок", которая настигла Принятый ответ, я дам более самостоятельный ответ, сосредоточенный на теме быстрого получения первоначального предположения, подходящего для суперлинейной итерации.

«Опрос», проведенный метамеристом ( Wayback link ), предоставил некоторые временные сравнения для различных комбинаций начального значения / итерации (включая методы Ньютона и Галлея). Его ссылки на работы W. Кахан , «Вычисление настоящего кубического корня» и К. Турковский , «Вычисление корня куба».

metamarist обновляет технику поиска битов эпохи DEC-VAX У. Кахана с помощью этого фрагмента, который «предполагает 32-разрядные целые числа» и использует удвоенный формат IEEE 754 для «генерации начальных оценок с точностью 5 битов»:

inline double cbrt_5d(double d) 
{ 
   const unsigned int B1 = 715094163; 
   double t = 0.0; 
   unsigned int* pt = (unsigned int*) &t; 
   unsigned int* px = (unsigned int*) &d; 
   pt[1]=px[1]/3+B1; 
   return t; 
} 

Код К. Турковского обеспечивает немного большую точность («приблизительно 6 бит») с помощью обычного масштабирования степеней двух на float fr, за которым следует квадратичное приближение к его корню куба за интервал [0,125,1,0) :

/* Compute seed with a quadratic qpproximation */
fr = (-0.46946116F * fr + 1.072302F) * fr + 0.3812513F;/* 0.5<=fr<1 */

и последующее восстановление показателя степени два (с поправкой на одну треть). Для извлечения и восстановления экспоненты / мантиссы используются математические вызовы из до frexp и ldexp.

Сравнение с другими приближениями "начального" корневого куба

Чтобы оценить эти приближения кубического корня, нам нужно сравнить их с другими возможными формами. Сначала критерии для оценки: мы рассматриваем аппроксимацию на интервале [1 / 8,1] и используем наилучшую (минимизирующую максимальную) относительную погрешность.

То есть, если f(x) является предлагаемым приближением к x^{1/3}, мы находим его относительную ошибку:

        error_rel = max | f(x)/x^(1/3) - 1 | on [1/8,1]

Самым простым приближением, конечно, было бы использование одной константы на интервале, и наилучшая относительная ошибка в этом случае достигается путем выбора f_0(x) = sqrt(2)/2, среднего геометрического значения в конечных точках. Это дает 1,27-битную относительную точность, быструю, но грязную отправную точку для итерации Ньютона.

Лучшим приближением будет лучший полином первой степени:

 f_1(x) = 0.6042181313*x + 0.4531635984

Это дает 4,12 бита относительной точности, большое улучшение, но меньше 5-6 битов относительной точности, обещанных соответствующими методами Кахана и Турковского. Но он находится в поле и использует только одно умножение (и одно сложение).

Наконец, что если мы позволим себе деление вместо умножения? Получается, что с одним делением и двумя «сложениями» мы можем иметь лучшую линейно-дробную функцию:

 f_M(x) = 1.4774329094 - 0.8414323527/(x+0.7387320679) 

, что дает 7,265 бит относительной точности.

На первый взгляд это кажется привлекательным подходом, но старое эмпирическое правило заключалось в том, чтобы рассматривать стоимость деления FP как три умножения FP (и в основном игнорировать сложения и вычитания). Однако в современных конструкциях FPU это нереально. Хотя относительная стоимость умножения на сложение / вычитание снизилась, в большинстве случаев до коэффициента два или даже равного, стоимость деления не снизилась, но часто в 7-10 раз превышала стоимость умножения. Поэтому мы должны быть скупы с нашими операциями подразделения.

Answer 3

Похоже, что вопрос оптимизации уже решен, но я хотел бы добавить улучшение к функции cubeRoot (), опубликованной здесь, для других людей, спотыкающихся на этой странице, ищущих быстрый алгоритм корня куба.

Существующий алгоритм работает хорошо, но за пределами диапазона 0-100 он дает неверные результаты.

Вот пересмотренная версия, которая работает с числами от - / + 1 квадриллиона (1E15).Если вам нужно работать с большими числами, просто используйте больше итераций.

static double cubeRoot( double num ){
    boolean neg = ( num < 0 );
    double x = Math.abs( num );
    for( int i = 0, iterations = 60; i < iterations; i++ ){
        x = ( ( 2 * x * x * x ) + num ) / ( 3 * x * x ); 
    }
    if( neg ){ return 0 - x; }
    return x;
}

Что касается оптимизации, я предполагаю, что оригинальный постер спрашивал, как предсказать минимальное количество итераций для точного результата, учитываяпроизвольный размер ввода.Но, похоже, что в большинстве общих случаев выигрыш от оптимизации не стоит дополнительной сложности.Даже при использовании описанной выше функции 100 итераций занимают менее 0,2 мс на среднем потребительском оборудовании.Если бы скорость была крайне важна, я бы подумал об использовании предварительно вычисленных таблиц поиска.Но это исходит от разработчика настольных систем, а не от инженера по встраиваемым системам.

Answer 4

static double cubeRoot(double num) {
    double x = num;

    if(num >= 0) {
        for(int i = 0; i < 10 ; i++) {
            x = ((2 * x * x * x) + num ) / (3 * x * x); 
        }
    } 
    return x;
}