Расчет с плавающей точкой Пауэрс (PHP/BCMath)

Question

Я пишу оболочку для расширения bcmath и bug #10116 относительно bcpow() особенно раздражает - это бросает $right_operand ($exp) к (родной PHP, а не произвольной длины) целое число, поэтому при попытке вычислить квадратный корень (или любой другой корень выше, чем 1) числа вы всегда в конечном итоге с 1 вместо правильного результата.

Я начал искать алгоритмы, которые позволили бы мне вычислить п-й корень из числа, и я found this answer, который выглядит довольно твердый, я на самом деле expanded the formula использованием WolframAlpha, и я был в состоянии улучшить его скорость примерно 5%, сохраняя при этом точность результатов.

Вот чистая реализация PHP имитируя мою реализацию BCMath и свои ограничения:

function _pow($n, $exp) 
{ 
    $result = pow($n, intval($exp)); // bcmath casts $exp to (int) 

    if (fmod($exp, 1) > 0) // does $exp have a fracional part higher than 0? 
    { 
     $exp = 1/fmod($exp, 1); // convert the modulo into a root (2.5 -> 1/0.5 = 2) 

     $x = 1; 
     $y = (($n * _pow($x, 1 - $exp))/$exp) - ($x/$exp) + $x; 

     do 
     { 
      $x = $y; 
      $y = (($n * _pow($x, 1 - $exp))/$exp) - ($x/$exp) + $x; 
     } while ($x > $y); 

     return $result * $x; // 4^2.5 = 4^2 * 4^0.5 = 16 * 2 = 32 
    } 

    return $result; 
} 

выше seems to work greatкроме случаев, когда 1/fmod($exp, 1) не дает целое. Например, если $exp является 0.123456, его обратное будет 8.10005 и исход pow() и _pow() будет немного другой (demo):

• pow(2, 0.123456) = 1.0893412745953
• _pow(2, 0.123456) = 1.0905077326653
• _pow(2, 1/8) = _pow(2, 0.125) = 1.0905077326653

Как достичь такого же уровня точности, используя «ручной» экспоненциальный расчет?

Answer 1

Используемый алгоритм, чтобы найти п ^я корня из (положительных) чисел a является алгоритмом Ньютона для нахождения нуля

f(x) = x^n - a. 

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

Вычисление мощности с показателем 0 < y < 1, где y не в форме 1/n с целым числом n является более сложным. Делая аналог, решая

x^(1/y) - a == 0 

снова включать вычисление мощности с нецелым показателем, сама проблема, которую мы пытаемся решить.

Если y = n/d рациональна с малым знаменателем d, проблема легко решается путем вычисления

x^(n/d) = (x^n)^(1/d), 

но для наиболее рационального 0 < y < 1, числитель и знаменатель имеют довольно большой, и промежуточный x^n будет огромным, поэтому вычисление будет использовать много памяти и займет (относительно) долгое время. (Для примера показателя 0.123456 = 1929/15625, это не так уж плохо, но 0.1234567 будет довольно руление.)

Один из способов расчета мощности для общего рационального 0 < y < 1 это написать

y = 1/a ± 1/b ± 1/c ± ... ± 1/q 

с целыми числами a < b < c < ... < q и для умножения/деления индивидуума x^(1/k). (Каждое разумное 0 < y < 1 имеет такие представления, а самые короткие такие представления, как правило, не включают в себя многие термины, например

1929/15625 = 1/8 - 1/648 - 1/1265625; 

используя только дополнение в разложении приводит к более длинным представлениям с большими знаменателями, например

1929/15625 = 1/9 + 1/82 + 1/6678 + 1/46501020 + 1/2210396922562500, 

, так что потребуется больше работы.)

Некоторое улучшение возможно путем смешивания подходов, сначала найдите близкое рациональное приближение к y с малым знаменателем v ia разложение продолжительной фракции y - для примера показатель 1929/15625 = [0;8,9,1,192] и с использованием первых четырех частичных коэффициентов дает приближение 10/81 = 0.123456790123... [обратите внимание, что 10/81 = 1/8 - 1/648, частичные суммы кратчайшего разложения на чистые фракции являются сходящимися], а затем разлагают остаток в чистую фракции.

Однако в целом этот подход приводит к вычислению п ^я корней при больших n, который также является медленной и памятью интенсивно, если требуемая точность конечного результата высока.

В общем, это, вероятно, проще и быстрее реализовать exp и log и использовать

x^y = exp(y*log(x))