Случай с числовой стабильностью - Алгоритмы

Есть два отрезка на плоскости, (X1,Y1)->(X2,Y2) и (X3,Y3)->(X4,Y4). Собственно, нас интересует точка (X3,Y3). Надо на первом отрезке вычислить координату X, соответствующую координате Y3.

Ну это элементарная пропорция:

K = (X2-X1) / (Y2-Y1)
X = (Y3-Y1) * K + X1

Проблема в том, что вычисления происходят на пределе числовой стабильности. Отрезки очень длинные, а расстояние (X2,Y2)->(X3,Y3) очень маленькое. В результате, вычисленная координата X может оказаться больше X3, хотя абсолютной точности представления в районе X,X3 все еще вполне достаточно.
Как бы так преобразовать выражение, учитывая специфику операций с плавающей точкой, чтобы на данном рисунке вычисленная X никогда не получалась больше X3?

Здравствуйте, McSeem2, Вы писали:

MS>Есть два отрезка на плоскости, (X1,Y1)->(X2,Y2) и (X3,Y3)->(X4,Y4).
А можно получить несколько тестовых примеров? Мне так шаманить проще...

... << RSDN@Home 1.1.4 beta 6a rev. 436>>

Здравствуйте, WolfHound, Вы писали:

MS>>Есть два отрезка на плоскости, (X1,Y1)->(X2,Y2) и (X3,Y3)->(X4,Y4).
WH>А можно получить несколько тестовых примеров? Мне так шаманить проще...

Это нелегко — оптимизатор мешается под ногами, так что его надо отключить. Иначе он все запихивает в регистры и считает в double. Понятно, что на реальном коде у него нет таких возможностей для оптимизации в силу того, что промежукточные значения хранятся в контейнерах, и типы там — float.

double xa1 = 1505707765.9675;
double ya1 = -4634090802.2879;
double xa2 = -1505707161.2032;
double ya2 = 4634091378.1953;
double xa3 = 325.11308288574;
double ya3 = 265.22279167175;

float x1 = (float)xa1; 
float y1 = (float)ya1; 
float x2 = (float)xa2; 
float y2 = (float)ya2; 
float x3 = (float)xa3; 
float y3 = (float)ya3; 
float dx = x2-x1;
float dy = y2-y1;
float k  = dx / dy;
float x  = (y3-y1) * k + x1;

Здесь отрезок — это маленький кусочек x1,y1,x2,y2 (не xa1...!). Ось Y идет вниз. Зеленая точка — (x3,y3). Красная — (x,y3), которая по идее должна лежать на отрезке. Ну или как минимум быть левее зеленой точки. На самом деле, с double все то же самое, но менее очевидно.

Здравствуйте, McSeem2, Вы писали:

MS>Здесь отрезок — это маленький кусочек x1,y1,x2,y2 (не xa1...!). Ось Y идет вниз. Зеленая точка — (x3,y3). Красная — (x,y3), которая по идее должна лежать на отрезке. Ну или как минимум быть левее зеленой точки. На самом деле, с double все то же самое, но менее очевидно.
Ну а чего ты хотел? Монтиса в float'е 23 бита десятичный логорифм 2^23 = 6,923689900 те 6 или с большой натяжкой 7 десятичных цифр, а у тебя их 14.
Для таких чисел придется использовать double ибо во float они не влезают физически.

... << RSDN@Home 1.1.4 beta 6a rev. 436>>

Во-первых, согласно приведенным данным левый отрезок должен иметь наклон примерно -1, т.е. в противоположную сторону от того, как на рисунке. Поэтому точка х может оказаться правее. Советую решить уравнение на пересечение прямых для проверки этого факта -- где будет лежать корень.
Плюс слишком много значащих цифр в данных. В float -- макс 7, в double -- 15.

Здравствуйте, McSeem2, Вы писали:

MS>Как бы так преобразовать выражение, учитывая специфику операций с плавающей точкой, чтобы на данном рисунке вычисленная X никогда не получалась больше X3?
На практике до таких задач не добрался. Из академических работ очень нравилась Practical Segment Intersection with Finite Precision Output. Пошарил здесь.
Еще несколько работ есть здесь.

Здравствуйте, WolfHound, Вы писали:

