Числа с плавающей точкой

← →
Пономарев Андрей (2005-08-31 15:19) [0]

Добрый день, господа.
Стал свидетелем удивительного явления! :)
Имеем код:

Result := R * arccos(
sin(PI*one.Latitude/180)*sin(PI*two.Latitude/180) +
cos(PI*one.Latitude/180)*cos(PI*two.Latitude/180)*
cos(PI*(one.Longtitude-two.Longtitude)/180)
);

(по этой формуле вычисляется расстояние между двумя точками на поверхности сферы)
Latitude и Longtitude имеют тип Real,
а встроенные PI, cos, и sin - Extended.

При попытке обсчитать довольно большой массив точек (около 7тыс), встретилось 4 пары, для которых расчет формулы привел к возникновению исключения (аргумент arccos > 1). Удивлению не было предела. Пересчитал на windows-калькуляторе - получил 1 под arccos.

Погрешности вычислений? Да, Latitude и Longtitude "всего лишь" real, но после первого же умножения (на PI) тип становится extended - самое точное, что возможно в Delphi.

Что делать - я придумал сам, пусть кривовато, но сойдет. Мучает вопрос - ПОЧЕМУ???

Спасибо за внимание.

← →
Digitman © (2005-08-31 15:24) [1]

> встретилось 4 пары

приведи значения этих пар ..

← →
wal © (2005-08-31 15:39) [2]

IMHO на (PI*one.Latitude/180) компилятор мог сначала деление выполнить, а потом умножение, при этом дополнительную потерю точности имеем.

С уважением.

← →
Пономарев Андрей (2005-08-31 15:47) [3]

Сложно сказать. FloatToStr для первой пары дает:
(52.033332824707; 59.716667175293)
(52.033332824707; 59.716667175293)
(координаты действительно одинаковые)

не знаю насколько FloatToStr честен.

← →
Пономарев Андрей (2005-08-31 15:52) [4]

2 wal: не уверен, что это имеет значение... не могу сообразить... ведь даже если ошибка округления там появляется, то результат становится идентичен полученному с какой-то немного отличающейся latitude - ну и что? формула должна работать с любой!..

← →
Юрий Зотов © (2005-08-31 16:06) [5]

Определяем константу:
const K: extended = PI/180;
После чего просто умножаем на нее - ИМХО, будет и точнее, и быстрее. Еще лучше - сразу хранить данные в радианах, тогда вообще ничего не надо будет преобразовывать.

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

← →
wal © (2005-08-31 16:10) [6]

> [4] Пономарев Андрей (31.08.05 15:52)
Если в одном случае мы точность потеряем, а в другом нет (не уверен, что такое возможно, но все же), то в разных частях формулы появятся разные (пусть и не на много) числа, а по логике должны быть одинаковые. В итоге вполне можем получить под арккосинусом и >1. Опять же ИМХО. Еще могу порекомендовать отдебагить в окне CPU-FPU, и каждый (значимый) шаг на листочек записать с результатами.

С уважением.

← →
Digitman © (2005-08-31 16:16) [7]

ситуацию воспроизвести не удалось

type TRec = packed record Latitude, Longtitude: Real; end; var R: Double = 1.0; Rs: Double; one: TRec = (Latitude: 52.033332824707; Longtitude: 59.716667175293); two: TRec = (Latitude: 52.033332824707; Longtitude: 59.716667175293); procedure TForm1.Button1Click(Sender: TObject); begin Rs := R * arccos( sin(PI*one.Latitude/180)*sin(PI*two.Latitude/180) + cos(PI*one.Latitude/180)*cos(PI*two.Latitude/180)* cos(PI*(one.Longtitude-two.Longtitude)/180) ); end;

покажи как у тебя вычисляются значения координат ...

← →
Jeer © (2005-08-31 16:43) [8]

Неверно обозвана долгота.
Надо Longitude, тогда все получиться:))

← →
palva © (2005-08-31 17:06) [9]

Я бы перед вычислением arccos x вставил проверку на выход x за пределы диапазона [-1, 1]. И если выходит, то делать x = 1 или x = -1. Вернее, даже arccos не вычислять а полагать сразу 0 или pi.

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

← →
Пономарев Андрей (2005-08-31 20:59) [10]

Всем большое спасибо за участие. За что и люблю форумы на этом сайте. :)

