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

Ответ 1

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

Ответ 2

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

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

metamarist обновляет методику битвыборки эпохи DEC-VAX W. Kahan с этим фрагментом, который "предполагает 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 раз больше стоимости умножения. Поэтому мы должны быть скупыми с нашими операциями деления.

Ответ 3

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;
}

Ответ 4

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

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

Здесь приведена пересмотренная версия, которая работает с числами между -/+ 1 quadrillion (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 мс на среднем потребительском оборудовании. Если бы скорость имела первостепенное значение, я бы подумал об использовании предварительно вычисленных таблиц поиска. Но это происходит от разработчика рабочего стола, а не от встроенного системного инженера.