pidigits.pas

{ The Computer Language Benchmarks Game
  http://benchmarksgame.alioth.debian.org

  contributed by Vincent Snijders
}

{$mode objfpc}


program pidigits;

type
  { TBigInt }

  DWord = 0..4294967295;
  PDWord = ^DWord;

  PBigInt = ^TBigInt;

  { TBigInt }

  TBigInt = class
  private
    Digit: pdword;
    FSize: integer;
    Capacity: integer;
    FNextFree: PBigInt; (* used to maintain the freelist *)
    procedure Clear;
    procedure Resize(NewSize: integer);
    function IsNegative: boolean; inline;
    function IsZero: boolean; inline;
    procedure Negate;
  public
    procedure Create(InitialSize: integer);
    procedure Destroy; override;
    function GetDigit(i: integer) : DWord; inline;
  end;

type
  TBigMatrix = array[1..2,1..2] of TBigInt;
  TIntMatrix = array[1..2,1..2] of integer;

var
  BigIntStack: PBigInt;
  BigIntStackLen: integer;
  BigIntTop: integer;
  FirstFreeBigInt: TBigInt;

{ BigInt memory management }

procedure FreeBigInts;
var
  Next: TBigInt;
begin
  while assigned(FirstFreeBigInt) do begin
    Next := FirstFreeBigInt.FNextFree;
    FirstFreeBigInt.Free;
    FirstFreeBigInt := Next;
  end;
end;

function GetBigInt(Size: integer; DoClear: boolean = true) : TBigInt;
var
  Current, Previous: TBigInt;
begin
  if assigned(FirstFreeBigInt) then begin
    Previous := nil;
    Current := FirstFreeBigInt;
    repeat
      if (Current.Capacity>=Size) then begin
        Result:=Current;
        Result.FSize:= Size;
        if DoClear then
          Result.Clear;
        if assigned(previous) then
          Previous.FNextFree := Current.FNextFree
        else
          FirstFreeBigInt := Current.FNextFree;
        exit;
      end;
      Previous := Current;
      Current := Current.FNextFree;
    until Current=nil;
    Result := FirstFreeBigInt;
    Result.Resize(Size);
    FirstFreeBigInt := FirstFreeBigInt.FNextFree;
  end
  else result := TBigInt.Create(Size);
end;

function GetBigInt(bi: TBigInt) : TBigInt; inline;
begin
  result := GetBigInt(bi.FSize, false);
  Move(bi.Digit^, Result.Digit^, bi.FSize*sizeof(dword));
end;

procedure FreeBigInt(bi: TBigInt);
begin
  bi.FNextFree := FirstFreeBigInt;
  FirstFreeBigInt := bi;
end;

{ TBigInt }

operator := (i: integer) : TBigInt; inline;
begin
  Result := GetBigInt(1);
  Result.Digit[0] := dword(i);
end;

constructor TBigInt.Create(InitialSize: integer);
begin
  FSize:= InitialSize;
  Capacity:= 2*FSize;
  GetMem(Digit, Capacity*sizeof(DWord));
  Clear;
end;

destructor TBigInt.Destroy;
begin
  FreeMem(Digit);
  inherited Destroy;
end;

procedure TBigInt.Clear;
begin
  FillChar(Digit[0], FSize*sizeof(DWord), 0);
end;

procedure TBigInt.Resize(NewSize: integer);
begin
  FreeMem(Digit);
  FSize:= NewSize;
  Capacity:= 2*FSize;
  GetMem(Digit, Capacity*sizeof(DWord));
  Clear;
end;

function TBigInt.IsNegative: boolean; inline;
begin
  result := (Digit[FSize-1] and $80000000)>0;
end;

function TBigInt.IsZero:boolean;inline;
begin
  result := (FSize=1) and (Digit^=0);
end;

procedure TBigInt.Negate;
var
  value: int64;
  valueparts : array[0..1] of dword absolute value;
  carry: integer;
  CurDigit: PDWord;
