Обработка малых чисел

Эта публикация родилась из обсуждения одной задачки в блоге Команда ФПИ по программированию.
Итак, недавно, столкнулся с двумя похожими ситуациями, когда заведомо известный результат не получается при вычислении его программным образом.
Ситуация 1: Excel не дает ноль для тригонометрической функции cos(Pi/2).
Ситуация 2. FreePascal не дает ноль для функции frac(log2(128)) Команда ФПИ по программированию — функция log2.

Я и раньше встречался с подобными ситуациями. Особенно опасны они ненахождением решения в ситуации проверки равенства двух вещественных чисел. Разберемся, почему всё это происходит? Очевидно, что по причине ограниченности разрядной сетки числа. Например, посмотрим про вещественные типы данных в справке Delphi:
Type		Range	                        Significant digits
Real48		2.9x10^-39 .. 1.7x0^38			11-12
Single		1.5x10^-45 .. 3.4x10^38			7-8
Double		5.0x10^-324 .. 1.7x10^308		15-16
Extended	3.6x10^-4951 .. 1.1x10^4932		19-20

И сделаем два вывода:
1) нельзя задать сколь угодно малое или сколь угодно большое вещественное число, так как есть ограничение на порядок числа;
2) точность представления числа ограничена мантиссой (количеством значащих цифр).

Отсюда важное для программиста заключение: сравнивать вещественные числа так: «x>y» или так: «x<y» можно, но так: «x=y» — опасно.
Когда-то давно решал я подобную проблему, связанную с заполнением до предела (типа задачи «о заполнении рюкзака», ru.wikipedia.org/wiki/Задача_о_ранце) и столкнулся с проблемой сравнения двух решений, представленных вещественными числами.
— — — — — — — — — — Для текущего изложения я придумал несколько примеров, достаточно простых, но демонстрирующих суть проблемы. Несколько итераций позволят нам выйти на приемлемое решение обсуждаемой проблемы.
program Zero_01;
{$APPTYPE CONSOLE}
uses
  SysUtils,Math;
var argE: double;
    argR,arg: single;
    argR48: Real48;
