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

Question

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

real s;  //sum of all previous terms
ulong k; //factorial

s += 1.0/ k;

после каждого шага k становится еще меньше, но после 10-го раунда результат больше не меняется и застрял на 2.71828

Answer 1

Типы с плавающей запятой с фиксированной точностью, те, которые изначально поддерживаются модулем с плавающей запятой вашего ЦП (float, double, real), не оптимальны для любых вычислений, которые требуют много цифр точности, например, в примереВы дали.

Проблема в том, что эти типы с плавающей запятой имеют конечное число цифр точности (фактически двоичные цифры), что ограничивает длину числа, которое может быть представлено таким типом данных.Тип float имеет ограничение приблизительно в 7 десятичных цифр (например, 3.141593);тип double ограничен 14 (например, 3.1415926535898);и тип real имеет аналогичное ограничение (чуть больше, чем double).

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

float a = 1.234567f, b = 0.0000000001234567
float c = a + b;

writefln("a = %f b = %f c = %f", a, b, c);

Оба a и b являются действительными значениями с плавающей запятой и сохраняют приблизительно 7 цифр прецизионной единицы в отдельности.Но при добавлении сохраняются только первые 7 цифр, потому что они возвращаются обратно в число с плавающей точкой:

1.2345670001234567 => 1.234567|0001234567 => 1.234567
                              ^^^^^^^^^^^
                         sent to the bit bucket

Таким образом, c в конечном итоге равняется a, поскольку более точные цифры от добавленияa и b получают удар.

Вот еще одно объяснение концепции , вероятно, намного лучше, чем у меня.

---

Ответ на эту проблемуявляется арифметикой произвольной точности.К сожалению, поддержка арифметики произвольной точности не в аппаратном обеспечении процессора;следовательно, это не (обычно) в вашем языке программирования.Однако есть много библиотек, которые поддерживают типы с плавающей точкой произвольной точности и математические вычисления, которые вы хотите выполнить для них.См. этот вопрос для некоторых предложений.Вероятно, вы не найдете сегодня для этой цели каких-либо специфичных для D библиотек, но существует множество библиотек C (GMP, MPFR и т. Д.), Которые должны быть достаточно простыми для использования в изоляции, и даже более того, если вы сможете найтиD привязки для одного из них.

Answer 2

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

auto data = real[N];
foreach(i, ref v; data) {
  v = Fn(i);
}

while(data.length > 1) {
  data.sort(); // IIRC .sort is deprecated but I forget what replaced it.
  data[1] += data[0];
  data = data[1..$];
}

return data[0];

(минутная куча сделает это немного быстрее.)

Answer 3

Как уже упоминалось, вам нужно использовать стороннюю арифметическую библиотеку с плавающей точкой с множественной точностью (я думаю, что у Танго или Фобоса есть только модуль для целочисленной арифметики произвольной длины). - это проект D, использующий MPFR.Вы должны найти там привязки.