Округление дробных величин при сложении и вычитании

Проблема в том, что эти ошибки — в районе 1e-309 на каждое сложение-вычитание, а я работал с величинами 1e-10 — 1e-40. Для человека разница очевидна, но как её задать для машины?

Бери пакет с символьными вычислениями, а не численными. Проблема зарыта в представлении чисел с плавающей точкой.

Так-то да, но вопрос о том, как делают в простых калькуляторах. Или есть легковесные символьные калькуляторы?

bc -l. Вроде символьный. Я последнее время только им и пользуюсь

Не символьный, а с произвольной точностью, насколько я знаю. Грубо говоря, не загоняет все числа в 8-байтное double, а считает в столбик, рисуя столько знакомест, сколько нужно.

Надоело отсутствие в нём ряда функций, которые во wcalc «просто есть» и удобно реализованы.

Не символьный, а с произвольной точностью, насколько я знаю.

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

В одной статье как раз на эту особенность bc обратили внимание (в ответ на вопрос читателя): «bc считает с заранее заданной точностью, но не гарантирует, что результат содержит желаемое кол-во корректных цифр». И посоветовали альтернативу: spigot. В нём сколько заказали, столько получили, не ломая голову, какую же надо задать точность, чтобы на выходе получить столько-то значащих цифр.

Вот, например, где bc срабатывает нежелательно, а wcalc и spigot, как хотелось бы:

$ echo "scale=10; s(0.3)*3" | bc -l
.8865606198

$ wcalc -r -P10 "sin(0.3)*3"
~= 0.8865606199

$ spigot -d 10 "sin(0.3)*3"
0.8865606199

Теперь интересно, может в spigot ещё больше функций не хватает, чем в wcalc?

Кстати, было бы интересно увидеть минимальный пример, иллюстрирующий проблему (в wcalc).

Кстати, было бы интересно увидеть минимальный пример, иллюстрирующий проблему (в wcalc).

-> 2.2205e-16
 = 2.2205e-16
-> 2.2204e-16
~= 0
-> \p 1
Precision = 1
 = 2.2e-16
-> \p -1
Precision = auto
~= 0
-> \cons
Now Not Using Conservative Precision

-> 2.2204e-16
~= 2.220399999e-16
->

То есть, если результат вычислений меньше определённого порога, включён режим «accuracy guards» (команда \cons), и точность выбирается автоматически (команда \p -1), то выводится результат округлённый до нуля. В данном случае операций никаких, просто вывод введённого числа, но для любых вычислений одинаково.

-> \cons
Now Not Using Conservative Precision

Т.е. этот режим можно отключить?

~e-16 - это же epsilon для double. Значит, wcalc - вообще без произвольной точности? И в man там нет слова arbitrary.

Насчёт примера: я имел ввиду минимальные вычисления, которые приводят к такому-то результату (а ожидался вот такой), иначе мне не совсем понятно, что нужно.

И пока нет контраргумента против spigot.

Теперь интересно, может в spigot ещё больше функций не хватает, чем в wcalc?

Сравнил http://www.chiark.greenend.org.uk/~sgtatham/spigot/spigot.html и http://w-calc.sourceforge.net/manual.html Ебилда spigot нет, поэтому пока не ставил. Некоторые функции не упомянуты в документации wcalc, например exp, я их пытался искать по автодополнению.

В spigot есть, во wcalc нет:

• Несколько способов вычислять остаток в одном выражении. wcalc требует включать-выключать режим «C-style mod».
• hypot — гипотенуза по катетам.
• atan2 — арктангенс для заданного вектора (x; y).
• Отдельные наборы тригонометрических функций для градусов и радианов, wcalc требует переключать режим.
• Отдельные функции для возведения в степень 2 и 10. Логарифмы по основаниям 2 и 10 есть в обоих.
• pow — степень как функция от 2 аргументов, как в Си. exp, ^ и ** есть в обоих.
• frac — дробная часть. ceiling и floor есть в обоих.
• expm1, log1p — exp(x) - 1 и log(1 + x).
• erf, erfc, Phi, norm — функция ошибок, она же -1, плотность нормального распределения, она же.
• erfinv, inverf, erfcinv, inverfc, Phiinv, norminv, invPhi, invnorm — обратные к ним.
• W и Wn — ветви W-функции Ламберта.
• Ei, En, E1, Ein — 1-е и n-е интегралы exp(x)/x с разными константами.
• li, Li(x) — интегралы 1/log(x) с разными константами.
• Si, si, Ci, Cin — интегралы sin(x)/x с разными константами, cos(x)/x и (1-cos(x))/x.
• FresnelS, FresnelC, UFresnelS, UFresnelC — интегралы Френеля: sin(π x^2/2) и cos(π x^2/2), нормализованные и ненормализованные.
• algebraic — найти корень многочлена на заданном интервале. Если верить мануалу, глючит, если корень не один, или если параметры не рациональные.
• phi — постоянная золотого сечения, (1+sqrt(5))/2.
• apery — постоянная Апери. Дзета-функция есть в обоих.
• Любое основание системы счисления, но не выше 36.