WH>Ну а чего ты хотел? Монтиса в float'е 23 бита десятичный логорифм 2^23 = 6,923689900 те 6 или с большой натяжкой 7 десятичных цифр, а у тебя их 14.
WH>Для таких чисел придется использовать double ибо во float они не влезают физически.

Ну и что, что не влазят. Пусть координаты концов представлены грубо — это не имеет значения. Все что я хочу — это "положить" вычисленную точку на отрезок с той точностью, которая существует в районе значений точки, а не концов отрезка. А точность там — хорошая. Но типа double у меня нет. Как быть?

Здравствуйте, FedorovStanislav, Вы писали:

FS>Во-первых, согласно приведенным данным левый отрезок должен иметь наклон примерно -1, т.е. в противоположную сторону от того, как на рисунке. Поэтому точка х может оказаться правее. Советую решить уравнение на пересечение прямых для проверки этого факта -- где будет лежать корень.

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

Здравствуйте, McSeem2, Вы писали:

MS>Отрезок здесь вообще не при чем. Возможно я сбил с толку этим вторым отрезком. Существует только левый отрезок и точка (X3,Y3). Посему, никаких точек пересечения прямых искать не надо — надо решить пропорцию.

Нет. Просто эта точка может лежать левее этого отрезка. Ты необоснованно утверждаешь обратное. Я понимаю, что для вычислений нужно знать только ординату y3, но потом ты результат сверяешь с x3 -- а оно то любым может быть, а ты говоришь, должно быть левее. С чего вдруг? Нарисуй оба отрезка в каком-нибуть тулсе.

Здравствуйте, FedorovStanislav, Вы писали:

FS>Нет. Просто эта точка может лежать левее этого отрезка. Ты необоснованно утверждаешь обратное. Я понимаю, что для вычислений нужно знать только ординату y3, но потом ты результат сверяешь с x3 -- а оно то любым может быть, а ты говоришь, должно быть левее. С чего вдруг? Нарисуй оба отрезка в каком-нибуть тулсе.

На самом деле именно это я и хочу определить. И точка 3 находится однозначно правее (проверено с double). Частичнм решением оказался выбор — с какого конца отсчитывать пропорцию. К какому концу отрезка y3 ближе, от того и отсчитывать. Но это не решает проблему, если y3 у нас примерно посередине.

Вот такой метод выдает правильный результат:

float x1 = (float)xa1; 
float y1 = (float)ya1; 
float x2 = (float)xa2; 
float y2 = (float)ya2; 
float x3 = (float)xa3; 
float y3 = (float)ya3; 

float xb1 = x1; 
float yb1 = y1; 
float xb2 = x2; 
float yb2 = y2; 
float x, y;
while((yb2 - yb1) > 0.001)
{
    x = xb1 / 2 + xb2 / 2;
    y = yb1 / 2 + yb2 / 2;
    if(y3 > y)
    {
        xb1 = x;
        yb1 = y;
    }
    else
    {
        xb2 = x;
        yb2 = y;
    }
}

Здесь отрезок продолжает "скакать" с той дискретностью, которая существует на концах отрезка (иначе и быть не может), но вычисленная точка всегда попадает на отрезок с точностью не хуже 0.001 (если это 0.001 не больше самой дикретности). Конечно, этот цикл будет жутко тормозить, но он доказывает, что задача принципиально решаема.

Здравствуйте, McSeem2, Вы писали:

MS>Есть два отрезка на плоскости, (X1,Y1)->(X2,Y2) и (X3,Y3)->(X4,Y4). Собственно, нас интересует точка (X3,Y3). Надо на первом отрезке вычислить координату X, соответствующую координате Y3.
MS>

MS>Ну это элементарная пропорция:
MS>

MS>K = (X2-X1) / (Y2-Y1)
MS>X = (Y3-Y1) * K + X1
MS>

MS>Проблема в том, что вычисления происходят на пределе числовой стабильности. Отрезки очень длинные, а расстояние (X2,Y2)->(X3,Y3) очень маленькое. В результате, вычисленная координата X может оказаться больше X3, хотя абсолютной точности представления в районе X,X3 все еще вполне достаточно.
MS>Как бы так преобразовать выражение, учитывая специфику операций с плавающей точкой, чтобы на данном рисунке вычисленная X никогда не получалась больше X3?

А не поможет здесь известный приём "сначала умножать, потом делить"? Т.е.: X = X1 + ( (Y3-Y1)*(X2-X1) ) / (Y2-Y1);

