Расчет чисел с плавающей точкой (PHP/BCMath)

Question

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

Я начал искать алгоритмы, которые позволили бы мне вычислить n-й корень числа, и я нашел этот ответ, который выглядит довольно солидно, я фактически расширил формулу, используя 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;
}

Выше, кажется, прекрасно работает, кроме случаев, когда 1 / fmod($exp, 1) не дает целое число Например, если $exp является 0.123456, его обратное будет 8.10005 и результат pow() а также _pow() будет немного по-другому ( демо):

• pow(2, 0.123456) знак равно 1.0893412745953
• _pow(2, 0.123456) знак равно 1.0905077326653
• _pow(2, 1 / 8) знак равно _pow(2, 0.125) знак равно 1.0905077326653

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

Answer 1

Используемый алгоритм нахождения n-^го корня (положительного) числа 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 с небольшим знаменателем через продолжение расширения дроби y - для примера степени 1929/15625 = [0;8,9,1,192] и использование первых четырех частных отношений дает приближение 10/81 = 0.123456790123... [Обратите внимание, что 10/81 = 1/8 - 1/648частичные суммы кратчайшего разложения на чистые дроби являются сходящимися], а затем разлагают остаток на чистые дроби.

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

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

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