Стабильный котангенс
Существует ли более стабильная реализация для котангентной функции, чем return 1.0/tan (x);?
Ответы
Ответ 1
cot(x) = cos(x)/sin(x)
должна быть более численно устойчивой вблизи π/2, чем cot(x) = 1/tan(x)
. Вы можете эффективно реализовать sincos
на своих платформах.
Другая возможность - cot(x) = tan(M_PI_2 - x)
. Это должно быть быстрее, чем указано выше (даже если sincos
доступно), но оно также может быть менее точным, так как M_PI_2
- это, конечно, только приближение трансцендентного числа π/2, поэтому разность M_PI_2 - x
будет не быть точным для полной ширины a double
мантиссы - на самом деле, если вам не повезло, у нее могут быть только несколько значимых бит.
Ответ 2
Если вы рассмотрите угол между двумя векторами (v
и w
), вы также можете получить котангенс следующим образом (используя Eigen :: Vector3d):
inline double cot(Eigen::Vector3d v, Eigen::Vector3d w) {
return( v.dot(w) / (v.cross(w).norm()) );
};
При тета-угле между v
и w
выше функция верна, потому что:
- | VxW | = | v |. | w |.sin (тета)
- против w = | v |. | w |.cos (тета)
- детская кроватка (тета) = cos (тета)/грех (тета) = (против w)/| vxw |
Ответ 3
TL; DR No.
Как правило, при поиске источников неточности следует в первую очередь заботиться о дополнениях и вычитаниях, которые могут привести к проблеме вычитания. Умножения и деления, как правило, безвредны для точности, кроме добавления дополнительной ошибки округления, но могут вызвать проблемы из-за переполнения и недостаточного значения при промежуточных вычислениях.
Ни один номер машины x
может быть достаточно близко к кратным числам π/2, чтобы вызвать переполнение tan(x)
, поэтому tan(x)
является четко определенным и конечным для всех кодировок с плавающей запятой для любой из чисел с плавающей запятой IEEE-754 форматы, и, соответственно, так же, как cot(x) = 1.0/tan(x)
.
Это легко продемонстрировать, выполнив исчерпывающий тест со всеми числовыми кодировками с float
, поскольку исчерпывающий тест с использованием double
невозможен, за исключением, возможно, самых больших суперкомпьютеров, существующих сегодня.
Используя математическую библиотеку с точной реализацией tan()
с максимальной ошибкой ~ = 0,5 ulp, мы обнаруживаем, что при вычислении cot(x) = 1.0/tan(x)
максимальная ошибка составляет менее 1,5 ulp, где дополнительная Ошибка по сравнению с самой tan()
вносится ошибкой округления деления.
Повторяя этот исчерпывающий тест для всех значений с float
с помощью cot(x) = cos(x)/sin(x)
, где sin()
и cos()
вычисляются с максимальной ошибкой ~ = 0,5 ulp, мы находим, что максимальная ошибка у cot()
меньше 2,0 ulps, поэтому немного больше. Это легко объяснить наличием трех источников ошибок вместо двух в предыдущей формуле.
Наконец, cot(x) = tan (M_PI_2 - x)
страдает от проблемы вычитания, упомянутой ранее, когда x
около M_PI_2, и от проблемы, которая в арифметике с плавающей точкой конечной точности, M_PI_2 - x == M_PI_2
когда x
достаточно мала по величине. Это может привести к очень большим ошибкам, из-за которых мы не получим правильных битов в результате.