MS>Ну это элементарная пропорция:
MS>

MS>K = (X2-X1) / (Y2-Y1)
MS>X = (Y3-Y1) * K + X1
MS>

double ratio(double a,double b,double x) { return (x-a)/(b-a); }

double line(double a,double b,double t) { return b*t+a*(1-t); }

line(X1,X2,ratio(Y1,Y2,Y3));

В первую функуию можно добавить ещё обрезание по [0,1].

Здравствуйте, Шахтер, Вы писали:

Ш>

Ш>double ratio(double a,double b,double x) { return (x-a)/(b-a); }

Ш>double line(double a,double b,double t) { return b*t+a*(1-t); }

Ш>line(X1,X2,ratio(Y1,Y2,Y3));
Ш>

Попробовал, все то же самое, увы. Дополнительная погрешность возникает в операциях сложения и вычитания.

double xa1 = 1505707765.9675;
double ya1 = -4634090802.2879;
double xa2 = -1505707161.2032;
double ya2 = 4634091378.1953;
double xa3 = 325.11308288574;
double ya3 = 265.22279167175;

float x1 = (float)xa1; 
float y1 = (float)ya1; 
float x2 = (float)xa2; 
float y2 = (float)ya2; 
float x3 = (float)xa3; 
float y3 = (float)ya3; 
float k  = (x2-x1) / (y2-y1);

float xb1 = x1; 
float yb1 = y1; 
float xb2 = x2; 
float yb2 = y2; 
float x, y;

while((yb2 - yb1) > 1000)
{
    x = xb1 / 2 + xb2 / 2;
    y = yb1 / 2 + yb2 / 2;
    if(y3 > y)
    {
        xb1 = x;
        yb1 = y;
    }
    else
    {
        xb2 = x;
        yb2 = y;
    }
}

x = (y3-yb1) * k + xb1;

Здесь цикл последовательно приближает точки отрезка к искомой методом деления пополам. После того, как мы уложились в рамки значений не более 1000, выражение работает вполне точно. Если же закоментировать цикл, то части "y3-yb1" и "... + xb1" дают дополнительную погрешность. Раскрытие скобок, "y3*k — yb1*k" не помогает. Но при этом, вычисление тангенса угла наклона, "k=(x2-x1)/(y2-y1)" или же "k=(xb2-xb1)/(yb2-yb1)" не влияет на результирующую точность.
Я опишу проблему более подробно в отдельной ветке.

Здравствуйте, McSeem2, Вы писали:

Все что я хочу — это "положить" вычисленную точку на отрезок с той точностью, которая существует в районе значений точки, а не концов отрезка. А точность там — хорошая. Но типа double у меня нет. Как быть?

float y = y3;
float a = y / (y1 - y2) - y2 / (y1 - y2);
float x = a * x1 + (1 - a)*x2;

На вашем примере (выше по ветке) отклонение -5.1

... << RSDN@Home 1.1.4 stable SR1 rev. 568>>

Здравствуйте, What, Вы писали:

W>

W>float y = y3;
W>float a = y / (y1 - y2) - y2 / (y1 - y2);
W>float x = a * x1 + (1 - a)*x2;
W>

W>На вашем примере (выше по ветке) отклонение -5.1

Спасибо, конечно, но подобный вариант уже предлагал Шахтер. Порядок значений в окрестностях точки — сотни. Отклонение порядка 5 — это в данном диапазоне далеко за пределами того, что хотелось бы иметь. То есть, там, где значения исчисляются миллиардами — это нормально (близко к концам), но посередине хотелось бы выбрать всю возможную точность. Метод последовательных приближений это позволяет сделать. Должен быть и прямой способ — я чувствую

Здравствуйте, McSeem2, Вы писали:

MS>Спасибо, конечно, но подобный вариант уже предлагал Шахтер. Порядок значений в окрестностях точки — сотни. Отклонение порядка 5 — это в данном диапазоне далеко за пределами того, что хотелось бы иметь. То есть, там, где значения исчисляются миллиардами — это нормально (близко к концам), но посередине хотелось бы выбрать всю возможную точность. Метод последовательных приближений это позволяет сделать. Должен быть и прямой способ — я чувствую

