19番会議室のユーザ談話室で話題になっていたペプシ算を計算するプログラムで
す。

「ペプシ算」（(C)1999 凛, 初出 nifty:FDELPHI/MES/19/1992）とは、次のよう
なものです。

  清涼飲料水のキャップのオマケとして、10 種類の清涼飲料水マンのキャッ
  プがついてくるキャンペーンがありました。そこで、豪気な凛さんは 24 本
  の箱買いを実行しました。24 本も買えば、10 種類そろっていそうな感じで
  すが、実際には（以下省略）そこで、そろわない場合の確率はどのくらいに
  なるのか？果たして凛さんは、全部の種類の清涼飲料水マンのキャップを
  持っているのだろうか？

10種類を24本購入した場合の、重複順列の場合の数は10Π24 = 10^24 通りで
す。以下のプログラムでは、10 種類全部そろっている場合の数を計算します。
これを 10^24 で割れば、10 種類全部そろっている確率が求められます。計算の
アルゴリズムは、Taka2 さんの nifty:FDELPHI/MES/19/2009 によります。

なお、計算の途中で桁あふれが発生するため、多倍長整数の演算を行っていま
す。多倍長計算の出典は、「コンピュータアルゴリズム事典」（1987 奥村晴彦,
技術評論社）です。より詳しい理論などは、Dr.Knuth の "The Art of Computer
Programming" などにあります。



program Pepsi;

uses SysUtils;

const
  Base = 100;
  BaseMinusOne = 99;
  Base2 = 10000;
  MaxLen = 500;

type
  TDigit = 0..BaseMinusOne;

  TANumber = record
    Len: 0..MaxLen;
    Num: array[0..MaxLen] of TDigit;
  end;

procedure LongAdd(var X, A, B: TANumber);
procedure LongSub(var X, A, B: TANumber);
procedure ShortMul(var X, A: TANumber; N: TDigit);
procedure LongMul(var X, A, B: TANumber);
procedure LongAssign(var X: TANumber; I: Integer);

// 多倍長の整数の加算 X = A + B
procedure LongAdd(var X, A, B: TANumber);

  procedure Aux(var X, A, B: TANumber);
  var
    I, T: Integer;
  begin
    T := 0;
    for I := 1 to B.Len do
    begin
      T := A.Num[I] + B.Num[I] + T;
      if T < Base then
      begin
        X.Num[I] := T;
        T := 0;
      end
      else
      begin
        X.Num[I] := T - Base;
        T := 1;
      end;
    end;
    I := B.Len + 1;
    while (I <= A.Len) and (T <> 0) do
    begin
      T := A.Num[I] + T;
      if T < Base then
      begin
        X.Num[I] := T;
        T := 0;
      end
      else
      begin
        X.Num[I] := T - Base;
        T := 1;
      end;
      Inc(I);
    end;
    for I := I to A.Len do
      X.Num[I] := A.Num[I];
    if T = 0 then
      X.Len := A.Len
    else if A.Len < MaxLen then
    begin
      X.Len := A.Len + 1;
      X.Num[X.Len] := 1;
    end
    else
      raise Exception.Create('桁あふれが発生しました');
  end;

begin
  if A.Len >= B.Len then
    Aux(X, A, B)
  else
    Aux(X, B, A);
end;

// 多倍長の整数の減算 X = A - B （ただし A >= B）
procedure LongSub(var X, A, B: TANumber);
var
  I, T: Integer;
begin
  if A.Len >= B.Len then
  begin
    T := 0;
    for I := 1 to B.Len do
    begin
      T := A.Num[I] - B.Num[I] - T;
      if T >= 0 then
      begin
        X.Num[I] := T;
        T := 0;
      end
      else
      begin
        X.Num[I] := T + Base;
        T := 1;
      end;
    end;
    I := B.Len + 1;
    while (I <= A.Len) and (T <> 0) do
    begin
      T := A.Num[I] - T;
      if T >= 0 then
      begin
        X.Num[I] := T;
        T := 0;
      end
      else
      begin
        X.Num[I] := T + Base;
        T := 1;
      end;
      Inc(I);
    end;
    for I := I to A.Len do
      X.Num[I] := A.Num[I];
    X.Len := A.Len;
    X.Num[0] := 1;
    while X.Num[X.Len] = 0 do
      X.Len := X.Len - 1;
  end
  else
    T := 1;
  if T <> 0 then
    raise Exception.Create('結果が負になりました');
