program mp_test;
{$mode objfpc}
{$modeswitch advancedrecords}
{$warn 5089 off} //suppress uninitialized warning
{$warn 5090 off} //suppress uninitialized warning
{$warn 5092 off} //suppress uninitialized warning
{$warn 5025 off} //suppress variable not used warning
{$warn 4055 off} //suppress not portable
uses
mp_types, btypes, mp_real, mp_base, mp_numth, mp_calc, sysutils, windows, math, strings, JwaWinBase;
const default_prec=125;
const print_prec=80;
type
mpfloat = record
public
mpf: mp_float;
strict private
class operator Initialize(var aRec: mpfloat);
class operator Finalize(var aRec: mpfloat);
class operator Copy(constref cFrom:mpfloat; var cTo:mpfloat);
public
class operator < (const a,b:mpfloat):Boolean;
class operator <= (const a,b:mpfloat):Boolean;
class operator > (const a,b:mpfloat):Boolean;
class operator >= (const a,b:mpfloat):Boolean;
class operator <> (const a,b:mpfloat):Boolean;
class operator = (const a,b:mpfloat):Boolean;
end;
class operator mpfloat.Initialize(var aRec: mpfloat);
begin
ReturnNilIfGrowHeapFails := True;
mpf_set_default_prec(round(default_prec*3.321928094887362));
mpf_init(aRec.mpf);
end;
class operator mpfloat.Finalize(var aRec: mpfloat);
begin
mpf_clear(aRec.mpf);
end;
class operator mpfloat.Copy(constref cFrom:mpfloat;var cTo:mpfloat);
begin
mpf_clear(cTo.mpf);
mpf_set_default_prec(round(default_prec*3.321928094887362));
mpf_init(cTo.mpf);
mpf_copy(cFrom.mpf, cTo.mpf);
end;
//----------------------------------------------------------------
procedure writelnq(const n:mpfloat; const digits : longint=80);
begin
if s_mpf_is_neg(n.mpf) then
begin
mpf_output_decimal(n.mpf, digits);
writeln;
end
else
begin
write(' ');
mpf_output_decimal(n.mpf, digits);
writeln;
end;
end;
procedure writeq(const n:mpfloat; const digits : longint=80);
begin
if s_mpf_is_neg(n.mpf) then
mpf_output_decimal(n.mpf, digits)
else
begin
write(' ');
mpf_output_decimal(n.mpf, digits);
end;
end;
operator := (const r : pchar8) z : mpfloat;
begin
mpf_read_decimal(z.mpf, r);
end;
operator := (const r : longint) z : mpfloat;
begin
mpf_set_int(z.mpf, r);
end;
operator := (const r : double) z : mpfloat;
begin
mpf_set_dbl(z.mpf, r);
end;
operator + (const x : mpfloat; const y : mpfloat) z : mpfloat;
begin
mpf_add(x.mpf, y.mpf, z.mpf);
end;
operator - (const x : mpfloat; const y : mpfloat) z : mpfloat;
begin
mpf_sub(x.mpf, y.mpf, z.mpf);
end;
operator * (const x : mpfloat; const y : mpfloat) z : mpfloat;
begin
mpf_mul(x.mpf, y.mpf, z.mpf);
end;
operator * (const x : mpfloat; const y : longint) z : mpfloat;
begin
mpf_mul_int(x.mpf, y, z.mpf);
end;
operator * (const x : longint; const y : mpfloat ) z : mpfloat;
var r:mpfloat;
begin
mpf_set_int(r.mpf, x);
mpf_mul(r.mpf, y.mpf, z.mpf);
end;
operator / (const x : mpfloat; const y : mpfloat) z : mpfloat;
begin
mpf_div(x.mpf, y.mpf, z.mpf);
end;
operator / (const x : mpfloat; const y : longint) z : mpfloat;
begin
mpf_div_int(x.mpf, y, z.mpf);
end;
operator / (const x : longint; const y : mpfloat ) z : mpfloat;
var r:mpfloat;
begin
mpf_set_int(r.mpf, x);
mpf_div(r.mpf, y.mpf, z.mpf);
end;
class operator mpfloat.< (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)<0);
end;
class operator mpfloat.<= (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)<=0);
end;
class operator mpfloat.> (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)>0);
end;
class operator mpfloat.>= (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)>=0);
end;
class operator mpfloat.<> (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)<>0);
end;
class operator mpfloat.= (const a,b:mpfloat):Boolean;
begin
result:=(mpf_cmp(a.mpf, b.mpf)=0);
end;
function sqrt(const x : mpfloat):mpfloat;
begin
mpf_sqrt(x.mpf, sqrt.mpf);
end;
function sin(const x : mpfloat): mpfloat;
begin
mpf_sin(x.mpf, sin.mpf);
end;
function cos(const x : mpfloat): mpfloat;
begin
mpf_cos(x.mpf, cos.mpf);
end;
function tan(const x : mpfloat): mpfloat;
begin
mpf_tan(x.mpf, tan.mpf);
end;
function arcsin(const x : mpfloat): mpfloat;
begin
mpf_arcsin(x.mpf, arcsin.mpf);
end;
function arccos(const x : mpfloat): mpfloat;
begin
mpf_arccos(x.mpf, arccos.mpf);
end;
function arctan(const x : mpfloat): mpfloat;
begin
mpf_arctan(x.mpf, arctan.mpf);
end;
function log(const x : mpfloat): mpfloat;
begin
mpf_ln(x.mpf, log.mpf);
end;
function exp(const x : mpfloat): mpfloat;
begin
mpf_exp(x.mpf, exp.mpf);
end;
function abs(const x : mpfloat): mpfloat;
begin
mpf_abs(x.mpf, abs.mpf);
end;
//=============================================================
const n = 50;
var
start_time,
end_time :int32;
delta, t :double;
b: array[0..n] of mpfloat;
c: array[0..n] of mpfloat;
one, tmp : mpfloat;
i, j : longint;
begin
//---------------------------------------------------------
// Akiyama–Tanigawa algorithm for Bn
//---------------------------------------------------------
t:=round((n+print_prec)*0.96);
if (round(t/8)*8)<t then
t:=(round(t/8)+1)*8;
if default_prec<t then
begin
writeln('for ',n,' bernoulli numbers at ',print_prec,' digits of precision you need ',t,' digits of precision');
end
else
begin
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;
For i := 0 To n do
begin
if i>2 then
begin
if (i and 1)=0 then
begin
write('B[',i:2,'] = '); writelnq(b[i], print_prec);
end;
end
else
begin
write('B[',i:2,'] = '); writelnq(b[i], print_prec);
end;
end;
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');
end;
end.
Recent
