Как поднять точность вычислений

С Новым Годом всех!

Проблема вот какая.
Процессоры Intel, и не только они, имеют паразитное свойство снижать общую точность. Например:

a[color=#808030; ]=[/color][color=#008c00; ]0.0[/color] do i[color=#808030; ]=[/color][color=#008c00; ]1,100[/color] a[color=#808030; ]=[/color]a[color=#808030; ]+[/color][color=#008c00; ]0.00001[/color] end do

Если на каждом шаге выводить значения "а", то уже через некоторое время получается что-то вроде 0.0001101. Более наглядно это видно если "а" имеет двойную точность. Так сказать ошибки представления чисел с плавающей запятой в результате работы сопроцессора. Как обычно от этого избавляются в настоящее время? Совершенно очевидно только то, что даже крупные коммерческие производители программного обеспечения этому внимания не уделяют. Как пример - реализация статистических функций в Excel или известный статистический пакет SPSS, результатам расчетов которых доверять нельзя. (Единственный статистический пакет, который работает корректно - это Statistica). Аналогично хамят и многие расчетные программы - например Elcat - отечественная конечно-элементная система моделирования, когда дело касается решения нестационарных задач при малых градиентах расчетных величин. ("Монстра" ANSYS я еще не проверял). Похоже надо как то производить коррекцию на ассемблере, но не хочется изобретать велосипед.
(Обычные методы увеличения точности, как то перестановка переменных, оптимизация алгоритма не шибко-то помогают).

Только что посчитал ваш пример в матлабе. Ответ: 0.00100000000000
Процессоры интел непричём. Представление чисел с плавающей точкой не позволяет представить все возможные числа. Простых выхода для вас два: считайте или в максимально длинном формате с плавающей точкой чтобы уменьшить ошибки округления, или в целых числах, хоть это и неудобно. Или к нам, в президиум Т.е. матлабом. Там математика написана грамотно. Можно написать функцию в матлабе и скомпилировать в ехе, или прямо из визуала использовать.

au

Возьмите любой компилятор фортрана от Microsoft, DEC или Compaq, любой компилятор Basic от PDS 7.1 до Visual Basic 6.0, на других не проверял, запретите использование эмулятора и посмотрите что получится.
С матлабом я не очень знаком. Кстати, там библиотеки в каком виде?

Это в догонку.
У меня все тестовые примеры на работе, так что 5-го выложу код и результат.
Кстати, вот еще какое дело. Некоторые мои расчетные программы были сделаны еще под ДВК, потом просто переписаны для х86. Так вот: один и тот же алгоритм дает разную достоверность результатов для ДВК и х86, причем для х86 хуже несмотря на вроде большую точность вычислений. Есть большое подозрение на механизм работы сопроцессора в х86, когда число с плавающей точкой, независимо от точности представления, преобразуется в 80-бит формат представления сопроцессора, а потом обратно.
Парадокс, но самые точные и достоверные результаты получались у меня (при простейшем расчете одномерных температурных полей, нестационарная задача) под старо-старой Искра-1256, в которой вся плавающая арифметика была реализована микропрограммно а не аппаратно.

1. Библиотеки у Матлаба в основном в исходниках (*.М - файлах), что очень удобно

2. Насколько я знаю, у Матлаба все числа по умолчанию хранятся как раз в формате сопроцессора (80-бит), так что лишних преобразований не производится.

3. Возможно, при переносе программ со старых машин необходимо явно объявлять переменные как 80-битные ("long double" в С, "extended" в Паскакале), иначе используется "старый" 32-64-битный формат.

Кстати, у меня в Матлабе пример (первый самый) выдал
A=0.001000000000000002

Точнее,
A = 1.0000000000000020000e-003

Э-э-э. немного неожиданно, прямо скажем. Неужели не у кого не было курса приближенных вычислений в ВУЗе?

Очень точные вычисления тут лежат. Вообще-то ни MS, ни Борланд ни другие "крупные" ребята никогда особенно не заботились о точности. Заботились те, для кого это важно было (смотри ссылку выше).

По памяти - книга Крылова вроде была не плоха. Есть разные специальные методы для повышения точности.

А еше гляньте сюдао.

[quote|GrayCat, 02.01.2004 00:14:47:] 2. Насколько я знаю, у Матлаба все числа по умолчанию хранятся как раз в формате сопроцессора (80-бит), так что лишних преобразований не производится.


[/QUOTE]
Весь сопроцессор реализован в соответствии со стандартом IEEE.

au, 01.01.2004 15:43:47:

Только что посчитал ваш пример в матлабе. Ответ: 0.00100000000000
Процессоры интел непричём. Представление чисел с плавающей точкой не позволяет представить все возможные числа. Простых выхода для вас два: считайте или в максимально длинном формате с плавающей точкой чтобы уменьшить ошибки округления, или в целых числах, хоть это и неудобно. Или к нам, в президиум Т.е. матлабом. Там математика написана грамотно. Можно написать функцию в матлабе и скомпилировать в ехе, или прямо из визуала использовать.

 

В Mathlabe та же проблема появиться если оперировать с числами порядка 1е-13

Ну, это не проблема Малтаба etc. - это проблема сопроцессора, и, следовательно, всех программ, на него полагающихся. Если для какой-то задачи требуется бОльшая точность, или если алгоритм предполагает такие "плохие" операции, как в примере, следует пользоваться программной реализацией математики, с юбой нужной точностью.

В данном же примере закавыка состоит в том, что внутреннее двоичное представление чисел не может точно отразить десятичное их представление. Например, 1e-5 сопр "видит" как 1.0000000000000001000e-005, ну и дальше, соответственно, накапливает ошибку.

Если же взять число, хорошо "ложащееся" на двоичное представление, например, 1/(2¹⁶)==1.525878906250e-005, то никаких ошибок нет вообще. Т.е. суммирование 100 раз даст ровно 1.5258789062500000000e-003.

Так что, "Наличие компьютера не избавляет от необходимости думать" © мой

Я понял. С температурными полями и прочим таким — в матлаб однозначно.

Матлаб все хранит в 64-х битном формате. Т.е. даже 80-бит нема. Вот.
Машинный нуль Матлаба ~ 2.2e-16.

Плавающая арифметика - полнейший отствой. Настоящие программисты её никогда не используют.

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

[quote|Mishka от 02.01.2004 02:56:12:]Э-э-э. немного неожиданно, прямо скажем. Неужели не у кого не было курса приближенных вычислений в ВУЗе?[/QUOTE]

Лично у меня не было. Обычно такие курсы идут у тех, кто обучается по специальности "Прикладная математика". В начале 1980-х в МИСиС такого не было. В МВМИ, который я кончал в 1986 уж подавно. Все изучал самостоятельно.

[quote|Mishka от 02.01.2004 02:56:12:]Очень точные вычисления тут лежат. Вообще-то ни MS, ни Борланд ни другие "крупные" ребята никогда особенно не заботились о точности. Заботились те, для кого это важно было (смотри ссылку выше).[/QUOTE]

NAG мне известен, поскольку их библиотеки входят по умолчанию в поставку современных компиляторов фортрана. Но их исходники - чисто классические.

[quote|au от 02.01.2004 15:20:34:]Плавающая точка на то и плавающая, чтобы можно было работать и с очень маленькими числами, и с очень большими, но не смешивая. Ваша задача требует этого?[/QUOTE]

Естественно. При расчете температурных полей конечно-разностным методом с одной стороны числа порядка 10³, с другой порядка 10^-6. В регрессионном анализе без них вообще никуда не денешься.

Может быть в трехтомнике "Искусстве программирования" что-нибудь есть? Сам я эту книгу не читал.

AidarM, 02.01.2004 13:21:07:

Матлаб все хранит в 64-х битном формате. Т.е. даже 80-бит нема. Вот.
Машинный нуль Матлаба ~ 2.2e-16.

 

Тьфу ты, черт, так и есть, 8 байт... Бли-и-ин, аЦтой...

>Тьфу ты, черт, так и есть, 8 байт... Бли-и-ин, аЦтой...

Вот именно. Пишите вставочки на Borland С. Visual С - такой же отстой. Точь-в-точь. Лично я стараюсь в целые числа переводить.

Ну давайте не отвлекаться на частности, типа VC, Borland C, Matlab и т.п. Любые вычисления с плавающей точкой в принципе обладают лишь ограниченной точностью. Тут правила простые: чем меньше операций, тем лучше; группируйте операции с числами близкими по абсолютному значению. Это простые и интуитивные правила, которые многие забывают.

Вот вы, Centuriones, привели пример , который является хрестоматийным примером того как НЕЛЬЗЯ работать с таким типом данных. Чем вам не понравилось такое решение (если вам надо получит список чисел с интервалом в 0.00001):

do i=1,100
a=(i-1)*0.00001
end do

В вашей версии при длине цикла порадка O(N) ошибка имеет порядок O(N), а в этом примере O(1) (то есть не зависит от длины цикла). Сравните эти два метода на одной и той же машине и вы увидите разницу. И в ответ на ваш вопрос - от ошибок в точности вичислений избавиться невозможно - можно их только снизить. И именно с помощью грамотной организации алгоритмов. Этим занимается целый раздел прикладной математики и без этого было бы невозможным решение многих рутинных числовых задач с приемлемой точностью, неважно на каком железе.

GrayCat, 02.01.2004 20:50:52:

AidarM, 02.01.2004 13:21:07 :

Матлаб все хранит в 64-х битном формате.  Т.е. даже 80-бит нема. Вот. 
Машинный нуль Матлаба ~ 2.2e-16.

 

Тьфу ты, черт, так и есть, 8 байт... Бли-и-ин, аЦтой...

 

Тут дело такое - внутренне сопроцессор всегда работает с 80 битами. Это называется временный формат. Его специально используют для того, чтобы точность у чисел с 64 битами повысить. 32 битные внутри тоже преобразуются в 80 битные. Если записать изначальный алгоритм на ассемблере, то ошибок не будет.

Rada, 02.01.2004 23:16:19:

Ну давайте не отвлекаться на частности, типа VC, Borland C, Matlab и т.п. Любые вычисления с плавающей точкой в принципе обладают лишь ограниченной точностью. Тут правила простые: чем меньше операций, тем лучше; группируйте операции с числами близкими по абсолютному значению. Это простые и интуитивные правила, которые многие забывают.

 

Ага, я про тоже. Могу только добавить, что особая аккуратность нужна с операциями сложежния/вычитания, умножение/деление намного более устойчивы. Отсюда вывод, что для очень маленьких и очень больших заводят разные аккумуляторы, которые тоже могут объединять (суммировать) позже. Не стоит забывать, что все числа храняться в нормализованном виде. Для * и / это нормально, а вот для + и - приходится числа денормализовывать (если у них порядок разный). Не стоит также забывать, что из 64 разрядов, 12 отводяться на характеристику (степень), 1 на знак и 52 на мантиссу (базу). Да-да 52. Т.к. числа храняться в нормлизованном виде, то одна единичка предпологается по умолчанию и в сетке не храниться. Таким образом, на представление числа отводиться 53 бита - последние 2или3 как правило не считаются надежными - в них очень быстро накапливаются ошибки. Отсюда вывод - если вы попытаетесь сложить два числа с разнице в порядке превышающим 2⁵³, то у вас большее число не измениться. 2⁵³ это примерно 10^16. Если брать 32 битное число, то там 8 бит степень, 1 знак и 23(24) мантисса. Т.е. для точности всего 23(24)бита. Или, если разница между величинами достигает 10⁷, то мы уже наблюдаем проблемы.

Вот, нашел, стандарт IEEE 754, вот тут можно посмотреть немного подробнее.

Вот вы, Centuriones, привели пример , который является хрестоматийным примером того как НЕЛЬЗЯ работать с таким типом данных. Чем вам не понравилось такое решение (если вам надо получит список чисел с интервалом в 0.00001):

do i=1,100
a=(i-1)*0.00001
end do

Ему здесь не понравилось то, что надо суммировать где-то результаты, простое суммирование дает эти странные ответы.

 

В вашей версии при длине цикла порадка O(N) ошибка имеет порядок O(N), а в этом примере O(1) (то есть не зависит от длины цикла). Сравните эти два метода на одной и той же машине и вы увидите разницу.

 

В качестве генерации констант Ваш метод очень не плох, но он не помогает в задаче автора

И в ответ на ваш вопрос - от ошибок в точности вичислений избавиться невозможно - можно их только снизить.

 

Золотые слова!

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

 

И развит он очень не плохо.

А вычисления с произвольной точностью - они не плохи, но очень медленны. Но иногда без них не обойтись.

Resurrector, 02.01.2004 15:06:19:

Плавающая арифметика - полнейший отствой. Настоящие программисты её никогда не используют.

 

А он не как программист ее использует, а как ученный. И перевод плавающей в целочисленную не помогает - остаются те же проблемы, если только не делать вычисления с произвольной точностью (но это отдельная песня).

GrayCat, 02.01.2004 04:47:12:

Ну, это не проблема Малтаба etc. - это проблема сопроцессора, и, следовательно, всех программ, на него полагающихся. Если для какой-то задачи требуется бОльшая точность, или если алгоритм предполагает такие "плохие" операции, как в примере, следует пользоваться программной реализацией математики, с юбой нужной точностью.

 

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

В данном же примере закавыка состоит в том, что внутреннее двоичное представление чисел не может точно отразить десятичное их представление. Например, 1e-5 сопр "видит" как 1.0000000000000001000e-005, ну и дальше, соответственно, накапливает ошибку.

 

Так эта проблема стоит для любого процессора с фиксированным представлением чисел.

Если же взять число, хорошо "ложащееся" на двоичное представление, например, 1/(2¹⁶)==1.525878906250e-005, то никаких ошибок нет вообще. Т.е. суммирование 100 раз даст ровно 1.5258789062500000000e-003.

 

А попробуй повторить операцию на простых (не двойных) плавающих 2¹⁸ раз и увидишь ошибку. Сумма перестанет возрастать.

Так что, "Наличие компьютера не избавляет от необходимости думать" © мой

 

О, это точно! Пускай комп считает, а программист думает!

[quote|Centuriones, 02.01.2004 18:01:29:] [quote|Mishka от 02.01.2004 02:56:12:]Э-э-э. немного неожиданно, прямо скажем. Неужели не у кого не было курса приближенных вычислений в ВУЗе?[/QUOTE]

Лично у меня не было. Обычно такие курсы идут у тех, кто обучается по специальности "Прикладная математика". В начале 1980-х в МИСиС такого не было. В МВМИ, который я кончал в 1986 уж подавно. Все изучал самостоятельно.


[/QUOTE]
Пардон, не знал. Я по старой привычке считал, что в Российском образовании это обязаловка.

Может быть в трехтомнике "Искусстве программирования" что-нибудь есть? Сам я эту книгу не читал.

 

Не припомню ничего такого особого. Приду домой посмотрю. У Кнута большей частью дискретная математика и анализ алгоритмов. Попробуй поискать Крылова. Мне казалось, что он прост. Есть статьи и монографии у Марчука, но они более суровы. Я бы поискал у физиков-математиков (те, которые применяют математику для рассчетов). Зайди в гугл и дай что-то вроде high precision calculation - я получил массу ссылок на реализацию конкретных алгоритмов с повышенной точностью - с них можно уже чему-то научиться. А вообще-то Maple, Mathematica - дают офигенную точность - но порой за счет скорости.

IEEE 754: "Ассоциативным можешь ты не быть, но коммутативным быть обязан":

(1.0e10 + -1.0e10) + 1.0e-10 not= 1.0e10 + (-1.0e10 + 1.0e-10)

но

1.0e10 + 1.0e-10 = 1.0e-10 + 1.0e10.

Кстати, может кто-нибудь подсказать: как реализован механизм плавающей арифметики на IBM 360/370 (т.е. ЕС ЭВМ) и на PDP-11 (LSI-11) (cоответственно наши "Электроника-60" и ДВК)?
Аппаратно или микропрограммно?