begin
        arg:=Pi;      // Pi - Returns 3.1415926535897932385
                      // всего 20 значащих цифр
                      // у arg - тип single - 7 значащих цифр
        argR48:=Pi;
        writeln('Pi=',#9,Pi);
        writeln('arg=',#9,arg);
        writeln('argR48=',#9,argR48);
        writeln('__________');
        writeln(' sin=');
{1}     writeln(#9,sin(Pi));

{2}     writeln(#9,sin(arg));
        argE:=sin(arg);  // argE: double;

{3}     writeln(#9,argE);

          argR:=argE/1E30;
          argR:=argR*1E30;
{4}     writeln(#9,argR);

          argR:=argE/1E39;
          argR:=argR*1E39;
{5}     writeln(#9,argR);

        readln;
end.

Вот результаты этой программы:
Обработка малых чисел
Первое что мы увидим — число Пи представлено 20-ю значащими цифрами. Но мы то знаем, что оно иррациональное. Отсюда уже будут возникать ошибки в расчете тригонометрических функций.
Далее видим, что присвоение Pi переменной может еще ухудшить ситуацию: в переменной типа single несовпадения с Pi начинаются с 8-й значащей цифра, а для Real48 — с 13-ой;
Соответственно и вычисления синуса от разных Пи будут отличаться. Но и само значение sin(Pi) окажется отличным от нуля во всех случаях. Однако, отличие от нуля для значения Pi, представленного константой проявляется порядком E-20, а для значения Pi, представленного переменной типа single — Е-08. То есть в подсчет через переменную дает результат в 10^12 раз менее точно.

Основываясь на полученных результатах пойдем дальше. Раз уж вещественные числа так осложняют нам расчеты, попробуем их искусственно упростить, то есть округлить (снизить точность расчетов).

program Zero_02;
{$APPTYPE CONSOLE}
uses
  SysUtils,Math;
var arg: single;

function comparing1(arg1: single; arg2: single): boolean;
begin
  result:= abs(arg1-arg2)<1E-6;
end;

function comparing2(arg1: single; arg2: single): boolean;
var temp: single;
begin
  temp:=abs(arg1-arg2)/1E39;
  result:= temp*1E39 = 0;
end;

begin
        arg:=Pi;
        writeln('Pi=',#9,Pi);
        writeln('arg=',#9,arg);
        writeln('__________');
        writeln(' sin=');
{1}     writeln(#9,sin(2*Pi));
{2}     writeln(#9,sin(2*arg));
{3}     if comparing1(sin(2*arg),0)
                then writeln(#9,'=0')
                else writeln(#9,'<>0');
{4}     if comparing2(sin(2*arg),0)
                then writeln(#9,'=0')
                else writeln(#9,'<>0');

        readln;
end.

Вот результаты этой программы:
Обработка малых чисел
Как видите, константа Pi и переменная всё так же отличаются, как и результаты посчитанные для величины 2*Pi. Результаты всё так же отличаются не только друг от друга, но, что более печально, и от нуля. Однако в последних двух строчках программы нули таки выводятся. Тут я применил два подхода, обозначенных ещё в самом начале этого изложения по поводу ограничений на порядок числа и на точность, заданной мантиссой числа.
Функции сравнения comparing1 и comparing2 разными путями приходят к одинаковому результату: если числа отличаются ненамного, значит их можно признать равными. Только первая функция напрямую сравнивает, можно сказать, порядок несовпадения чисел, а вторая — отсекает все значащие цифры из мантиссы делением/умножением. Деление на порядок 39 выбрано, конечно, не случайно — надо ведь дотянуть до границы в Single 1.5x10^-45.
Итак, этот подход можно применять для работы с малыми числами или для сравнения двух малоотличающихся чисел.

Пойдём дальше и попробуем применить подход к какой-нибудь задаче. Я особо времени не тратил на изобретение задач, а взял задачу с заведомо известным результатом (для того чтобы заранее знать ответ).
Пусть нужно найти точки пересечения двух функций: cos(x/2) и sin(2*x), в диапазоне изменения x от 2,5 до 3,5. Нетрудно догадаться, что при x=Pi значения функций равны. Но только не для вычислителя — см. рисунок:
Обработка малых чисел
Обратите внимание, что косинус и синус для обозначенного значения в Excel не равны нулю.
Так же следует учитывать, что мы считаем значения с определенной дискретностью (шаг на рис. равен 0,01), то есть мы можем просто не попасть точно в точку пересечения функций. Как же быть?
Проведем первый эксперимент, имитируя Excel — просто запустим цикл с определенным, фиксированным шагом 1E-06:
program Zero_03;
{$APPTYPE CONSOLE}
uses SysUtils,Math;

function comparing1(arg1: single; arg2: single): boolean;
begin
  result:= abs(arg1-arg2)<1E-5;
end;

var low,hi,i,z1,z2: single;
    k: word;

begin
  low:=2; hi:=4;
  i:=low; k:=0;
  repeat
    z1:=cos(i/2);
    z2:=sin(2*i);
    if z1=z2 then writeln(#9,i);
    if comparing1(z1,z2)
              then
              begin
                writeln(#9,'z1=',z1);
                writeln(#9,'z2=',z2);
                inc(k);
              end;
    i:=i+1E-06;
  until i>hi;
  writeln(#9,k,' solutions');
  readln;
end.

Обратите внимание — я оставил строку «if z1=z2 then writeln(#9,i);» для того, чтобы удостовериться, что прямое сравнение вещественных чисел не даст желаемого результата.
Вот результаты работы программы:
Обработка малых чисел
Итак, прямой перебор всего диапазона даст нам 8 решений с заданной точностью. Видимо такой подход не оптимален, особенно по причине большого перебора, зависящего от величины диапазона и точности вычисления. Классическим способом преодоления указанной проблемы считается использование метода половинного деления с постепенным сужением диапазона.
program Zero_04;
{$APPTYPE CONSOLE}
uses SysUtils,Math;
const error=1E-6;
var low,hi,i,z1,z2: single;
    k: word;

function comparing1(arg1: single; arg2: single): boolean;
begin
  result:= abs(arg1-arg2)<error;
end;

begin
  low:=2; hi:=4; k:=0;
  // в начале на нижней границе  z1>z2
  // в начале на верхней границе z1<z2
  i:=(low+hi)/2;
  z1:=cos(i/2);
  z2:=sin(2*i);
  while not comparing1(z1,z2) do
  begin
    inc(k);
    i:=(low+hi)/2;
    z1:=cos(i/2);
    z2:=sin(2*i);
    if z1>z2
      then low:=i
      else hi:=i;
  end;
                writeln('Pi=',#9,Pi);
                writeln('i=',#9,i);
                writeln('z1=',#9,z1);
                writeln('z2=',#9,z2);
                writeln('k=',#9,k);                
  readln;
end.

Вот результаты работы программы:
Обработка малых чисел
Что видим: решение (точка пересечения двух функций) с заданной точностью (error) найдено за 21 шаг для аргумента близкого к Пи (отличия проявляются лишь в 8-ой значащей позиции мантиссы аргумента).

==========================================
summary
1) вещественные числа сравнивать операцией "=" не правильно;
2) переменные вещественного типа дают результат с ограниченной точностью;
3) если в задаче не указано с какой точностью найти решение, надо самому об этом позаботиться.
==========================================
P.S.1. Я мог что-то не заметить, что-то упустить — желающие могут оставить свои комментарии…
P.S.2. Программы апробировались на Delphi 7 и Delphi 2010.

Комментарии (1) свернуть  |  развернуть

  • avatar
  • belan
  • 09 декабря 2011, 20:15
+1
Конечно, было бы странным, что в такой продвинутой среде разработке приложений как Delphi не было бы предусмотрено функции сравнения вещественных чисел отличающихся совсем немного или сравнения очень маленького вещественного числа с нулем. Да, такая функция есть — SameValue, она переподгружаемая:
Delphi syntax:
function SameValue(const A, B: Single; Epsilon: Single = 0): Boolean; overload;
function SameValue(const A, B: Double; Epsilon: Double = 0): Boolean; overload;
function SameValue(const A, B: Extended; Epsilon: Extended = 0): Boolean; overload;

Первые два аргумента сравниваемые величины, третий аргумент — необязательный, он обозначает максимальную допустимую разницу между сравниваемыми величинами, до которой они признаются равными.
Используем её для реализации той же задачи, что представлена в последнем примере, заменив ею мою функцию comparing1:
program P_Eq_Zero_func_pribl_SameValue;
{$APPTYPE CONSOLE}
uses SysUtils,Math;
const error=1E-6;
var low,hi,i,z1,z2: single;
    k: word;

begin
  low:=2; hi:=4; k:=0;
  // в начале на нижней границе  z1>z2
  // в начале на верхней границе z1<z2
  i:=(low+hi)/2;
  z1:=cos(i/2);
  z2:=sin(2*i);
  while not SameValue(z1,z2,error) do
  begin
    inc(k);
    i:=(low+hi)/2;
    z1:=cos(i/2);
    z2:=sin(2*i);
    if z1>z2
      then low:=i
      else hi:=i;
  end;
                writeln('Pi=',#9,Pi);
                writeln('i=',#9,i);
                writeln('z1=',#9,z1);
                writeln('z2=',#9,z2);
                writeln('k=',#9,k);                
  readln;
end.

Вот результаты работы этой программы:Видно, что результаты один в один совпадают с программой где использовалась моя функция comparing1.
Пойдем дальше… Очевидно, что необязательный параметр, потому и необязательный, что его необязательно указывать, проведем эксперимент — я не буду опять всю программу вставлять, а только укажу измененную строку:
while not SameValue(z1,z2) do
Изменятся ли результаты работы программы?Видно, что количество шагов поиска решения сократилось до 11 и числа признаны равными с большей погрешностью — интересно почему?
Посмотрим в модуле Math как устроена функция SameValue
const
  FuzzFactor = 1000;
  ExtendedResolution = 1E-19 * FuzzFactor;
  DoubleResolution   = 1E-15 * FuzzFactor;
  SingleResolution   = 1E-7 * FuzzFactor;

function SameValue(const A, B: Extended; Epsilon: Extended): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := Max(Min(Abs(A), Abs(B)) * ExtendedResolution, ExtendedResolution);
  if A > B then
    Result := (A - B) <= Epsilon
  else
    Result := (B - A) <= Epsilon;
end;

function SameValue(const A, B: Double; Epsilon: Double): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := Max(Min(Abs(A), Abs(B)) * DoubleResolution, DoubleResolution);
  if A > B then
    Result := (A - B) <= Epsilon
  else
    Result := (B - A) <= Epsilon;
end;

function SameValue(const A, B: Single; Epsilon: Single): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := Max(Min(Abs(A), Abs(B)) * SingleResolution, SingleResolution);
  if A > B then
    Result := (A - B) <= Epsilon
  else
    Result := (B - A) <= Epsilon;
end;
Выбор подгружаемой функции и, соответсвенно погрешность, зависит от типа входных значений.
— — — — — — — Кстати, там же находим функцию сравнения числа с нулем IsZero, аналогично устроеную:
function IsZero(const A: Extended; Epsilon: Extended): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := ExtendedResolution;
  Result := Abs(A) <= Epsilon;
end;

function IsZero(const A: Double; Epsilon: Double): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := DoubleResolution;
  Result := Abs(A) <= Epsilon;
end;

function IsZero(const A: Single; Epsilon: Single): Boolean;
begin
  if Epsilon = 0 then
    Epsilon := SingleResolution;
  Result := Abs(A) <= Epsilon;
end;
Может пригодиться…

Прокомментировать