[C++] реализация чисел с плавающей точкой на int'ах

GMP - то, что доктор прописал.

Он таки написан на C и таки gnu'сным образом. Как у него с переносимостью в винды?

> Он таки написан на C

http://www.gnu.org/software/gmp/manual/html_node/C---Class-Interface.html#C++...

>Он таки написан на C

Там есть С++-интерфейс, как другой анон показал

таки gnu'сным образом

Вы говорите так, будто это что-то плохое

Как у него с переносимостью в винды?

Учитывая, что он используется в gcc (емнип), который в свою очередь прекрасно пашет под виндой, думаю, что нормально.

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

строки из интов, сюрприз?

Насколько помню, в ПЛ/1 было что-то подобное :)

Правда? И как часто происходит это «обычно»? Интересует также скорость подобных вычислений по сравнению с традиционными численными.

Интересно, что это за ограничения?

А ещё, может возможностей long double хватит? Под MSVC, правда, не работает, зато Intelовский компилер отлично с ним справляется.

А если нужна бОльшая точность - GMP подойдет.

Объясните кто-нибудь, в чем разница между GMP и MPFR?

Ну и в тему больших длинных типов, есть ещё __float128 в gcc и аналогичный _Quad в icc. С ними могут быть некоторые проблемы (математика стандартная не факт, что захочет с ними работать), но это не так и страшно. А MSVC снова в пролете.

В первую очередь имеется в виду вычитание двух близких по значению double'ов. Ну и численное интегрирование с маленьким шагом тоже меня несколько удивило. А long double отсутствует в стандарте или MSVC? Потому как для винды я использую mingw, который его разумеется умеет.

long double есть в стандарте C99, и поддерживается в gcc и icc. Возможно, он есть в новом стандарте C++0x, слету не скажу, надо смотреть.

В Микрософте же решили выделиться, поэтому в MSVC long double - синоним double.

В микрософте не решили выделяться, а просто забили на С99, в чем честно признаются.

Если long double не подойдет, то можно попробовать __float128 (gcc) / _Quad (icc). Причем для _Quad интеловский компилятор умеет всю стандартную математику, что несомненно плюс.

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

В микрософте не решили выделяться, а просто забили на С99, в чем честно признаются.

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

Впрочем, на многих задачах double действительно достаточно.

Алгоритмы должны быть максимально human readable, чтобы их мог проверять на семантическую корректность человек далёкий от double'ов.

Это все хорошо, конечно. Но давайте я два примера приведу.

• Надо два double на равенство сравнить, оба получены в результате каких-то вычислений. Способ a == b конечно хорош, и отлично читаем человеком, далеким от этих сложных вычислений. Только не работает. Недавно тут один товарищ на эту тему жаловался.
• Надо было как-то раз получить номер ячейки (все одинакового размера и начинаются с 0), в которой лежит заданная точка. Метод

  int num = (int)floor(x / size_x);
  

  или даже

  int num = (int)(floor(x / size_x) + 0.5);
  

  очень прост для восприятия и понимания. Только вот не работает, из-за чего возникают странные баги. К примеру, число x=3 в памяти в результате вычислений представлялось как 2,99...94, и при вычислении (размер ячейки - единица) радостно превращалось в 2, что слегка ожиданиям не соответствует.

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

>Поэтому собственно вариантов несколько: либо можно обратиться к умным математикам...

Лучше не обращаться. Будет только хуже. Умные математики далеки от земных проблем.

>Как у него с переносимостью в винды?

Гугл может тебе много рассказать как под виндой GMP использовать.

Вот одна из ссылок.

http://cs.nyu.edu/exact/core/gmp/

C99 не требует условие long double > double

Я когда то для своих целей делал через int_32 и int_64:
/**Currency representation*/
typedef struct
{
int64_t data;
decimals_t decimals;
} currency_t;

/**Float number representation*/
typedef struct
{
int32_t data;
decimals_t decimals;
} float_t;

и соотв. функции арифметические, сравнения, округления и т.п. с контролем переполнения.
В течении 14 дней можно взять тут http://www.sendspace.com/file/230rtt

> К примеру, число x=3 в памяти в результате вычислений представлялось как 2,99...94, и при вычислении (размер ячейки - единица) радостно превращалось в 2, что слегка ожиданиям не соответствует. floor - отбрасывает целую часть, т.е. floor(2,999) = 2, Скобку в другое место ставь: int num = (int)(floor(x / size_x + 0.5); - получишь floor(2,999 + 0,5) = 3, т.е. то что ты ожидал

ещё есть фиксированная запятая http://ru.wikipedia.org/wiki/Фиксированная_запятая http://en.wikipedia.org/wiki/Fixed-point_arithmetic но у неё есть проблемы с погрешностью, и для интегрирования точности скорее всего не хватит

Скобку в другое место ставь: int num = (int)(floor(x / size_x + 0.5); - получишь floor(2,999 + 0,5) = 3, т.е. то что ты ожидал

Ну и сразу вылезет floor(2,6 + 0,5) = 3, т.е. 2,6 вместо двух превратится в 3, что тоже довольно странно.

Помню одну книгу по выпуклому анализу, так там эпсилон-окрестности чуть ли не на каждой странице... В одном графическом редакторе я использую floor для определения участка перерисовки, и ничего, работает :)

что странного? математическое округление же. 2,6 округляется до 3.

p.s. особых проблем с double не имел, если четко сформулирована задача.

Ну в графическом редакторе все и должно работать. :)