Я думаю, прямого способа нет.
Давайте проанализируем, где теряется точность.
При сложении чисел с плавающей точкой с одинаковой экспонентой абсолютная погрешность результата равна единице в последнем знаке мантиссы, умноженной на экспоненту. Если результат примерно равен нулю, то в итоге получается чудовищная относительная погрешность. Поэтому огромная погрешность возникает просто при подсчёте середины отрезка, когда его концы далеко, а середина в районе нуля:

#include <iostream>

double xa1 = 1505707765.9675;
double ya1 = -4634090802.2879;
double xa2 = -1505707161.2032;
double ya2 = 4634091378.1953;
double xa3 = (xa1 + xa2) / 2.;
double ya3 = (ya1 + ya2) / 2.;

float x1 = (float)xa1; 
float y1 = (float)ya1; 
float x2 = (float)xa2; 
float y2 = (float)ya2; 
float x3 = (x1 + x2) / 2.f;
float y3 = (y1 + y2) / 2.f;

int main()
{
    std::cout << (xa3 - x3) << "    " << (ya3 - y3) << std::endl;
    return 0;
}

возвращает

-17.6179    31.9537

... << RSDN@Home 1.1.4 stable SR1 rev. 568>>

W>Давайте проанализируем, где теряется точность.
W>При сложении чисел с плавающей точкой с одинаковой экспонентой абсолютная погрешность результата равна единице в последнем знаке мантиссы, умноженной на экспоненту. Если результат примерно равен нулю, то в итоге получается чудовищная относительная погрешность. Поэтому огромная погрешность возникает просто при подсчёте середины отрезка, когда его концы далеко, а середина в районе нуля.

Нет, проблема не в этом. Точность теряется еще при преобразовании из double во float. Пусть теряется — это неизбежно. А вот при нахождении середины — уже ничего не теряется. Проверим — возьмем снова double и присвоим им значения из float:

#include <iostream>

double xa1 = 1505707765.9675;
double ya1 = -4634090802.2879;
double xa2 = -1505707161.2032;
double ya2 = 4634091378.1953;
double xa3 = (xa1 + xa2) / 2.;
double ya3 = (ya1 + ya2) / 2.;

float x1 = (float)xa1; 
float y1 = (float)ya1; 
float x2 = (float)xa2; 
float y2 = (float)ya2; 
float x3 = (x1 + x2) / 2.f;
float y3 = (y1 + y2) / 2.f;

double xb1 = x1; 
double yb1 = y1; 
double xb2 = x2; 
double yb2 = y2; 
double xb3 = (xb1 + xb2) / 2.;
double yb3 = (yb1 + yb2) / 2.;

int main()
{
    std::cout << (xa3 - x3) << "    " << (ya3 - y3) << std::endl;
    std::cout << (xb3 - x3) << "    " << (yb3 - y3) << std::endl;
    return 0;
}

Выдает:

-17.6179    31.9537
0    0

Таким образом, при делении пополам никаких потерь не возникает. В данном примере это происходит потому, что суммы x1+x2, y1+y2 имеют маленькое значение по модулю, таким образом, абсолютное значение точности у нас выше, чем на концах. Но если оба x1 и x2 — положительны или отрицательны, операция сложения тоже может приводить к потерям. Но можно раскрыть скобки: x1/2+x2/2. Деление пополам — операция абсолютно точная (мантисса не трогается вообще), а при сложеним мы уже оперируем меньшими по модулю числами.

Задача заключается в том, как бы нам с такой же точностью найти произвольную пропорцию, а не только "пополам"?

Здравствуйте, McSeem2, Вы писали:

W>>При сложении чисел с плавающей точкой с одинаковой экспонентой абсолютная погрешность результата равна единице в последнем знаке мантиссы, умноженной на экспоненту. Если результат примерно равен нулю, то в итоге получается чудовищная относительная погрешность. Поэтому огромная погрешность возникает просто при подсчёте середины отрезка, когда его концы далеко, а середина в районе нуля.

MS>Нет, проблема не в этом. Точность теряется еще при преобразовании из double во float. Пусть теряется — это неизбежно. А вот при нахождении середины — уже ничего не теряется. Проверим — возьмем снова double и присвоим им значения из float:
MS>


[skip]

MS>{
MS>    std::cout << (xa3 - x3) << "    " << (ya3 - y3) << std::endl;
MS>    std::cout << (xb3 - x3) << "    " << (yb3 - y3) << std::endl;
MS>    return 0;
MS>}
MS>