end;

// 小さな整数と多倍長整数の乗算 X = A * N
procedure ShortMul(var X, A: TANumber; N: TDigit);
var
  I, T: Integer;
begin
  if N > 0 then
  begin
    T := 0;
    for I := 1 to A.Len do
    begin
      T := A.Num[I] * N + T;
      X.Num[I] := T mod Base;
      T := T div Base;
    end;
    if T = 0 then
      X.Len := A.Len
    else if A.Len = MaxLen then
      raise Exception.Create('桁あふれが発生しました')
    else
    begin
      X.Len := A.Len + 1;
      X.Num[X.Len] := T;
    end;
  end
  else
    X.Len := 0;
end;

// 多倍長整数の乗算 X = A * B
procedure LongMul(var X, A, B: TANumber);
var
  I, J, L, T, U: Integer;
begin
  if (A.Len <> 0) and (B.Len <> 0) then
  begin
    if A.Len + B.Len > X.Len then
      L := A.Len + B.Len
    else
      L := X.Len + 1;
    if L <= MaxLen then
    begin
      for I := X.Len + 1 to L do
        X.Num[I] := 0;
      for J := 1 to B.Len do
      begin
        U := B.Num[J];
        if U <> 0 then
        begin
          T := 0;
          for I := J to A.Len + J - 1 do
          begin
            T := X.Num[I] + A.Num[I - J + 1] * U + T;
            X.Num[I] := T mod Base;
            T := T div Base;
          end;
          I := A.Len + J;
          T := X.Num[I] + T;
          while T >= Base do
          begin
            X.Num[I] := T - Base;
            Inc(I);
            T := X.Num[I] + 1;
          end;
          X.Num[I] := T;
        end;
      end;
      if X.Num[L] = 0 then
        X.Len := L - 1
      else
        X.Len := L
    end
    else
      raise Exception.Create('桁あふれが発生しました')
  end
end;

// 整数を多倍長整数に代入 X = I
procedure LongAssign(var X: TANumber; I: Integer);
begin
  if I < 0 then
    raise Exception.Create('負の数は扱えません')
  else
    with X do
    begin
      Len := 0;
      while I <> 0 do
      begin
        Inc(Len);
        Num[Len] := I mod Base;
        I := I div Base;
      end;
    end;
end;

// ペプシ算の本体
function f(k, n:TDigit):TANumber;
var
  A,B,C:TANumber;
begin
  if k = 1 then
    LongAssign(Result, 1)
  else if n = 1 then
    LongAssign(Result, 0)
  else
  begin
    A := f(k, n - 1);
    ShortMul(B, A, k);
    A := f(k - 1, n - 1);
    ShortMul(C, A, k);
    LongAdd(Result, B, C);
  end;
end;

// プログラムメイン（10種類24本でペプシ算を行い、結果を出力）

var
  I: Integer;

begin
  with f(10, 24) do
  begin
    if Len > 0 then
    begin
      write(Num[Len] : 1);
      for I := Len - 1 downto 1 do
        write(Num[I] div 10 : 1, Num[I] mod 10 : 1);
    end
    else
      write(0:1);
  end;
end.

                  99/4/14(Wed) 06:00am [AirCraft開発] PFF01344 DUDE

Original document by DUDE            氏 ID:(PFF01344)

ここにあるドキュメントは NIFTY SERVEの Delphi Users' Forum の16番会議室「玉石混淆みんなで作るSample蔵」に投稿されたサンプルです。これらのサンプルはボーランド株式会社がサポートする公式のものではありません。また、必ずしも動作が検証されているものではありません。これらのサンプルを使用したことに起因するいかなる損害も投稿者、およびフォーラムスタッフはその責めを負いません。使用者のリスクの範疇でご使用下さい。