А у был очень баг из-за этого флоата: при создании больших графов (~10^7 ребер) в некоторых случаях возникала описанная ошибка, что приводило к тому, что группа ребер из графа ошибочно выбрасывалась, после чего на этом графе вычисления были неправильными. А на маленьких графах ничего не возникало. Был в восторге от процесса отладки.

Кстати, по поводу сравнения double'ов - вполне устроил вот этот пост, может кому пригодится.

что странного?

Прочитайте исходный пост. Речь об округлении не шла. И формулировка там предельно четкая (п. 2).

Когда-то давно держал в руках учебник по интервальной математике, так там утверждалось, что довольно приличное количество численных методов при доказанной сходимости в реальности могут даже не сходится из-за неточных вычислений, т.е. погрешности! Другими словами, теоретическая сходимость без учета интервалов не дает еще гарантии. Некорректные методы получаются. Вот такие дела. Вообще, сложная тема. Очень матанная.

ок, просто надо понимать что пишете, пост читал бегло, floor(x + .5) - это стандартное математическое округление, оно вам не к чему. в вашем случае нужно добавлять небольшое число, к примеру floor(x + .00001)

p.s. пишу опять же бегло, но смысл я думаю ясен - добавляем небольшое значение, которое подходит под условия задачи

Да-да, согласен.

Кстати, а можно ссылку на учебник? :) Или хотя бы название и автора?

Прошло уже одинадцать лет. Но название примерно такое: «Интервальная математика» или «Интервальные вычисления». Переплет твердый, хороший. Коричневого цвета. Кажется, авторы иностранные. Уже не помню.

А так, во многих случаях спасает регуляризация или прочие хитрые методы. И как говорил один мой приятель, который занимался краевыми задачами в аспирантуре: «ни за что не сяду в тот самолет, который будет спроектирован с учетом моих результатов!» ;)

хотел ведь написать «одиннадцать» :)

ок, просто надо понимать что пишете, пост читал бегло, floor(x + .5) - это стандартное математическое округление

Вы, пожалуйста, определитесь, с кем хотите спорить: со мной или с анонимусом (который тоже бегло прочитал, не разобрался и что-то не то сказал). И заодно определитесь, что именно Вы хотите оспорить.

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

нужно добавлять небольшое число, к примеру floor(x + .00001)

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

Ок, спасибо, постараюсь найти.

Ага, а самолет вот назло возьмет и взлетит. :)

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

lol. мсье теоретик? себе руки отрывайте, когда свои задачи решить не можете.

во-первых - я вам предложил не DBL_EPSILON а конкретную величину, во-вторых я указал

добавляем небольшое значение, которое подходит под условия задачи

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

как правило точности хватает. программы как-то пишутся, задачи как-то решаются... мистика...

если нужна сверх-точность - то встают вопросы как автора темы.

во-первых - я вам предложил не DBL_EPSILON а конкретную величину

Вы в курсе, чем отличается «относительный» от «абсолютный»?

Прочитайте ещё раз про слона в посудной лавке. Ваша 0.0001 мало чем от DBL_EPSILON в этом случае не отличается.

не отличается = отличается же.

спорьте сами с собой дальше.

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

как правило точности хватает. программы как-то пишутся, задачи как-то решаются... мистика...

> Ага, а самолет вот назло возьмет и взлетит. :)

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

А у меня все проще. Упала или подвисла программа - клиент ее просто перезапустил :)

> если нужна сверх-точность - то встают вопросы как автора темы.

Это еще надо смотреть, что за задача. А вдруг она некорректная?

согласен )

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

Не поможет. Ну будет тебе польза от ответа 0.1±1e+10? Я как раз сейчас такой алгоритм реализовывал, там одна из стадий — суммирование бинома типа (1-x)^n при х близком к 1. Представляешь, какая там скорость накопления ошибок? Помогает только арифметика высокой точности (реализация была нужна для проверки точности численного многомерного интегрирования, которое оказалось быстрее и точнее аналитического как раз из-за неустойчивости аналитических формул)

И какую реализацию арифметики высокой точности использовали, если не секрет?

Использовал Maple, но можно было легко реализовать на GMP. В общем-то, было три уровня: рекуррентное соотношение с начальным условием вроде логарифма (устойчивая рекурсия), потом сбор значений промежуточных интегралов (суммировал два значения рекурсивной функции), потом подстановка полученных значений в двойное суммирование с биномиальными коэффициентами. На последней стадии неустойчивость была такая, что 50-й член последовательности, будучи суммируем с 50-значной же точностью, давал всего 10 правильных цифр ответа... Оценил асимптотику — число верных знаков уменьшалось со скоростью 4^n, где n — номер искомого интеграла, из-за потери точности на вычитаниях.

Пардон, наврал: на 50-значной точности было только 25 знаков правильного ответа для 50-ого члена последовательности

«стандартная» же арифметика на 15-значной точности давала ответ на 7 порядков больший правильного, т.е. реализовался как раз случай 0.1±10⁷. Понятно, что мантисса НИКАКОГО отношения к правильному ответу даже не имела. Первая правильная значащая цифра появилась при 25-значной арифметике

Если интересно, потом на почту могу послать ссылку на соотв. заметку в блоге (только напишу ее) :)

а чего на почту-то? сюда кидай