Вычисление расстояния между двумя точками в полярных координатах с использованием чисел с плавающей запятой

Когда координаты двух точек на плоскости представлены в полярной форме как r1, a1 и r2, a2, где r1, r2, a1, a2 — числа с плавающей запятой, и цель состоит в том, чтобы вычислить расстояние между двумя точками в виде числа с плавающей запятой, может возникнуть соблазн использовать следующую математическую формулу:

D = sqrt(r1*r1 + r2*r2 - 2*r1*r2*cos(a1-a2))

Эту формулу можно найти в этом и нескольких других ответах на StackOverflow. Ссылка на источник, указанный в связанном ответе, мертва, но вы найдете эту формулу во многих математических ресурсах, например этого веб-сайта.

Как формула с плавающей запятой, эта формула имеет ряд нежелательных свойств, перечисленных ниже в порядке уменьшения серьезности. Короче говоря, мой вопрос: какая формула лучше для вычисления расстояния D в условиях, указанных выше?

A/ sqrt применяется к отрицательному числу, когда r1 = r2 и r1*r1 субнормальны

Когда r1*r1 ниже нормы, а r1 = r2, округление r1*r1 отличается от округления 2*r1*r2. На одном полюсе этого спектра контрпримеров r1*r1 равно нулю, а 2*r1*r2 отлично от нуля. С другой стороны, r1*r1 является субнормальным и округляется несколько резко вниз, а 2*r1*r2 является нормальным и округляется менее резко. В любом случае выражение внутри sqrt оказывается отрицательным, а результатом формулы с плавающей запятой для D является NaN.

Примеры значений для двойной точности:

  double r1 = sqrt(DBL_MIN * DBL_EPSILON) * sqrt(0.45);
  double r2 = r1;

Запустите его в Compiler Explorer.

B/ sqrt применяется к отрицательному числу, когда r1 близко к r2, а отличается

Когда r1 и r2 очень близки, расстояние между математическими произведениями r1*r1 и r1*r2 очень близко к расстоянию между математическими произведениями r1*r2 и r2*r2. Когда это общее расстояние соответствует небольшому нечетному числу половинных ULP, может возникнуть ситуация, когда умножение с плавающей запятой r1*r1 округляется в меньшую сторону, r1*r2 округляется в большую сторону и r2*r2 снова округляется в меньшую сторону. В этих условиях, опять же, часто встречающаяся формула берет квадратный корень из отрицательного числа.

Примеры значений для двойной точности:

  double r1 = 0x1.1c71c71c71d3fp0;
  double r2 = 0x1.1c71c71c71dcbp0;

Запустите его в Compiler Explorer.

C/ Когда r1 близко к r2, вычитание увеличивает приближения, которые имеют место при умножении

На самом деле это более мягкий симптом тех же основных причин проблем A/ и B/. Когда r1 очень близко к r2, возникает явление, известное как катастрофическая отмена. Относительная точность вычитания ужасна (например, если a1 = a2 и r1 близки к r2, округление при умножении может сделать результат вычитания равным 0,0, хотя оптимальный ответ, fabs(r1 - r2), был представим и бесконечно более точен в относительных единицах. Примечание. что абсолютная точность результата порядка ULP r1 и r2, и это все еще может быть хорошо.

D/Переполнение, когда величина r1 или r2 слишком велика

Если переполняется только r1*r1 или r2*r2, результат вычисляется как +inf, что может быть не лучшим представимым приближением математического расстояния.

Если r1*r2 переполняется, то результатом вычитания является NaN, и поэтому расстояние вычисляется как NaN.

Проблемы A/ и B/ делают бессмысленным результат, которого не должно быть, и могут быть решены путем вычисления dq > 0 ? sqrt(dq) : 0 вместо dq. Для входных данных, вызвавших их, это изменение делает ответ 0.0 вместо этого. Этот результат имеет бесконечную относительную ошибку, как и результаты для других входных данных, потому что это не решает проблему C/.

Проблема D/ может быть решена путем масштабирования, если программист ожидает, что вычисление будет использоваться в условиях, которые его запускают. Если на то пошло, проблема A/ также может быть решена масштабированием, но это не решит проблему B/.

Вполне вероятно, что полное решение, которое решает все проблемы A-D, потребует больше вычислений. Может существовать несколько золотых точек, которые решают только некоторые проблемы или решают проблему C/ более или менее тщательно, вычисляя расстояния с точностью до 10 ULP, или 3, или 1. Любое решение, которое улучшается по сравнению с исходной точкой, является достойный ответа.

Гийом Мелькион уже указывал за пределами сайта, что приведенная ниже формула эквивалентна оригиналу в математических терминах и, очевидно, избегает вопросов A/ и B/, поскольку аргумент квадратного корня представляет собой сумму неотрицательных членов:

D = sqrt((r1-r2)*(r1-r2) + 2*r1*r2*(1 - cos(a2-a1)))

В этом решении катастрофическая отмена происходит в 1 - cos(a2-a1), поэтому некоторые аспекты проблемы С/ остаются (хотя расчет по этой формуле оптимален для a1=a2 и r1, близких к r2, так как тогда r1-r2 и cos(a2-a1) точны). Ситуация с проблемой D/ улучшилась, но остались случаи, когда результаты представлялись в виде конечного значения, а формула вычисляла +inf.