{$mode objfpc}
{$MODESWITCH NESTEDCOMMENTS+}
{$warn 6058 off}
{$warn 5025 off}
program speed_test;
uses sysutils
, strutils
, strings
, math
, Multi_Int
;
const prec = 50;
const n = 50;
type Rational = record
num : Multi_Int_XV;
den : Multi_Int_XV;
end;
var
big_int_size,
start_time,
end_time :int32;
delta :double;
b: array[0..n] of Rational;
c: array[0..n] of Rational;
one, tmp : Rational;
i, j : longint;
function GCD(const n : Multi_Int_XV; const m : Multi_Int_XV) : Multi_Int_XV;
var
a, b : Multi_Int_XV;
begin
a:=abs(n);
b:=abs(m);
While ((a <> 0) And (b <> 0)) do
If a > b Then a := a Mod b Else b := b Mod a;
GCD := a + b;
End;
function lcm(const a : Multi_Int_XV; const b : Multi_Int_XV) : Multi_Int_XV;
begin
lcm := (a * b) div gcd(a, b);
end;
procedure RationalSimplify( var r : Rational);
var
divisor : Multi_Int_XV;
begin
divisor := GCD(r.num, r.den);
r.num := r.num div divisor;
r.den := r.den div divisor;
if (r.den < 0) then
begin
r.num := -r.num;
r.den := -r.den;
end;
end;
function RationalAdd(const r1 : Rational; const r2 : Rational) : Rational;
var
res : Rational;
begin
res.num := r1.num * r2.den + r2.num * r1.den;
res.den := r1.den * r2.den;
RationalSimplify(res);
RationalAdd := res;
end;
function RationalSubtract(const r1 : Rational; const r2 : Rational) : Rational;
var
res : Rational;
begin
res.num := r1.num * r2.den - r2.num * r1.den;
res.den := r1.den * r2.den;
RationalSimplify(res);
RationalSubtract := res;
end;
function RationalMultiply(const r1 : Rational; const r2 : Rational) : Rational;
var
res : Rational;
begin
res.num := r1.num * r2.num;
res.den := r1.den * r2.den;
RationalSimplify(res);
RationalMultiply := res;
end;
function RationalDivide(const r1 : Rational; const r2 : Rational) : Rational;
var
res : Rational;
begin
res.num := r1.num * r2.den;
res.den := r1.den * r2.num;
RationalSimplify(res);
RationalDivide := res;
end;
function RationalDivide(const r1 : Rational; const r2 : LongInt) : Rational;
var
res : Rational;
begin
res.num := r1.num;
res.den := r1.den * r2;
RationalSimplify(res);
RationalDivide := res;
end;
operator := (const r : ansistring) z : Rational;
var
lhs, rhs : ansistring;
c, i, n : longint;
num, den : Multi_Int_XV;
begin
c:=0;
n:=length(r);
for i:=1 to n do
begin
if r[i]='/' then
begin
c:=i;
break;
end;
end;
lhs:=leftstr(r, c-1);
rhs:=rightstr(r, n-c);
z.num:=lhs;
z.den:=rhs;
RationalSimplify(z);
end;
operator := (const r : longint) z : Rational;
begin
z.num:=r;
z.den:=1;
end;
operator + (const x : Rational; const y : Rational) z : Rational;
begin
z:=RationalAdd(x, y);
end;
operator - (const x : Rational) z : Rational;
begin
z:=x;
z.num*=(-1);
end;
operator - (const x : Rational; const y : Rational) z : Rational;
begin
z:=RationalSubtract(x, y);
end;
operator * (const x : Rational; const y : Rational) z : Rational;
begin
z:=RationalMultiply(x, y);
end;
operator * (const x : Rational; const y : LongInt) z : Rational;
begin
z:=x;
z.num*=y;
RationalSimplify(z);
end;
operator / (const x : Rational; const y : Rational) z : Rational;
begin
z:=RationalDivide(x, y);
end;
operator / (const x : Rational; const y : LongInt) z : Rational;
begin
z:=RationalDivide(x, y);
end;
operator < (const a,b:Rational):Boolean;
var
n, m : Multi_Int_XV;
begin
n:=a.num*b.den;
m:=b.num*a.den;
result:=(n<m);
end;
operator <= (const a,b:Rational):Boolean;
var
n, m : Multi_Int_XV;
begin
n:=a.num*b.den;
m:=b.num*a.den;
result:=(n<=m);
end;
operator > (const a,b:Rational):Boolean;
var
n, m : Multi_Int_XV;
begin
n:=a.num*b.den;
m:=b.num*a.den;
result:=(n>m);
end;
operator >= (const a,b:Rational):Boolean;
var
n, m : Multi_Int_XV;
begin
n:=a.num*b.den;
m:=b.num*a.den;
result:=(n>=m);
end;
operator = (const a,b:Rational):Boolean;
var
n, m : Multi_Int_XV;
begin
n:=a.num*b.den;
m:=b.num*a.den;
result:=(n=m);
end;
function abs(const x : Rational): Rational;
var
a : Rational;
begin
a:=x;
if a.num<0 then a.num*=-1;
abs:=a;
end;
function RationalToStr(const n:Rational): ansistring;
var
s:ansistring;
begin
s:='';
if n.num>0 then s:=' ';
s:=s+n.num.toStr;
if n.den>1 then s:=s+'/'+n.den.toStr;
RationalToStr:=s;
end;
begin
big_int_size:= ceil( prec / (32*0.3010299956639812));
Multi_Int_Set_XV_Limit(big_int_size);
Multi_Int_Initialisation(big_int_size);
//---------------------------------------------------------
// Akiyama–Tanigawa algorithm for Bn
//---------------------------------------------------------
start_time:= GetTickCount64;
one := 1;
For i := 0 To n do
begin
c[i] := one / (i+1);
For j := i downTo 1 do
begin
tmp:=c[j - 1]-c[j];
c[j - 1] := tmp * j;
end;
b[i]:=c[0];
end;
end_time:= GetTickCount64;
b[1]:=-b[1]; // the algorithm is wrong, Bernoulli[1] = -1/2
For i := 0 To n do
begin
if i>2 then
begin
if (i and 1)=0 then
begin
writeln('B[',i:2,'] = ', RationalToStr(b[i]));
end;
end
else
begin
writeln('B[',i:2,'] = ', RationalToStr(b[i]));
end;
end;
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');
end.
Recent
