Почему Math.cbrt(1728) дает более точный результат, чем Math.pow(1728, 1/3)?

Ответ 1

Несколько общих замечаний впереди:

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

• Стандарты для компьютерных языков обычно не гарантируют особая точность математических функций или одинаковые границы ошибок между различными математическими функциями, такими как cbrt() или pow(). Но математические библиотеки, которые обеспечивают корректно округленные результаты для данного точность существует, например CRlibm.

В этом случае, однако, cbrt(x) будет давать более точные результаты, чем pow(x,1.0/3.0), даже если обе функции правильно округлены для всех входов.

Проблема заключается в том, что 1.0/3.0 невозможно представить точно как число с плавающей запятой, будь то в двоичном или десятичном формате. Число двойной точности IEEE-754, наиболее близкое к одной трети, равно 3.3333333333333331e-1 (или 0x1.5555555555555p-2, если оно выражено в шестнадцатеричном формате с плавающей запятой C/С++). Относительная репрезентативная ошибка: -5.5511151231257827e-17 (-0x1.0000000000000p-54), что означает, что наилучшее представление с двойной точностью 1/3 несколько меньше желаемого математического значения.

Эта начальная ошибка на одном из входов pow() не только передается на выход, но и увеличивается из-за свойства увеличения ошибки при экспонировании. В результате pow(x,1.0/3.0), как правило, будет давать результаты, которые слишком малы по сравнению с желаемым корнем куба, даже если pow() обеспечивает корректно округленные результаты. Для примера в вопросе корректно округленные результаты

cbrt(1728.0)        = 1.2000000000000000e+1  (0x1.8000000000000p+3)
pow(1728.0,1.0/3.0) = 1.1999999999999998e+1  (0x1.7ffffffffffffp+3)

то есть результат от pow() равен ulp меньше результата из cbrt(). Для аргументов большой величины различие будет намного больше. Например, если x равно 2 ¹⁰²² соответствующие результаты отличаются на 94 ulps:

x              = 4.4942328371557898e+307  (0x1.0000000000000p+1022)
cbrt(x)        = 3.5553731598732904e+102  (0x1.965fea53d6e3dp+340)
pow(x,1.0/3.0) = 3.5553731598732436e+102  (0x1.965fea53d6ddfp+340)

Относительная погрешность результата pow() в этом примере равна 1.3108e-14, демонстрируя увеличение относительной ошибки, упомянутой выше.

По соображениям точности и производительности, математические библиотеки, которые реализуют cbrt(), поэтому обычно не отображают cbrt(x) в pow(x,1.0/3.0), а используют альтернативные вычислительные схемы. В то время как реализации будут отличаться, обычно используемый подход должен начинаться с начального приближения с низкой точностью, за которым следует один или несколько шагов метода Halley, который имеет кубическую конвергенцию.

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