{$define STANDALONE}
{$ifdef STANDALONE}
program EXTARITH; {demonstration of extended arithmetic}
{$N+,R-}
{$else STANDALONE}

unit EXTARITH;

INTERFACE
{$endif STANDALONE}

const digits=200; notsig=-99;  exact=20000; digitset=['0'..'9'];
type number=record
       sigdigits :integer;
       sign      :-1..1; {need special flag for zero}
       decpoint  :integer;
       s         :array[1..digits] of char;
       end;

var one,zero :number;

{$ifndef STANDALONE}
function NumString(var a :number):string;
function NumValue(var a :number):double;
procedure StringToNum({not var!}v :string; var a :number);
procedure IntToNum(x :integer; var a :number);
function Compare(var a,b :number):shortint;
function Equal(var a,b :number):boolean;
procedure Add(a,b :number; var c :number);
procedure Subtract(a,b :number; var c :number);
procedure Multiply(a,b :number; var c :number);
procedure IntegerPower(a,n :number; var c :number);
procedure Factorial(a :number; var c :number);
procedure FastMultiply(a,b :number; var c :number);
procedure Divide(a,b :number; var c :number);
procedure SqRoot(x :number; var y :number);

IMPLEMENTATION
{$endif STANDALONE}

procedure Trim(var s :string);
          begin
          while (s[1]=' ') do Delete(s,1,1);
          while (s[length(s)]=' ') do dec(s[0]);
          end;

function AtoI(a :char):shortint;
         begin
         if (a>'9') then AtoI:=notsig
         else AtoI:=ord(a)-48;
         end;

function max(a,b :integer):integer;
         begin
         if (a>b) then max:=a
         else max:=b;
         end;

function min(a,b :integer):integer;
         begin
         if (a>b) then min:=b
         else min:=a;
         end;

function NumValue(var a :number):double;
         var result     :double;
             i          :integer;
             ws,expstr  :string;
         begin with a do begin
         Str(decpoint,expstr);
         result:=0;
         i:=1;
         if (sign=1) then ws:='.'
         else ws:='-.';
         repeat
           ws:=ws+s[i];
           inc(i);
         until (i>16) or (s[i]>'9');
         ws:=ws+'E'+expstr;
         Val(ws,result,i);
         if (i<>0) then
           if (pos('x',ws)<=3) then result:=0
           else writeln('Error in evaluating number.');
         NumValue:=result;
         end;end;

function NumString(var a :number):string;
         var result :string;
             i      :integer;
         begin with a do begin
         result[0]:=char(succ(digits));
         if (decpoint<0) then begin
            result[1]:='.';
            for i:=1 to -decpoint do result[1+i]:='0';
            for i:=1 to digits+decpoint-1 do result[i-decpoint+1]:=a.s[i];
            end
         else begin
           for i:=1 to decpoint do result[i]:=s[i];
           result[succ(decpoint)]:='.';
           for i:=succ(decpoint) to digits do result[succ(i)]:=s[i];
           end;
         if (pos('.',result)>0) then
           while (result[length(result)]>'9') or (result[length(result)]<'0') do dec(result[0])
         else
           while (result[length(result)]>'9') or (result[length(result)]<='0') do dec(result[0]);
         if (pos('.',result)>0) then while (result[length(result)]='0') do dec(result[0]);
         if (sign=1) then NumString:=result
         else NumString:='-'+result;
         end;end;

