Я попытаюсь привести более-менее реальный пример законного, значимого и полезного тестирования на равенство с плавающей точкой.
#include <stdio.h>
#include <math.h>
/* let's try to numerically solve a simple equation F(x)=0 */
double F(double x) {
return 2*cos(x) - pow(1.2, x);
}
/* I'll use a well-known, simple&slow but extremely smart method to do this */
double bisection(double range_start, double range_end) {
double a = range_start;
double d = range_end - range_start;
int counter = 0;
while(a != a+d) // <-- WHOA!!
{
d /= 2.0;
if(F(a)*F(a+d) > 0) /* test for same sign */
a = a+d;
++counter;
}
printf("%d iterations done\n", counter);
return a;
}
int main() {
/* we must be sure that the root can be found in [0.0, 2.0] */
printf("F(0.0)=%.17f, F(2.0)=%.17f\n", F(0.0), F(2.0));
double x = bisection(0.0, 2.0);
printf("the root is near %.17f, F(%.17f)=%.17f\n", x, x, F(x));
}
Я бы не стал объяснять метод деления пополам Использовал сам, но акцентировал внимание на условии остановки.Он имеет точно обсуждаемую форму: (a == a+d), где обе стороны являются плавающими: a - это наше текущее приближение корня уравнения, а d - наша текущая точность.Учитывая предварительное условие алгоритма - что должен быть корнем между range_start и range_end - мы гарантируем на каждой итерации, что корень остается между a и a+d в то время как d уменьшается вдвое на каждом шаге, сужая границы.
И затем, после ряда итераций, d становится настолько маленьким , что при добавлении с a онокругляется до нуля!То есть a+d оказывается ближе к a, а затем к любому другому значению с плавающей запятой ;и поэтому FPU округляет его до ближайшего значения: до a.Это может быть легко проиллюстрировано расчетом на гипотетической вычислительной машине;пусть у него будет 4-значная десятичная мантисса и большой диапазон экспонент.Тогда какой результат машина должна дать 2.131e+02 + 7.000e-3?Точный ответ - 213.107, но наша машина не может представить такое число;это должно округлить это.И 213.107 намного ближе к 213.1, чем к 213.2 - поэтому округленный результат становится 2.131e+02 - небольшое слагаемое исчезло, округленное до нуля.Точно то же самое гарантированно произойдет на некоторой итерации нашего алгоритма - и в этот момент мы больше не можем продолжать.Мы нашли корень с максимально возможной точностью.
Поучительный вывод, по-видимому, состоит в том, что поплавки хитры.Они настолько похожи на реальные числа, что каждый программист испытывает желание думать о них как о реальных числах.Но это не так.У них свое поведение, немного напоминающее реальных , но не совсем то же самое.Вы должны быть очень осторожны с ними, особенно при сравнении на равенство.
Обновление
Пересматривая ответ через некоторое время, я также заметил интересный факт: в приведенном выше алгоритмеодин не может фактически использовать "некоторое небольшое число" в состоянии остановки.Для любого выбора номера будут входы, которые будут считать ваш выбор слишком большим , вызывающим потерю точности, и будут входы, которые будут считать ваш выбор тожемаленький , вызывая лишние итерации или даже вход в бесконечный цикл.Далее следует подробное обсуждение.
Возможно, вы уже знаете, что в исчислении нет понятия «маленькое число»: для любого действительного числа вы можете легко найти бесконечно много даже меньших.Проблема в том, что одним из этих «еще меньших» может быть то, что мы на самом деле ищем;это может быть корнем нашего уравнения.Хуже того, для разных уравнений могут быть отдельные корни (например, 2.51e-8 и 1.38e-8), оба из которых будут аппроксимированы одним и тем же числомесли бы наше условие остановки было бы похоже на d < 1e-6.Какой бы «маленький номер» вы ни выбрали, многие корни, которые были бы правильно найдены с максимальной точностью при условии остановки a == a+d, будут испорчены из-за того, что «эпсилон» будет слишком большим .
Это правда, однако, что в числах с плавающей запятой показатель степени имеет ограниченный диапазон, так что вы можете найти наименьшее ненулевое положительное число FP (например, 1e-45 denorm для FP с одинарной точностью IEEE 754),Но это бесполезно!while (d < 1e-45) {...} будет зацикливаться вечно, предполагая одинарную точность (положительное ненулевое значение) d.
SettiПомимо этих патологических краевых случаев, любой выбор "малого числа" в условии остановки d < eps будет слишком маленьким для многих уравнений. В тех уравнениях, где корень имеет достаточно высокий показатель степени, результат вычитания двух мантисс, отличающихся только наименьшей значащей цифрой, легко превысит наш «эпсилон». Например, с 6-значными мантиссами 7.00023e+8 - 7.00022e+8 = 0.00001e+8 = 1.00000e+3 = 1000, что означает, что наименьшая возможная разница между числами с показателем +8 и 5-значными мантиссами составляет ... 1000! Который никогда не будет вписываться, скажем, 1e-4. Для этих чисел с (относительно) высоким показателем у нас просто недостаточно точности, чтобы когда-либо увидеть разницу в 1e-4.
Моя реализация выше учитывала и эту последнюю проблему, и вы можете видеть, что d уменьшается вдвое за каждый шаг, а не пересчитывается как разница (возможно, огромная по показателю) a и b. Поэтому, если мы изменим условие остановки на d < eps, алгоритм не будет застревать в бесконечном цикле с огромными корнями (это очень хорошо может быть с (b-a) < eps), но все равно будет выполнять ненужные итерации во время сжатия d ниже точность a.
Такое рассуждение может показаться чрезмерно теоретическим и излишне глубоким, но его цель - снова проиллюстрировать хитрость поплавков. При написании арифметических операторов вокруг них нужно быть очень осторожным с их конечной точностью.