Как правильно разделить крошечные числа с двойной точностью без ошибок точности?

Я пытаюсь диагностировать и исправить ошибку, которая сводится к тому, что X/Y дает нестабильный результат, когда X и Y малы:

введите здесь описание изображения

В этом случае и cx, и patharea плавно увеличиваются. Их отношение представляет собой гладкую асимптоту при больших числах, но неустойчивую для «маленьких» чисел. Очевидная первая мысль состоит в том, что мы достигаем предела точности с плавающей запятой, но сами числа далеки от него. Типы ActionScript «Числа» представляют собой числа с плавающей запятой двойной точности IEE 754, поэтому должны иметь точность 15 десятичных цифр (если я правильно понял).

Некоторые типичные значения знаменателя (patharea):

0.0000000002119123
0.0000000002137313
0.0000000002137313
0.0000000002155502
0.0000000002182787
0.0000000002200977
0.0000000002210072

И числитель (сх):

0.0000000922932995
0.0000000930474444
0.0000000930582124
0.0000000938123574
0.0000000950458711
0.0000000958000159
0.0000000962901528
0.0000000970442977
0.0000000977984426

Каждое из них монотонно возрастает, но соотношение, как видно выше, хаотично.

При больших числах она сводится к гладкой гиперболе.

Итак, мой вопрос: как правильно работать с очень маленькими числами, когда вам нужно разделить одно на другое?

Я думал заранее умножить числитель и/или знаменатель на 1000, но не смог.

Фактически рассматриваемый код — это функция recalculate() здесь. Он вычисляет центроид многоугольника, но когда полигон крошечный, центроид беспорядочно прыгает вокруг места и может оказаться на большом расстоянии от многоугольника. Приведенные выше ряды данных являются результатом перемещения одного узла многоугольника в постоянном направлении (вручную, поэтому он не идеально гладкий).

Это Adobe Flex 4.5.


person Steve Bennett    schedule 06.05.2012    source источник
comment
Что произошло, когда вы умножили на 1000, разделили, а затем разделили результат на 1000?   -  person K2xL    schedule 06.05.2012
comment
Ну, ничего хорошего :) В одном текущем воплощении я получил частное с точностью всего до двух знаков. В этот момент мне казалось, что я на самом деле кодирую методом проб и ошибок, поэтому моя цель — получить немного знаний о том, как делать что-то правильно.   -  person Steve Bennett    schedule 06.05.2012
comment
Кроме того, вы бы умножили числитель или знаменатель? Я полагаю, последнее?   -  person Steve Bennett    schedule 06.05.2012


Ответы (2)


Я считаю, что проблема, скорее всего, вызвана следующей строкой в ​​​​вашем коде:

sc = (lx*latp-lon*ly)*paint.map.scalefactor;

Если ваш полигон очень маленький, то lx и lon почти одинаковы, как и ly и latp. Они оба очень большие по сравнению с результатом, поэтому вы вычитаете два числа, которые почти равны.

Чтобы обойти это, мы можем использовать тот факт, что:

x1*y2-x2*y1 = (x2+(x1-x2))*y2 - x2*(y2+(y1-y2))
            = x2*y2 + (x1-x2)*y2 - x2*y2 - x2*(y2-y1)
            = (x1-x2)*y2 - x2*(y2-y1)

Итак, попробуйте следующее:

dlon = lx - lon
dlat = ly - latp
sc = (dlon*latp-lon*dlat)*paint.map.scalefactor;

Значение математически то же самое, но члены на порядок меньше, поэтому ошибка также должна быть на порядок меньше.

person Jeffrey Sax    schedule 06.05.2012
comment
Гений. Абсолютно решает проблему. Вы заслуживаете большего, чем мой скудный плюс и галочка «правильный ответ». - person Steve Bennett; 06.05.2012
comment
Кстати, не могли бы вы (или кто-то) уточнить, что вы подразумеваете под ошибкой? Дело в том, что, уменьшая разницу в величине между значениями и их различиями, мы каким-то образом лучше используем ограниченное пространство 64-битного числа с плавающей запятой? Это происходит на уровне компилятора или на уровне чипа? - person Steve Bennett; 09.05.2012
comment
О, на самом деле я понял. lx и latp составляют около 100. Перемножив их вместе, вы получите около 10 000 (две значащие цифры десятичной точности потеряны). Теперь я понимаю, почему умножение на 1000 ничего не меняет. - person Steve Bennett; 09.05.2012

Джеффри Сакс правильно определил основную проблему - потерю точности из-за объединения терминов, которые (намного) больше, чем окончательный результат. Предлагаемое переписывание устраняет часть проблемы — по-видимому, достаточно для реального случая, учитывая положительный ответ.

Однако вы можете обнаружить, что если полигон снова станет (намного) меньше и/или дальше от начала координат, неточность снова проявится. В переписанной формуле члены все еще немного больше, чем их разница.

Кроме того, в алгоритме есть еще одна проблема «объединения больших и сравнимых чисел с разными знаками». Различные значения 'sc' в последующих циклах итерации по краям многоугольника эффективно объединяются в конечное число, которое (намного) меньше, чем отдельные sc(i). (если у вас есть выпуклый многоугольник, вы обнаружите, что существует одна непрерывная последовательность положительных значений и одна непрерывная последовательность отрицательных значений, в невыпуклых многоугольниках отрицательные и положительные значения могут быть переплетены).

Фактически алгоритм вычисляет площадь многоугольника путем добавления площадей треугольников, натянутых на ребра и начало координат, где некоторые члены отрицательны (всякий раз, когда ребро проходится по часовой стрелке, если рассматривать его из начала координат) и некоторый позитив (против часовой стрелки пройти через край).

Вы избавляетесь от ВСЕХ проблем с потерей точности, определяя начало координат в одном из углов многоугольника, скажем (lx,ly), а затем добавляя треугольные поверхности, натянутые на края и этот угол (итак: преобразование lon в ( lon-lx) и latp to (latp-ly) — с дополнительным бонусом, заключающимся в том, что вам нужно обрабатывать на два треугольника меньше, потому что, очевидно, ребра, которые ссылаются на выбранный исходный угол, дают нулевые поверхности.

Для области-части это все. Для центроидной части вам, конечно, придется «преобразовать» результат в исходный кадр, то есть добавить (lx, ly) в конце.

person Bert te Velde    schedule 04.09.2012
comment
Спасибо за это дополнительное объяснение. - person Steve Bennett; 05.09.2012