begin
  if IsZero then exit;
  CurDigit:= @Digit[FSize-1];
  repeat
    CurDigit^:= not CurDigit^;
    dec(CurDigit);
  until CurDigit<Digit;
  carry := 1;
  CurDigit:=Digit;
  while (carry>0) do begin
    value := CurDigit^;
    inc(value);
    CurDigit^ := valueparts[0];
    carry := valueparts[1];
    inc(CurDigit);
  end;
end;

function TBigInt.GetDigit(i: integer): DWord; inline;
begin
  if (i<FSize) then
    result := Digit[i]
  else if IsNegative then
    result := $FFFFFFFF
  else
    result := 0;
end;

{ BigInt Calculation }

procedure PushBigInt(bi: TBigInt);
begin
  inc(BigIntTop);
  if (BigIntTop=BigIntStackLen) then
    RunError(1025); // not implemented, too complicated calculation
  BigIntStack[BigIntTop]:=bi;
end;

procedure PushBigIntByValue(bi: TBigInt);
begin
  inc(BigIntTop);
  if (BigIntTop=BigIntStackLen) then
    RunError(1025); // not implemented, too complicated calculation
  BigIntStack[BigIntTop]:= GetBigInt(bi);
end;

function PopBigInt: TBigInt;
begin
  result:=BigIntStack[BigIntTop];
  dec(BigIntTop);
end;

procedure BigIntAdd;
var
  a, b: TBigInt;
  bSignExtend: dword;
  Result: TBigInt;
  carry: integer;
  sum: int64;
  maxsize, minsize, i: integer;
  sumparts : array[0..1] of integer absolute sum;
  aDigit, bDigit, ResultDigit: PDWord;
begin
  if BigIntStack[BigIntTop-1].FSize<BigIntStack[BigIntTop].FSize then begin
    a:= BigIntStack[BigIntTop];
    b:= BigIntStack[BigIntTop-1];
  end
  else begin
    a:= BigIntStack[BigIntTop-1];
    b:= BigIntStack[BigIntTop];
  end;
  if b.IsZero then
    Result := a
  else begin
    maxsize:=a.FSize;
    minsize:=b.FSize;
    Result := GetBigInt(maxsize+1);
    carry := 0;
    aDigit:= a.Digit; bDigit:= b.Digit; ResultDigit:= Result.Digit;
    for i:= 0 to minsize-1 do begin
      sum := int64(aDigit^) + int64(bDigit^) + carry;
      carry := sumparts[1];
      ResultDigit^ := sumparts[0];
      inc(aDigit); inc(bDigit); inc(ResultDigit);
    end;
    if b.IsNegative then
      bSignExtend := $FFFFFFFF
    else
      bSignExtend := 0;
    for i:= minsize to maxsize do begin
      sum := int64(a.GetDigit(i)) + bSignExtend + carry;
      carry := sumparts[1];
      ResultDigit^ := sumparts[0];
      inc(ResultDigit);
    end;
    while (Result.FSize>1) and (Result.Digit[Result.FSize-1]=0)
      and (Result.Digit[Result.FSize-2] and $80000000=0) do
      dec(Result.FSize);
    while (Result.FSize>1) and
       (Result.Digit[Result.FSize-1]=$FFFFFFFF) and
       (Result.Digit[Result.FSize-2] and $80000000>0) do
      dec(Result.FSize);
    FreeBigInt(a);
  end;
  FreeBigInt(b);
  dec(BigIntTop);
  BigIntStack[BigIntTop]:=Result;
end;

procedure BigIntMulInt(int: integer);
type
  TWordPart = record
    w1, w2: word;
  end;
var
  mcarry: dword;
  value: qword;
  valueparts : array[0..1] of dword absolute value;
  BiNeg, IntNeg: boolean;
  i:Integer;
  TopBi, Result: TBigInt;
  TopBiDigit, ResultDigit: PDWord;