procedure StringToNum({not var!}v :string; var a :number);
          var i,zeros,exponent,ta :integer;
              ss                  :string;
          begin with a do begin
          Trim(v);
          i:=pos('E',v);
          if (i=0) then exponent:=0
          else begin
            ss:=copy(v,succ(i),length(v)-i);
            if (ss[1]='+') then Delete(ss,1,1);
            Val(ss,exponent,ta);
            v:=copy(v,1,pred(i));
            end;
          if (v[1]='-') then begin sign:=-1; Delete(v,1,1); end
          else sign:=1;
          while (v[1]='0') and (length(v)>0) do Delete(v,1,1);
          decpoint:=pos('.',v);
          if (decpoint=0) then decpoint:=length(v)
          else begin Delete(v,decpoint,1); dec(decpoint); end;
          for i:=1 to length(v) do s[i]:=v[i];
          for i:=succ(length(v)) to digits do s[i]:='0';
          if (decpoint=0) then begin
            zeros:=1;
            while (zeros<=digits) and (s[zeros]='0') do inc(zeros);
            dec(zeros);
            if (zeros=digits) then begin a:=zero; exit; end
            else if (zeros>0) then begin
              for i:=1 to digits-zeros do s[i]:=s[i+zeros];
              decpoint:=-zeros;
              end;
            end
          else zeros:=0;
          for i:=succ(length(v))-zeros to digits do s[i]:='0';
          decpoint:=decpoint+exponent;
          a.sigdigits:=exact; {! I don't know why this is necessary, but you lose accuracy if you change it.}
          end;end;

procedure IntToNum(x :integer; var a :number);
          var v :string;
          begin
          Str(x,v);
          StringToNum(v,a);
          end;

procedure DoubleToNum(x :double; var a :number);
          var v :string;
          begin
          Str(x,v);
          StringToNum(v,a);
          end;

procedure Normalize(var c :number);
          var i,ci :integer;
              r :char;
          begin
          i:=1;
          while (c.s[i]='0') and (i<digits) do inc(i);
          if (i=digits) then begin c:=zero; exit; end;
          dec(i);
          if (c.sigdigits=exact) then r:='0' else r:='x';
          if (i>0) then begin
            dec(c.decpoint,i);
            if (c.sigdigits=exact) then
            else dec(c.sigdigits,i);
            for ci:=1 to digits-i do c.s[ci]:=c.s[ci+i];
            for ci:=succ(digits-i) to digits do c.s[ci]:=r;
            end;
          if (abs(c.sign)>1) then
            writeln('Sign greater than one error.');
          end;

function Compare(var a,b :number):shortint;
         var i :integer;
         begin
         if (a.sign>b.sign) then begin Compare:=1; exit; end;
         if (b.sign>a.sign) then begin Compare:=-1; exit; end;
         if (a.decpoint>b.decpoint) then begin Compare:=a.sign; exit; end;
         if (b.decpoint>a.decpoint) then begin Compare:=-a.sign; exit; end;
         i:=1;
         repeat
           if (a.s[i]>b.s[i]) then begin Compare:=a.sign; exit; end;
           if (b.s[i]>a.s[i]) then begin Compare:=-a.sign; exit; end;
           inc(i);
         until (i>digits) or (a.s[i]='x') or (b.s[i]='x');
         Compare:=0;
         end;

function Equal(var a,b :number):boolean;
         begin
         Equal:=(Compare(a,b)=0);
         end;

procedure AddDigits(var x,y,z :char; var carry :integer);
          begin
          if (x>'9') or (y>'9') then begin z:='x'; carry:=0; exit; end;
          z:=char(ord(x)+ord(y)-48+carry);
          if (z>'9') then begin
            carry:=1; z:=char(ord(z)-10); end
          else carry:=0;
          end;

procedure MultiplyDigits(var x,y,z :char; var carry :integer);
          var result :integer;
          begin
          if (x>'9') or (y>'9') then begin z:='x'; carry:=0; exit; end;
          result:=(ord(x)-48)*(ord(y)-48)+carry;
          carry:=result div 10;
          z:=char(48+(result mod 10));
          end;

procedure SubtractDigits(var x,y,z :char; var borrow :integer);
          begin
          if (x>'9') or (y>'9') then begin z:='x'; borrow:=0; exit; end;
          z:=char(48+ord(x)-ord(y)-borrow);
          if (z<'0') then begin
            borrow:=1; z:=char(ord(z)+10); end
          else borrow:=0;
          end;

{$ifdef STANDALONE}
procedure Subtract(a,b :number; var c :number); forward;
{$endif STANDALONE}

procedure NumberTimesDigit(var a :number; b :char; var c :number);
          var ai,bi,ci,carry :integer;
              as,bs,cs       :char;
          begin
          if (b='0') or (a.sign=0) then begin c:=zero; exit; end;
          if (b='1') then begin c:=a; exit; end;
          carry:=0;
          MultiplyDigits(a.s[digits],b,cs,carry);
          if (a.sigdigits=exact) and (cs='0') then c.sigdigits:=exact
          else c.sigdigits:=a.sigdigits;
          ci:=digits;
          for ai:=digits-1 downto 1 do begin
            MultiplyDigits(a.s[ai],b,c.s[ci],carry);
            dec(ci);
            end;
          if (carry=0) then begin
            for ci:=1 to pred(digits) do c.s[ci]:=c.s[succ(ci)];
            c.s[digits]:=cs;
            c.decpoint:=a.decpoint;
            end
          else begin
             c.s[1]:=char(carry+48);
             c.decpoint:=succ(a.decpoint);
             end;
          c.sign:=a.sign;
          end;

procedure Add(a,b :number; var c :number);
          var ai,bi,ci :integer;
              as,bs    :char;
              carry    :integer;
          begin
          if (a.sign=0) then begin c:=b; exit; end;
          if (b.sign=0) then begin c:=a; exit; end;
          if (a.sign>b.sign) then begin b.sign:=1; Subtract(a,b,c); b.sign:=-1; exit; end;
          if (b.sign>a.sign) then begin a.sign:=1; Subtract(b,a,c); a.sign:=-1; exit; end;
          if (a.sign=-1) and (b.sign=-1) then begin
              a.sign:=1; b.sign:=1;
              Add(a,b,c);
              a.sign:=-1; b.sign:=-1; c.sign:=-1;
              exit;
              end;
          {from here we can assume a and b both positive}
          c.sign:=1;
          c.decpoint:=max(a.decpoint,b.decpoint);
          if (a.decpoint=b.decpoint) and (AtoI(a.s[1])+AtoI(b.s[1])>=9) then inc(c.decpoint);
          carry:=0;
          ai:=digits+a.decpoint-c.decpoint;
          bi:=digits+b.decpoint-c.decpoint;
          for ci:=digits downto 1 do begin
            if (ai>a.sigdigits) then as:='x'
            else if (ai>digits) then begin
              if (a.sigdigits=exact) then as:='0' else as:='x';
              end
            else if (ai<1) then as:='0'
            else as:=a.s[ai];
            if (bi>b.sigdigits) then bs:='x'
            else if (bi>digits) then begin
              if (b.sigdigits=exact) then bs:='0' else bs:='x';
              end
            else if (bi<1) then bs:='0'
            else bs:=b.s[bi];
            AddDigits(as,bs,c.s[ci],carry);
            dec(ai); dec(bi);
            end;
          if (carry>0) then begin
            for ci:=digits downto 2 do c.s[ci]:=c.s[pred(ci)];
            c.s[1]:=char(48+carry);
            inc(c.decpoint);
            end;
          Normalize(c);
          if (a.sigdigits=exact) and (b.sigdigits=exact) then c.sigdigits:=exact
          else c.sigdigits:=c.decpoint+min(a.sigdigits-a.decpoint,b.sigdigits-b.decpoint);
          end;

procedure Subtract(a,b :number; var c :number);
          var i,stopi,
              ai,bi,ci   :integer;
              as,bs,r    :char;
              borrow     :integer;
          begin
          if (b.sign=0) then begin c:=a; exit; end;
          if (a.sign=0) then begin c:=b; c.sign:=-c.sign; exit; end;
          if (a.sign>b.sign) then begin b.sign:=1; Add(a,b,c); b.sign:=-1; exit; end;
          if (b.sign>a.sign) then begin a.sign:=1; Add(a,b,c); a.sign:=-1; c.sign:=-1; exit; end;
          if (Compare(a,b)=0) then begin c:=zero; exit; end;
          if (a.sign=-1) and (b.sign=-1) then begin
              a.sign:=1; b.sign:=1;
              Subtract(a,b,c);
              a.sign:=-1; b.sign:=-1; c.sign:=-c.sign;
              exit;
              end;
          case Compare(a,b) of
             0 : begin
                 c:=zero;
                 if (a.sigdigits=exact) and (b.sigdigits=exact) then exit;
                 c.sigdigits:=min(a.sigdigits,b.sigdigits)-a.decpoint;
                 end;
            -1 : begin
                 Subtract(b,a,c);
                 c.sign:=-c.sign;
                 exit;
                 end;
            end; {case}
          {from here we can assume a and b both positive, and a>b}
          c.sign:=1;
          c.decpoint:=max(a.decpoint,b.decpoint);
          borrow:=0;
          ai:=digits+a.decpoint-c.decpoint;
          bi:=digits+b.decpoint-c.decpoint;
          for ci:=digits downto 1 do begin
            if (ai>a.sigdigits) then as:='x'
            else if (ai>digits) then begin
              if (a.sigdigits=exact) then as:='0' else as:='x';
              end
            else if (ai<1) then as:='0'
            else as:=a.s[ai];
            if (bi>b.sigdigits) then bs:='x'
            else if (bi>digits) then begin
              if (b.sigdigits=exact) then bs:='0' else bs:='x';
              end
            else if (bi<1) then bs:='0'
            else bs:=b.s[bi];
            SubtractDigits(as,bs,c.s[ci],borrow);
            dec(ai); dec(bi);
            end;
          Normalize(c);
          if (a.sigdigits=exact) and (b.sigdigits=exact) then c.sigdigits:=exact
          else c.sigdigits:=c.decpoint+min(a.sigdigits-a.decpoint,b.sigdigits-b.decpoint);
          end;

procedure Multiply(a,b :number; var c :number);
          var part :number;
              i    :integer;
          begin
          if (a.sign=0) or (b.sign=0) then begin c:=zero; exit; end;
          if (a.sign>b.sign) then begin
              b.sign:=1; Multiply(a,b,c); b.sign:=-1;
              if (Compare(c,zero)=1) then c.sign:=-1; exit;
              end;
          if (b.sign>a.sign) then begin
              a.sign:=1; Multiply(a,b,c); a.sign:=-1;
              if (Compare(c,zero)=1) then c.sign:=-1; exit;
              end;
          if (a.sign=-1) and (b.sign=-1) then begin
              a.sign:=1; b.sign:=1;
              Multiply(a,b,c);
              a.sign:=-1; b.sign:=-1;
              exit;
              end;
          {from here we can assume a and b both positive}
          c:=zero;
          for i:=1 to digits do if (b.s[i]>'0') and (b.s[i]<='9') then begin
            NumberTimesDigit(a,b.s[i],part);
            part.decpoint:=part.decpoint+b.decpoint-i;
            Add(c,part,c);
            end;
          {if (c.s[1]='0') then begin
            for i:=1 to pred(digits) do c.s[i]:=c.s[succ(i)];
            c.s[digits]:='0';
            dec(c.decpoint);
            end;}
          c.sigdigits:=min(a.sigdigits,b.sigdigits);
          Normalize(c);
          end;

procedure FastMultiply(a,b :number; var c :number);
          var part                          :number;
              aa,bb,cc                      :longint;
              as,bs,cs                      :string;
              i,j,ii,jj,adigits,bdigits,err :integer;
          begin
          if (a.sign=0) or (b.sign=0) then begin c:=zero; exit; end;
          if (a.sign>b.sign) then begin
              b.sign:=1; FastMultiply(a,b,c); b.sign:=-1;
              if (Compare(c,zero)=1) then c.sign:=-1; exit;
              end;
          if (b.sign>a.sign) then begin
              a.sign:=1; FastMultiply(a,b,c); a.sign:=-1;
              if (Compare(c,zero)=1) then c.sign:=-1; exit;
              end;
          if (a.sign=-1) and (b.sign=-1) then begin
              a.sign:=1; b.sign:=1;
              FastMultiply(a,b,c);
              a.sign:=-1; b.sign:=-1;
              exit;
              end;
          {from here we can assume a and b both positive}
          adigits:=min(digits,a.sigdigits);
          bdigits:=min(digits,b.sigdigits);
          c:=zero;
          i:=1;
          repeat
            while (i<=adigits) and (a.s[i]='0') and (a.s[i]<='9') do inc(i);
            if (i<adigits) and (a.s[i] in digitset) then begin
              ii:=i;
              as:=a.s[i]; inc(i);
              if (i<adigits) and (a.s[i] in digitset) then begin
                as:=as+a.s[i]; inc(i);
                if (i<adigits) and (a.s[i] in digitset) then begin
                   as:=as+a.s[i]; inc(i);
                   end;
                end; {maximum of three digits}
              Val(as,aa,err); if (err<>0) then writeln('Numerical conversion error in Multiply: "',as,'"');
              j:=1;
              repeat
                while (j<=bdigits) and (b.s[j]='0') and (b.s[j]<='9') do inc(j);
                if (j<adigits) and (b.s[j] in digitset) then begin
                  jj:=j;
                  bs:=b.s[j]; inc(j);
                  if (j<bdigits) and (b.s[j] in digitset) then begin
                    bs:=bs+b.s[j]; inc(j);
                    if (j<bdigits) and (b.s[j] in digitset) then begin
                       bs:=bs+b.s[j]; inc(j);
                       end;
                    end; {maximum of three digits}
                  Val(bs,bb,err); if (err<>0) then writeln('Numerical conversion error in Multiply: "',as,'"');
                  cc:=aa*bb;
                  Str(cc,cs);
                  StringToNum(cs,part);
                  part.decpoint:=part.decpoint+a.decpoint+b.decpoint-ii-jj-4;
                  Add(c,part,c);
                  end;
              until (j>=bdigits) or (b.s[j]>'9');
              end;
          until (i>=adigits) or (a.s[i]>'9');
          c.sigdigits:=min(a.sigdigits,b.sigdigits);
          Normalize(c);
          end;

procedure IntegerPower(a,n :number; var c :number);
   {computes the nth power of a, where n is an integer; otherwise, rounds n to nearest integer.}
          var i,nn :integer;
          begin
          nn:=round(NumValue(n));
          c:=a;
          for i:=1 to nn-1 do
             FastMultiply(a,c,c);
          end;

procedure Factorial(a :number; var c :number);
          var i,n :integer;
              b   :number;
          begin
          n:=abs(round(NumValue(a)));
          if (n=0) then begin c:=one; exit; end;
          c:=a;
          for i:=2 to pred(n) do begin
            IntToNum(i,b);
            Multiply(c,b,c);
            end;
          if (n<0) then c.sign:=-1;
          end;

procedure Divide(a,b :number; var c :number);
          var aa,apart,cpart  :number;
              ra,rb,rpart     :double;
              ta,lastdecpoint :integer;
              done            :boolean;
              sigdigs         :integer;
          begin
          if (a.sign=0) then begin c:=zero; exit; end;
          if (b.sign=0) then begin writeln('Divide by zero error.'); c:=zero; exit; end;
          if (a.sign>b.sign) then begin b.sign:=1; Divide(a,b,c); b.sign:=-1; c.sign:=-1; exit; end;
          if (b.sign>a.sign) then begin a.sign:=1; Divide(a,b,c); a.sign:=-1; c.sign:=-1; exit; end;
          if (a.sign=-1) and (b.sign=-1) then begin
              a.sign:=1; b.sign:=1;
              Divide(a,b,c);
              a.sign:=-1; b.sign:=-1;
              exit;
              end;
          {from here we can assume a and b both positive}
          sigdigs:=min(a.sigdigits,b.sigdigits);
          sigdigs:=min(sigdigs,digits);
          aa:=a;
          c:=zero;
          rb:=NumValue(b);
          repeat
            ra:=NumValue(aa);
            rpart:=ra/rb;
            DoubleToNum(rpart,cpart); cpart.sigdigits:=exact;
            Add(c,cpart,c);
            FastMultiply(cpart,b,apart);
            lastdecpoint:=aa.decpoint;
            Subtract(aa,apart,aa);
            done:=(Compare(aa,zero)=0);
          until (done) or (aa.decpoint<=a.decpoint-sigdigs+3) or (lastdecpoint-aa.decpoint<=0);
          if (done) then c.sigdigits:=digits
          else begin
            c.sigdigits:=min(a.sigdigits,b.sigdigits);
            c.sigdigits:=min(c.sigdigits,digits);
            c.sigdigits:=min(c.sigdigits,a.decpoint-aa.decpoint);
            end;
          c.sign:=1;
          end;

procedure SqRoot(x :number; var y :number);
          var result,w,z,two :number;
              r            :double;
          begin
          IntToNum(2,two);
          r:=sqrt(NumValue(x));
          DoubleToNum(r,y);
          repeat
            Divide(x,y,z);
            Add(y,z,w);
            Divide(w,two,y);
            Multiply(y,y,w);
            Subtract(w,x,z);
          until (Equal(z,zero)) or (x.decpoint-z.decpoint>min(digits,x.sigdigits)-8);
          end;

{$ifdef STANDALONE}
(* This converges much too slowly...
procedure ComputePi; {uses arcsin(x)=x + 1/2 x^3 + 1/2 * 3/4 x^5/5 + 1/2 * 3/4 * 5/6 x^7/7 + ...}
          var piby2,pi,n,frac,term,wn :number;
              i,ni                  :integer;
          begin
          piby2:=one;
          ni:=1;
          frac:=one;
          for i:=1 to 100 do begin
             inc(ni,2);
             IntToNum(ni-2,wn);
             Multiply(frac,wn,frac);
             IntToNum(ni-1,wn);
             Divide(frac,wn,frac);
             IntToNum(ni,wn);
             Divide(frac,wn,term);
             Add(piby2,term,piby2);
             end;
          Add(piby2,piby2,pi);
          writeln('pi = ',NumString(pi));
          end;

procedure ComputePi; {uses arctan(x)=x - 1/3 x^3 + 1/5 * x^5 - 1/7 * x^7 + ... for x=1}
          var piby4,piby2,pi,n,term,wn :number;
              i,ni                  :integer;
          begin
          piby4:=one;
          ni:=1;
          for i:=1 to 1000 do begin
             inc(ni,2);
             IntToNum(ni,wn);
             Divide(one,wn,term);
             Subtract(piby4,term,piby4);
             inc(ni,2);
             IntToNum(ni,wn);
             Divide(one,wn,term);
             Add(piby4,term,piby4);
             end;
          Add(piby4,piby4,piby2);
          Add(piby2,piby2,pi);
          writeln('pi = ',NumString(pi));
          end;
*)
var fout:text;
procedure ComputePi; {uses arctan(x)=x - 1/3 x^3 + 1/5 * x^5 - 1/7 * x^7 + ... for x=û3/3, arctan=ã/6}
          var piby6,six,pi,n,frac,term,wn,x,x2,xn :number;
              i,ni                                :integer;
          begin
          IntToNum(3,wn); {xý=1/3}
          Divide(one,wn,x2);
          SqRoot(x2,x); {x=û3/3}
          piby6:=x;
          ni:=1;
          xn:=x;
          for i:=1 to 200 do begin
             inc(ni,2);
             IntToNum(ni,n);
             Multiply(xn,x2,xn);
             Divide(xn,n,term);
             Subtract(piby6,term,piby6);
             inc(ni,2);
             IntToNum(ni,n);
             Multiply(xn,x2,xn);
             Divide(xn,n,term);
             Add(piby6,term,piby6);
             end;
          IntToNum(6,six);
          Multiply(piby6,six,pi);
          writeln('pi = ',NumString(pi));
          end;

procedure ComputePiMoreAccurately; {uses arctan(x)=x - x^3/3 + x^5/5 -x^7/7 + ... for x=û((2-û3)/(2+û3)), arctan=ã/12}
          var piby12,twelve,pi,n,frac,term,
              cos,wn,x,x2,xn,num,denom           :number;
              i,ni                               :integer;
          begin
          StringToNum('.75',wn); {xý=3/4}
          SqRoot(wn,x); {x=û3/2}
          Subtract(one,x,num);
          Add(one,x,denom);
          Divide(num,denom,x2); {x= tan(ã/12}
          SqRoot(x2,x);
          piby12:=x;
          ni:=1;
          xn:=x;
          for i:=1 to 100 do begin
             inc(ni,2);
             IntToNum(ni,n);
             Multiply(xn,x2,xn);
             Divide(xn,n,term);
             Subtract(piby12,term,piby12);
             inc(ni,2);
             IntToNum(ni,n);
             Multiply(xn,x2,xn);
             Divide(xn,n,term);
             Add(piby12,term,piby12);
             end;
          IntToNum(12,twelve);
          Multiply(piby12,twelve,pi);
          writeln('pi = ',NumString(pi));
          end;

procedure FastPiComputation; {uses pi= (16/5 - 4/239) - (1/3) (16/5^3 - 4/239^3) + (1/5) (16/5^5 - 4/239^5) - ...}
          var a,b,a2,b2,num,denom,pi,knum,term :number;
              i,k                               :integer;
          begin
          IntToNum(16,num); IntToNum(5,denom); {a=16/5}
          Divide(num,denom,a);
          Multiply(denom,denom,denom); Divide(one,denom,a2); {a2=1/5ý}
          IntToNum(4,num); IntToNum(239,denom); {b=1/239}
          Divide(num,denom,b);
          Multiply(denom,denom,denom); Divide(one,denom,b2); {b2=1/239ý}
          Subtract(a,b,pi);
          k:=1;
          for i:=1 to 100 do begin
             inc(k,2); IntToNum(k,knum);
             Multiply(a,a2,a); Multiply(b,b2,b);
             Subtract(a,b,term);
             Divide(term,knum,term);
             Subtract(pi,term,pi);
             inc(k,2); IntToNum(k,knum);
             Multiply(a,a2,a); Multiply(b,b2,b);
             Subtract(a,b,term);
             Divide(term,knum,term);
             Add(pi,term,pi);
             end;
          writeln('pi = ',NumString(pi));
          end;

procedure ComputeE;
          var e,n,nfact,term :number;
              i              :integer;
              fout           :text;
          begin
          StringToNum('2.5',e);
          StringToNum('2',nfact);
          for i:=3 to 30 do begin
            IntToNum(i,n);
            Multiply(nfact,n,nfact);
            Divide(one,nfact,term);
            Add(e,term,e)
            end;
          writeln; writeln('e = ',NumString(e));
          end;

procedure Test;
          var a,b,c :number;
          begin
          writeln; writeln('DEMONSTRATION');
          {StringToNum('-7.3E-2',a);
          writeln(NumString(a),' = ',NumValue(a):8:6);
          StringToNum('-0.023E4',a);
          writeln(NumString(a),' = ',NumValue(a):8:6);
          StringToNum('123.45',a);
          writeln(NumString(a),' = ',NumValue(a):8:6);}
          StringToNum('123',a);
          StringToNum('8.75',b);
          Add(a,b,c);
          writeln(NumString(a),' + ',NumString(b),' = ',NumString(c));
          StringToNum('.123',a);
          StringToNum('.895',b);
          Add(a,b,c);
          writeln(NumString(a),' + ',NumString(b),' = ',NumString(c));
          Subtract(c,b,a);
          writeln(NumString(c),' - ',NumString(b),' = ',NumString(a));
          StringToNum('793.45',a);
          StringToNum('9.288',b);
          Add(a,b,c);
          writeln(NumString(a),' + ',NumString(b),' = ',NumString(c));
          Subtract(c,b,a);
          writeln(NumString(c),' - ',NumString(b),' = ',NumString(a));
          StringToNum('-923.66',a);
          StringToNum('875.28',b);
          Add(a,b,c);
          writeln(NumString(a),' + ',NumString(b),' = ',NumString(c));
          StringToNum('1.2345',a);
          StringToNum('1.2346',b);
          Subtract(a,b,c);
          writeln(NumString(a),' - ',NumString(b),' = ',NumString(c));
          StringToNum('-.003',a);
          StringToNum('1.234',b);
          Add(a,b,c);
          writeln(NumString(a),' + ',NumString(b),' = ',NumString(c));
          StringToNum('2',a);
          StringToNum('3',b);
          Multiply(a,b,c);
          writeln(NumString(a),' * ',NumString(b),' = ',NumString(c));
          Divide(c,a,b);
          writeln(NumString(c),' / ',NumString(a),' = ',NumString(b));
          StringToNum('.5',a);
          StringToNum('3',b);
          Divide(a,b,c);
          writeln(NumString(a),' / ',NumString(b),' = ',NumString(c));
          StringToNum('720',b);
          Divide(one,b,c);
          writeln(NumString(one),' / ',NumString(b),' = ',NumString(c));
          StringToNum('1.2345',a);
          StringToNum('-7.89',b);
          Multiply(a,b,c);
          writeln(NumString(a),' * ',NumString(b),' = ',NumString(c));
          StringToNum('6',b);
          Divide(one,b,c);
          writeln(NumString(one),' / ',NumString(b),' = ',NumString(c));
          StringToNum('144',a);
          StringToNum('-36',b);
          Divide(a,b,c);
          writeln(NumString(a),' / ',NumString(b),' = ',NumString(c));
          StringToNum('1',a);
          StringToNum('4',b);
          Divide(a,b,c);
          writeln(NumString(a),' / ',NumString(b),' = ',NumString(c));
          StringToNum('1',a);
          StringToNum('3',b);
          Divide(a,b,c);
          writeln(NumString(a),' / ',NumString(b),' = ',NumString(c));
          StringToNum('2',a);
          SqRoot(a,b);
          writeln('û',NumString(a),' = ',NumString(b));
          Multiply(b,b,c);
          Subtract(c,a,c);
          writeln('(û2)ý - 2 = ',NumString(c));
          StringToNum('1',a);
          StringToNum('89',b);
          Divide(a,b,c);
          writeln('1/89=',NumString(c));
          end;

procedure TestFastMultiply;
          var a,b,c,cc :number;
              i        :integer;
          begin
          writeln; writeln('DEMONSTRATION');
          StringToNum('123456789',a);
          StringToNum('9.87654321E-4',b);
          writeln(123456789*9.87654321E-4:17:12);
          writeln('Slow:');
          for i:=1 to 10000 do Multiply(a,b,c);
          writeln(NumString(c));
          writeln('Fast:');
          for i:=1 to 10000 do FastMultiply(a,b,cc);
          writeln(NumString(cc));
          readln;
          end;

procedure TestIntegerPower;
          var a,b,c :number;
          begin
          StringtoNum('2',a); StringtoNum('100',b);
          IntegerPower(a,b,c);
          writeln('2 to the 100th power = ',NumString(c));
          end;

procedure Temp;
          var a,b,c,nine,two,p :number;
              i       :integer;
              s       :string;
          begin
          StringtoNum('3',p); StringtoNum('9',nine);
          StringtoNum('2',two);
          for i:=1 to 100 do begin
            Multiply(p,nine,p);
            Subtract(p,one,a);
            Divide(a,two,b);
            SqRoot(b,c);
            writeln(i:2,' ',NumString(c));
            end;
          end;

{$endif STANDALONE}

begin
StringToNum('0',zero);
zero.sign:=0; {need this, else comparing with zero produces errors}
StringToNum('1',one);
{$ifdef STANDALONE}
{TestFastMultiply;}
{Temp;}
Test;
ComputeE;
FastPiComputation;
ComputePiMoreAccurately;
ComputePi;
TestIntegerPower;
{$endif}
end.