Есть во wcalc, нет в spigot:

• Переключение режимов градусы-радианы, приходится заменять все функции :)
• round — округление. ceiling и floor есть в обоих.
• rand, irand — случайные числа.
• sinc — кардинальный синус.
• K — постоянная Каталана.
• Порядка 40 физических констант. g, G, k, F, а.е.м., масса дейтрона...
• Возможность записывать эти константы греческими буквами в UTF.
• Сохранение переменных в файл конфигурации. spigot, если я правильно понял, умеет писать-читать только 1 число на файл.

Реально из перечисленного мне бывали нужны только функция ошибок и константы. algebraic и hypot могут пригодиться. Сохранение переменных тоже. Вывод: стоит попробовать поставить.

Т.е. этот режим можно отключить?

Да, проблема только в том, что он по умолчанию. Поэтому я спросил, как её обычно решают в калькуляторах.

~e-16 - это же epsilon для double.

Как это отражается на вычислениях, где все аргументы ~e-20? А если идёт умножение 1e3 на 1e-23?

Значит, wcalc - вообще без произвольной точности? И в man там нет слова arbitrary.

Похоже. И я не говорил, что он произвольной точности.

И пока нет контраргумента против spigot.

Долго писал сравнение :) И нет ебилда. Не найду — напишу.

~e-16 - это же epsilon для double.

Как это отражается на вычислениях, где все аргументы ~e-20? А если идёт умножение 1e3 на 1e-23?

Есть правило: при сравнение двух чисел с плавающей запятой, если их абсолютная разница <= машинного эпсилон, то их уже можно считать равными. Если я не ошибаюсь, можно сравнить с динамическим диапазоном: в одной операции имеет смысл оперировать числами в его пределах, не более. Разумеется, результат операции округлять до нуля в таком случае не нужно. А вот при использовании двух чисел с бОльшим динамическим диапазоном, меньшее всё равно будет что 0. И в таком случае калькулятор на double не подходит, если порядок операций нельзя поменять так, чтобы исключить операции с сильно различающимися по размеру числами.

Долго писал сравнение :) И нет ебилда. Не найду — напишу.

Спасибо за развёрнутое сравнение! В зависимостях у spigot из особо важного только GMP. Автор отзывчивый, если что-то интересное прояснишь, пиши.

Проблема в том, что эти ошибки — в районе 1e-309 на каждое сложение-вычитание, а я работал с величинами 1e-10 — 1e-40.

Тут написано что-то странное. Хотя 1e-309 приблизительно равно минимальному числу, которое представимо в виде double, это не значит, что погрешность будет всегда такой маленькой. В лучшем случае, она будет порядка 1e-17 от тех чисел, с которыми ты работаешь. То есть для 1e-10 это будет порядка 1e-27.

Чтобы найти эту погрешность, попробуй прибавлять небольшие числа к одному выбранному, уменьшая эти числа, пока результат не перестанет меняться:

>>> a = 1e-10
>>> b = a
>>> while a + b != a:
...   b = b * 0.999
... 
>>> b
6.457899767711049e-27
>>> a
1e-10
>>> b/a
6.457899767711049e-17
>>> 

Проблема в том, что эти ошибки — в районе 1e-309 на каждое сложение-вычитание, а я работал с величинами 1e-10 — 1e-40.

Тут написано что-то странное

Оговорился. Уже не помню, как я хотел это сформулировать.

Ошибки округления — последний бит мантиссы, То есть для 8-байтного с плавающей запятой 2^(-53)~1e-16 = 0.00000000000001% от входных данных. А в режиме «консервативной» точности округляет до нуля входные данные, которые по абсолютной величине меньше 2e-16. Для человека может быть очевидно, какой результат ожидать, но калькулятор округляет, как если бы всегда ожидал получить величины близкие к 1.

Хм, порог, при котором округляет до нуля, оказался очень близок к 2^(-52).

Используют рациональные числа, где только возможно → https://godoc.org/robpike.io/ivy