begin
  TopBi := BigIntStack[BigIntTop];
  if (int=0) or (TopBi.IsZero) then begin
    TopBi.FSize := 1;
    TopBi.Digit[0]:=0;
  end
  else begin
    BiNeg := TopBi.IsNegative;
    if BiNeg then
      TopBi.Negate;
    IntNeg := int<0;
    if IntNeg then
      int := -int;
    Result := GetBigInt(TopBi.FSize+1, false);
    mcarry := 0;
    TopBiDigit := TopBi.Digit;
    ResultDigit := Result.Digit;
    if (int and $FFFF0000)=0 then
      for i:= 0 to Result.FSize-2 do begin
        {this is what I want to do, but to get to the carry fpc compiles it into
        an expensive qword*qword mulitplication:   }
        {value := qword(TopBiDigit^) * int + mcarry;}
        value := TWordPart(TopBiDigit^).w1 * word(int) +
          qword(TWordPart(TopBiDigit^).w2 * word(int)) shl 16 + mcarry;
        ResultDigit^ := valueparts[0];
        mcarry := valueparts[1];
        inc(TopBiDigit); inc(ResultDigit);
      end
    else
    // this branch is less often taken, so no hand code dword * dword multiplication
      for i:= 0 to Result.FSize-2 do begin
        value := qword(TopBiDigit^) * int + mcarry;
        ResultDigit^ := valueparts[0];
        mcarry := valueparts[1];
        inc(TopBiDigit); inc(ResultDigit);
      end;
    ResultDigit^ := mcarry;
    while (Result.FSize>1) and (Result.Digit[Result.FSize-1]=0) and
      ((Result.Digit[Result.FSize-2] and $80000000)=0) do
      dec(Result.FSize);
    if (BiNeg<>IntNeg) then
      Result.Negate;
    FreeBigInt(TopBi);
    BigIntStack[BigIntTop]:=Result;
  end;
end;

function BigIntDivIntResult: integer;
var
  dividend: TBigInt;
  divisor: TBigInt;
  carry: dword;
  diff: int64;
  diffparts: array[0..1] of dword absolute diff;
  i: integer;
  DividendDigit: PDWord;
  DivisorDigit: PDWord;

  function DividendIsSmallerThanDivisor : boolean; inline;
  var
    i: integer;
  begin
    while (Dividend.FSize>1) and (Dividend.Digit[Dividend.FSize-1]=0)
      and (Dividend.Digit[Dividend.FSize-2] and $80000000=0) do
      dec(Dividend.FSize);
    if dividend.FSize=divisor.FSize then begin
      i := dividend.FSize-1;
      while (i>=0) and (dividend.Digit[i]=divisor.Digit[i]) do
        dec(i);
      Result:= (i>=0) and (dividend.Digit[i]<divisor.Digit[i]);
    end
    else
      Result:=dividend.FSize<divisor.FSize;
  end;

begin
  dividend := BigIntStack[BigIntTop-1];
  divisor := BigIntStack[BigIntTop];
  Result:=0;
  while not DividendIsSmallerThanDivisor do begin
    inc(Result);
    carry := 0;
    DividendDigit := Dividend.Digit;
    DivisorDigit := Divisor.Digit;
    for i:= 0 to divisor.FSize-1 do begin
      diff := int64(dividendDigit^) - (divisorDigit^ + carry);
      carry := diffparts[1] and $1;
      dividendDigit^ := diffparts[0];
      inc(DividendDigit); inc(DivisorDigit);
    end;
    for i:= divisor.FSize to dividend.FSize-1 do begin
      diff := int64(dividendDigit^) - (divisor.GetDigit(i) + carry);
      carry := diffparts[1] and $1;
      dividendDigit^ := diffparts[0];
      dividend.Digit[i] := diffparts[0];
      inc(DividendDigit);
    end;
  end;
  FreeBigInt(dividend);
  FreeBigInt(divisor);
  dec(BigIntTop,2);
end;

procedure Init;
begin
  BigIntStackLen := 8;
  GetMem(BigIntStack, BigIntStackLen * sizeof(TBigInt));
  BigIntTop := -1;
  FirstFreeBigInt := nil;
end;

procedure Finalize;
begin
  Freemem(BigIntStack);
  FreeBigInts;
end;

{ Matrix manipulation }

procedure FreeBigIntMatrix(a: TBigMatrix); inline;
begin
  FreeBigInt(a[1,1]);
  FreeBigInt(a[1,2]);
  FreeBigInt(a[2,1]);
  FreeBigInt(a[2,2]);
end;

function DotProduct(a1,a2: TBigInt; b1,b2: integer; FreeBigInt: boolean) : TBigInt; inline;
begin
  if FreeBigInt then
    PushBigInt(a1)
  else
    PushBigIntByValue(a1);
  BigIntMulInt(b1);
  if FreeBigInt then
    PushBigInt(a2)
  else
    PushBigIntByValue(a2);
  BigIntMulInt(b2);
  BigIntAdd;
  Result:= PopBigInt;
end;

operator * (a: TBigMatrix; b : TIntMatrix) : TBigMatrix;
begin
  result[1,1] := DotProduct(a[1,1],a[1,2], b[1,1], b[2,1], false);
  result[1,2] := DotProduct(a[1,1],a[1,2], b[1,2], b[2,2], true);
  result[2,1] := DotProduct(a[2,1],a[2,2], b[1,1], b[2,1], false);
  result[2,2] := DotProduct(a[2,1],a[2,2], b[1,2], b[2,2], true);
end;

operator * (a: TIntMatrix; b : TBigMatrix) : TBigMatrix;
begin
  result[1,1] := DotProduct(b[1,1],b[2,1],a[1,1],a[1,2], false);
  result[1,2] := DotProduct(b[1,2],b[2,2],a[1,1],a[1,2], false);
  result[2,1] := DotProduct(b[1,1],b[2,1],a[2,1],a[2,2], true);
  result[2,2] := DotProduct(b[1,2],b[2,2],a[2,1],a[2,2], true);
end;

function InitBigMatrix(a,b,c,d: integer): TBigMatrix;
begin
  result[1,1] := a;
  result[1,2] := b;
  result[2,1] := c;
  result[2,2] := d;
end;

function InitIntMatrix(a,b,c,d: integer): TIntMatrix; inline;
begin
  result[1,1] := a;
  result[1,2] := b;
  result[2,1] := c;
  result[2,2] := d;
end;

{ calculating pidigits}

procedure PrintPiDigits(const NumDigits: integer);
var
  n: integer = 0;
  k: integer = 0;
  z: TBigMatrix;
  x,p: TIntMatrix;
  Digit: integer;

  function Extract(x:integer) : integer;
  begin
    PushBigIntByValue(z[1,1]);
    BigIntMulInt(x);
    PushBigIntByValue(z[1,2]);
    BigIntAdd;
    PushBigIntByValue(z[2,1]);
    BigIntMulInt(x);
    PushBigIntByValue(z[2,2]);
    BigIntAdd;
    result := BigIntDivIntResult;
  end;
  
  function GetDigit : integer;
  begin
    result := Extract(3);
  end;

  function IsSafe : boolean;
  begin
    result := Digit = Extract(4);
  end;

  procedure Produce;
  begin
    p[1,2] := -10 * digit;
    z := p * z;
  end;

  procedure Consume;
  begin
    inc(k);
    x[1,1] := k;
    x[1,2] := 4*k+2;
    x[2,2] := 2*k+1;
    z:= z * x;
  end;

begin
  z := InitBigMatrix(1, 0, 0, 1);
  p := InitIntMatrix(10, 0, 0, 1);
  x[2,1] := 0;
  while (n<NumDigits) do begin
    Digit := GetDigit;
    while not IsSafe do begin
      Consume;
      Digit:= GetDigit;
    end;
    Produce;
    write(Digit);
    inc(n);
    if (n mod 10)=0 then
      writeln(#9':', n);
  end;
  FreeBigIntMatrix(z);
end;

var
  n: integer;
  errorcode: integer;

begin
  Init;
  if (ParamCount=1) then begin
    val(ParamStr(1), n, errorcode);
    PrintPiDigits(n);
  end;
  Finalize;
end.