2 Digiman:
То, что ситуацию смоделировать таким способом не удастся - почти не сомневался. Сделал все так же, как Вы, но на другой паре. Считает. Расписал real по байтам:
1) так получается в реальной программе:
Широта: 0 0 0 160 170 170 75 64
Долгота: 0 0 0 128 136 40 77 64
2) так получается, если забить результат FloatToStr руками в исходник:
Широта: 3 0 0 160 170 170 75 64
Долгота: 5 0 0 128 136 40 77 64
С "точки зрения" FloatToStr значение одно и то же! It"s a kind of magic!
В программе значения координат забираются у Interbasе"a.

2 Юрий Зотов:
Да, так, конечно, быстрее, и точнее... наверное... но ошибка все равно проявляется.

2 Jeer:
Сразу даже не поверил - уж больно привык к такому написанию. Слазил в Мюллера. Черт возьми, Вы правы! :) Спасибо! Одним заблуждением меньше.

2 palva:
Ну да, так и получается. Только я запихал вычисление в try...except. Так, конечно, медленнее получается, но меня устраивает.

Нда... печально. Как я понимаю, единственный универсальный совет в этом случае - "бди!".

2 Пономарев Андрей (31.08.05 20:59) [10]

>Нда... печально. Как я понимаю, единственный универсальный совет в этом случае - "бди!".

Ну, не знаю насколько печально, но данный совет универсален не только в этом случае!
Кстати в "Потрепаться" есть ветка по поводу этого самого "бди!"

> Расписал real по байтам

что-то я ума не приложу, откуда у тебя взялись такие значения байт ...

вот так это должно выглядеть для "проблемных" real-значений:

(расписано в порядке от старшего байта к младшему)

one.Latitude (FloatToStr() = 52.033332824707)
Hex: 40 4A 04 44 3F FF FF FC
Dec: 64 74 4 68 63 255 255 252

one.Longtitude (FloatToStr() = 59.716667175293)
Hex: 40 4D DB BB C0 00 00 04
Dec: 64 77 219 187 192 0 0 4

← →
Пономарев Андрей (2005-09-01 09:48) [13]

> Digitman © (01.09.05 08:58) [12]
Прошу прощения, что запутал, но я, кажется, оговорился, что пробовал на другой паре из числа проблемных.

uses SysUtils, Math; type Z = record case byte of 0 : (a : real); 1 : (b : array [0..sizeof(real)-1] of byte); end; rec = record Latitude : Z; Longtitude : Z; end; var one : rec = (Latitude:(a:55.3333320617676); Longtitude:(a:58.3166656494141)); two : rec = (Latitude:(a:55.3333320617676); Longtitude:(a:58.3166656494141)); procedure printRepresentation(value : Z); var i : integer; begin for i := 0 to sizeof(real)-1 do write(value.b[i]," "); end; procedure printAll; begin writeln(FloatToStr(one.Latitude.a),"; ",FloatToStr(one.Longtitude.a)); printRepresentation(one.Latitude); write("; "); printRepresentation(one.Longtitude); end; const PI_180 : extended = PI / 180; one_ : rec = (Latitude : (b : (0, 0, 0, 160, 170, 170, 75, 64)); Longtitude : (b : (0, 0, 0, 128, 136, 40, 77, 64))); two_ : rec = (Latitude : (b : (0, 0, 0, 160, 170, 170, 75, 64)); Longtitude : (b : (0, 0, 0, 128, 136, 40, 77, 64))); var Result : extended; begin printAll; writeln; Result := 6400 * arccos( sin(PI_180*one.Latitude.a)*sin(PI_180*two.Latitude.a) + cos(PI_180*one.Latitude.a)*cos(PI_180*two.Latitude.a)* cos(PI_180*(one.Longtitude.a-two.Longtitude.a)) ); one := one_; two := two_; printAll; writeln; try Result := 6400 * arccos( sin(PI_180*one.Latitude.a)*sin(PI_180*two.Latitude.a) + cos(PI_180*one.Latitude.a)*cos(PI_180*two.Latitude.a)* cos(PI_180*(one.Longtitude.a-two.Longtitude.a)) ); except writeln("EXCEPTION!"); readln; end; readln; end.

( one_ и two_ - это побайтовый "слепок" тех самых проблемных значений)

Т.е. наблюдаем именно то, чему было посвящено мое сообщение под номером 10: FloatToStr() результат возвращает одинаковый, но реально числа там разные, причем эта разница в данном случае оказывается критичной.

По простой аналогии, мне кажется, что если в Вашем примере Вы младшую 4 руками замените на 0, то эффект будет тот же.

Числа с плавающей точкой Найти похожие ветки