MS>Выдает:
MS>

MS>-17.6179    31.9537
MS>0    0
MS>

MS>Таким образом, при делении пополам никаких потерь не возникает. В данном примере это происходит потому, что суммы x1+x2, y1+y2 имеют маленькое значение по модулю, таким образом, абсолютное значение точности у нас выше, чем на концах. Но если оба x1 и x2 — положительны или отрицательны, операция сложения тоже может приводить к потерям. Но можно раскрыть скобки: x1/2+x2/2. Деление пополам — операция абсолютно точная (мантисса не трогается вообще), а при сложеним мы уже оперируем меньшими по модулю числами.

Итого. Потерями точности при преобразовании double -> float в нашем случае можно пренебречь. Таким образом, точность теряется при вычислениях. Их там всего два — сложение и деление пополам. Действительно, при делении пополам потери точности нет. (Кстати, поэтому всё равно, как вычислять, (x1 + x2) / 2 или x1 / 2 + x2 / 2). А вот при сложении и теряется вся точность. Получается, что из-за того, что результат близок к нулю, а слагаемые достаточно большие по модулю, большая часть первых знаков мантиссы сокращаяется и в мантиссе результата большинство знаков неверные. Из-за этого и возникает огромная относительная погрешность (по сути, для чисел с плавающей точкой относительная погрешность эквивалентна количеству неверных знаков в мантиссе). А вот если бы оба слагаемых были одного знака, все знаки мантиссы остались бы верными, и тогда были бы лишь незначительные (относительно числа) потери точности, связанные с размером мантиссы.
Вот пример в десятичных цифрах того, что происходит при вычислении середины:

Здесь ? обозначены нули, которые являются неверными знаками.
double x1:       +1.14 * 10^2
double x2:       -1.10 * 10^2
double x1 + x2:  +4.?? * 10^0
float  x1:       +1.1  * 10^2
float  x2:       -1.1  * 10^2
float  x1 + x2:  +?.?? * 10^0

MS>Задача заключается в том, как бы нам с такой же точностью найти произвольную пропорцию, а не только "пополам"?
Я понимаю, что в задаче нужно найти произвольную пропорцию. Я просто говорю, что мы не можем достаточно точно решить во float'ах даже частный случай нахождения середины (разумеетя, я имею в виду прямым методом, а не итеративным). А лучшее, что можно получить прямым методом, имхо, это вариант вида x1 * a + x2 * (1 — a).
Более того, чем ближе значение выражения x1 * a + x2 * (1 — a) к нулю и чем дальше от нуля x1 и x2, тем больше будут потери точности.

... << RSDN@Home 1.1.4 stable SR1 rev. 568>>

Здравствуйте, What, Вы писали:

W>Итого. Потерями точности при преобразовании double -> float в нашем случае можно пренебречь. Таким образом, точность теряется при вычислениях. Их там всего два — сложение и деление пополам.

Нет, как раз точность теряется при преобразовании double -> float, что вполне естественно и ожидаемо. Все значения во Float — это "реальность, данная нам в ощущениях", то есть, другого у нас нет. Но нахождение среднего не сбоит — результат получается в точности такой же, как и в вычислениях с double, если мы берем исходные данные из float (xb1... взяты из x1...). То есть, формула (a+b)/2 в данном случае работает точно. Более того, мы можем выпролнять эти пополамные деления многократно без потери точности.

От:	McSeem2	http://www.antigrain.com
Дата:	04.02.06 17:16
Оценка:

	От:	WolfHound
	Дата:	04.02.06 19:00
	Оценка:

От:	McSeem2	http://www.antigrain.com
Дата:	04.02.06 20:47
Оценка:

	От:	WolfHound
	Дата:	04.02.06 21:48
	Оценка:

	От:	FedorovStanislav
	Дата:	04.02.06 22:31
	Оценка:

От:	FreshMeat	http://www.rsdn.org
Дата:	04.02.06 23:51
Оценка:

От:	WinterMute	http://yarrr.ru
Дата:	06.02.06 12:34
Оценка:

	От:	What
	Дата:	07.02.06 13:43
	Оценка:

	От:	What
	Дата:	07.02.06 22:06
	Оценка:

	От:	What
	Дата:	08.02.06 19:03
	Оценка: