Двухминутка дзена

[Paul Kellerman]

Вот казалось бы, есть известная формулка для вычисления вероятности
наступление события "хотя бы один раз" в серии из m независимых испы-
таний с заданной вероятностью p наступления этого события: 1-(1-p)^m.
Да и запрограммировать формулу тоже казалось бы - одна строчка кода.
p имеет тип 80-битного вещественного, m имеет тип 64-битного целого.

Но не тут-то было. А что если p очень маленькие, например 10^(-50) при
которой вычисление 1-(1-p)^m в 80-битной арифметике (тип p: extended)
с плавающей запятой даст 0 в силу невозможности представления разно-
сти (1-10^(-50)) в 80-битном вещественном типе и округлением ее до 1.

Тут конечно сразу вспоминаем арифметику многократной точности, где
вещественные числа аппроксимируются рациональными числами, в кот.
в числителе и знаменателе программно обрабатываемые сверхдлинные
целые числа, содержащие миллионы, миллиарды и более битов. Однако,
как-то прикручивать в программе библиотеку с функциями арифметики
многократной точности ради одной единственной формулы не очень как-
то хочется. Да и еще есть нюанс: если, например, число p = 10^(-1000),
а m = 10^18 (m - 64-битной целое число), то там будет получаться раци-
ональное число, содержащее целые числа с (10^21) значащими цифрами.
И для хранения каждого целого понадобится ~ 3,32 * 10^21 бит памяти.
Как-то многовато... ~ 377659500 терабайт памяти на каждое из чисел.
Пример, конечно, запредельный, но все же вполне допустимый при 80-
битном вещественном p и 64-битном целом m. Так что думаем дальше.

Вспоминаем матанализ... 1-(1-p)^m = 1-exp(m*ln(1-p)) и в случае когда
p очень малое, ln(1-p) ~ -p, и, соответственно, 1-(1-p)^m ~1-exp(-m*p).
Но надо заметить, что когда само произведение m*p слишком малое, а
именно менее 10^(-19), то в 80-битной арифметики вычисление экспо-
ненты тоже начинает округляться до 1, и в итоге снова 1-exp(-m*p) = 0.
Но при очень малых самих m*p опять же вспоминаем матанализ: exp(-x) ~
1-x+(x^2)/2. Но тут тогда уж обратиться к самой 1-(1-p)^m и раскрыть
скобки, вытащив оттуда первые два (достаточно) слагаемых: 1-(1-p)^m ~
1-(1-m*p+m*((m-1)/2)*p^2) = m*p - m*((m-1)/2)*p^2. Теперь остается
понять, что считать "очень малым" в 80-битной арифметике, сравнивая с
кот. p и m*p, мы будем решать, когда по какой формуле лучше считать.
Опытным путем я подобрал, что "очень малым" считать все что < 10^(-9).
Это та граница, при которой обеспечивается наиболее плавный переход
результатов вычислений от формуле к формуле с мин. погрешностями.

В итоге вот такой код я нарисовал (Delphi) для вычисления 1-(1-p)^m:

unit Unit2;

interface

uses
  Math;
const
  fl_zero: Extended = 0;
  fl_one: Extended = 1;
  fl_eps: Extended = 1E-9;

function Calc(m: Int64; p: Extended): Extended;

implementation

function Calc;
begin
  if (p = fl_zero) then Result:= fl_zero;
  if (p = fl_one) then Result:= fl_one;
  if ((p > fl_zero) and (p < fl_one)) then
    begin
      if (m = 1) then Result:= p
        else if (p > fl_eps)
               then Result:= fl_one - Power((fl_one - p), m)
               else if (p * m > fl_eps)
                      then Result:= fl_one - exp(-p*m)
                      else Result:= p*m - p*p*m*(m-1)/2;
    end;
end;

end.

[mitek1989]

Сообщение от Paul Kellerman

если, например, число p = 10^(-1000),
а m = 10^18

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