双曲線正割関数のテイラー展開の展開係数と定義されている、Ｅｋがオイラー数 詳しいことは、インターネットで検索してください。
オイラー数 その１
ベルヌーイ数の計算を取り上げたので、オイラー数の計算も追加しました。
双曲線正割関数のテイラー展開の展開係数と定義されている、Ｅｋがオイラー数
詳しいことは、インターネットで検索してください。
とにかく、オイラー数の計算式をインターネットで検索してプログラムを組んでみることにしました。
プログラムは、Bigintegerで計算をしています。
Bigintegerの組み込みはベルヌーイ数 その４を参照して下さい。
左図は、７種類の計算式を取り上げてみたプログラムの実行画面です。
オイラー数の計算を実行させると、選んだボタンの計算式が左下に表示されるので参考にしてください。
オイラー数W～5迄の計算は非常に計算が遅いのでオイラー数Noを200に制限しています。
計算式の中に二項係数の計算があり、この計算の繰り返しにより、計算がおそくなっています。
階乗の計算は、テーブルを用意して、計算を省略していますが、整数割り算、DIVによる影響です。
二項係数の高速演算をインターネットで検索すると、Long Long
(64ビット)整数で、素数を利用して割り算を避けるプログラムが紹介されていて、ｎ，ｋの値が結構大きく設定できるようになってますが、正しい値を返すのは、32迄で、それを超えるとLong
Long (64ビット)の値を超えてしまいます。
２０００位迄計算しようとすると、６００桁を超えるので、Bigintegerの計算が必要となります。
これほど大きくなると、必要な素数を割り出すだけで、非常に時間がかかるので、二項係数の高速演算は使用していません。
７種類の計算の中で一番計算が速かったのは、オイラー数６の計算で、他のものと計算式を比べると一番単純です。
式は再帰を利用するような形になっていますが、ここでは For Next .ループで計算しています。
ｎ＝1000の計算で、約１秒程の計算時間となり、２０００にすると十数秒の時間となりますが、ベルヌーイ数の計算がｎ＝２０００で１秒で計算できるので、可なり時間が掛かっています。
プログラムは、1000迄の計算に制限していますが、もっと大きな値を指定する場合は、const の Nk =
1000 の値を変更します。
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.Buttons, system.Diagnostics,
Vcl.ExtCtrls, Vcl.Imaging.pngimage;
type
TForm1 = class(TForm)
Memo1: TMemo;
BigInteger3Btn: TBitBtn;
BigIntegerwBtn: TBitBtn;
BigInteger4Btn: TBitBtn;
BigInteger1Btn: TBitBtn;
BigInteger2Btn: TBitBtn;
Image1: TImage;
BigInteger5Btn: TBitBtn;
LabeledEdit1: TLabeledEdit;
BigInteger6Btn: TBitBtn;
LabeledEdit2: TLabeledEdit;
procedure BigInteger3BtnClick(Sender: TObject);
procedure BigIntegerwBtnClick(Sender: TObject);
procedure BigInteger4BtnClick(Sender: TObject);
procedure BigInteger1BtnClick(Sender: TObject);
procedure BigInteger2BtnClick(Sender: TObject);
procedure BigInteger5BtnClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure BigInteger6BtnClick(Sender: TObject);
private
{ Private 宣言 }
procedure loadimage(n: integer);
function Ninput(var n: integer): boolean;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses Velthuis.BigIntegers, Velthuis.BigDecimals;
const
Nj = 200; // 最大オイラー w ～ 5
Nk = 1000; // 最大オイラー 6
var
fact: array[0..Nk] of BigInteger; // 階乗テーブル配列
// 数式表示
procedure TForm1.loadimage(n: integer);
begin
case n of
0: Image1.Picture.LoadFromFile('Wikipedia.png');
1: Image1.Picture.LoadFromFile('euler1.png');
2: Image1.Picture.LoadFromFile('euler2.png');
3: Image1.Picture.LoadFromFile('euler3.png');
4: Image1.Picture.LoadFromFile('euler4.png');
5: Image1.Picture.LoadFromFile('euler5.png');
6: Image1.Picture.LoadFromFile('euler6.png');
end;
end;
// オイラーNo入力処理
function TForm1.Ninput(var n: integer): boolean;
var
c : integer;
begin
result := true;
val(labelededit1.Text, n, c);
if c <> 0 then begin
memo1.Text := 'オイラーNo に間違いがあります。';
result := false;
exit;
end;
if n > Nj then begin
memo1.Text := 'オイラーNo が大きすぎます' + intTostr(Nj) + '迄です。';
result := false;
end;
end;
// 二項係数 BigInteger
function binomialBig(n, k: integer): BigInteger;
var
nn, kn, knm : BigInteger;
begin
nn := fact[n];
kn := fact[k];
knm := fact[n - k];
result := nn div (kn * knm);
end;
// オイラー数 BigInteger
procedure TForm1.BigIntegerwBtnClick(Sender: TObject);
var
n : integer;
i, l , k : integer;
s, sf, mone, two, il : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
memo1.Clear;
if not Ninput(n) then exit;;
StopWatch := TStopwatch.StartNew;
mone := BigInteger.Negate(BigInteger.One); // -1
two := 2; // 2
for k := 0 to n div 2 do begin
s := 0;
for i := 1 to 2 * k do begin // 偶数のみ計算
sf := 0;
for l := 0 to 2 * i do begin
il := i - l;
sf := sf + BigInteger.Pow(mone, l) * binomialBig(2 * i, l) * BigInteger.pow(il, 2 * k);
end;
s := s + sf div BigInteger.Pow(two, i) * BigInteger.Pow(mone, i);
end;
if k = 0 then s := 1;
if s < 0 then mstr := ' = '
else mstr := ' = ';
if (n = k * 2) or (k <= 20) then begin
if (n = k * 2) and (k > 20) then Memo1.Lines.Append('');
Memo1.Lines.Append('E' + intTostr(2 * k) + mstr + s.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(0);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger1BtnClick(Sender: TObject);
var
i : integer;
n, k, j : integer;
s, sf, mone, two, jj : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not Ninput(i) then exit;;
memo1.Clear;
StopWatch := TStopwatch.StartNew;
mone := BigInteger.Negate(BigInteger.One); // -1
two := 2; // 2
for n := 0 to i div 2 do begin
s := 0;
for k := 1 to n do begin // 偶数のみ計算
sf := 0;
for j := 1 to k do begin
jj := j;
sf := sf + BigInteger.Pow(mone, j) * binomialBig(2 * k, k - j) * BigInteger.Pow(jj, 2 * n)
end;
s := s + sf div BigInteger.Pow(two, k - 1);
end;
if n = 0 then s := 1;
if s < 0 then mstr := ' = '
else mstr := ' = ';
if (i = n * 2) or (n <= 20) then begin
if (i = n * 2) and (n > 20) then Memo1.Lines.Append('');
Memo1.Lines.Append('E' + intTostr(n * 2) + mstr + s.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(1);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger2BtnClick(Sender: TObject);
var
i : integer;
n, k, j : integer;
s, sf, mone, two, nmk : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not Ninput(i) then exit;;
memo1.Clear;
StopWatch := TStopwatch.StartNew;
mone := BigInteger.Negate(BigInteger.One); // -1
two := 2; // 2
for n := 0 to i div 2 do begin
s := 0;
for k := 0 to n - 1 do begin // 偶数のみ計算
sf := 0;
for j := 0 to k do begin
sf := sf + binomialBig(2 * n - 2 * j, k - j) * BigInteger.Pow(two, j);
end;
nmk := n - k;
s := s + sf * BigInteger.Pow(mone, n - k) * BigInteger.Pow(nmk, 2 * n);
end;
if n = 0 then s := 1
else s := s div BigInteger.Pow(two, n - 1);
if s < 0 then mstr := ' = '
else mstr := ' = ';
if (i = n * 2) or (n <= 20) then begin
if (i = n * 2) and (n > 20) then Memo1.Lines.Append('');
Memo1.Lines.Append('E' + intTostr(n * 2) + mstr + s.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(2);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger3BtnClick(Sender: TObject);
var
n : integer;
s, sf, mone, two : BigInteger;
J2P1 : BigInteger;
j, k, i, l : integer;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not Ninput(n) then exit;;
memo1.clear;
StopWatch := TStopwatch.StartNew;
mone := BigInteger.Negate(BigInteger.One); // -1
two := 2; // 2
for l := 0 to n div 2 do begin
i := l * 2; // 偶数のみ計算
s := BigInteger.Zero; // s = 0
for k := 1 to i + 1 do begin
sf := BigInteger.Zero; // sf = 0
for j := 0 to k - 1 do begin
J2P1:= 2 * j + 1;
sf := sf + BigInteger.Pow(mone, j) * BigInteger.Pow(J2P1, i);
end;
s := s + sf * binomialBig(i + 1, k);
end;
s := s div BigInteger.Pow(two, i);
if s < 0 then mstr := ' = '
else mstr := ' = ';
if (i = n) or (i <= 40) then begin
if (i = n) and (i > 40) then Memo1.Lines.Append('');
Memo1.Lines.Append('E' + intTostr(i) + mstr + s.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(3);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger4BtnClick(Sender: TObject);
var
i : integer;
n, k, j : integer;
s, sf, mone, two, j2p1 : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not Ninput(i) then exit;;
memo1.Clear;
StopWatch := TStopwatch.StartNew;
mone := BigInteger.Negate(BigInteger.One); // -1
two := 2; // 2
for n := 0 to i div 2 do begin
s := 0;
for k := 0 to 2 * n do begin // 偶数のみ計算
sf := 0;
for j := 0 to k do begin
j2p1 := (2 * j + 1);
sf := sf + BigInteger.Pow(mone, j) * binomialBig(k, j) * BigInteger.Pow(j2p1, 2 * n);
end;
s := s + sf div BigInteger.Pow(two, k);
end;
if s < 0 then mstr := ' = '
else mstr := ' = ';
if (n = i div 2) or (n <= 20) then begin
if (n = i div 2) and (n > 20) then Memo1.Lines.Append('');
Memo1.Lines.Append('E' + intTostr(n * 2) + mstr + s.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(4);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger5BtnClick(Sender: TObject);
var
N :integer;
En : array of BigInteger;
i, k : integer;
Sn, Sd : BigInteger; // 分子,分母
Tn, Td : BigInteger; // 分子,分母
two, one : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not Ninput(n) then exit;;
memo1.Clear;
StopWatch := TStopwatch.StartNew;
setlength(En, N + 1);
one := BigInteger.One;
two := 2;
En[0] := 1;
memo1.Lines.Append('E' + intTostr(0) + ' = ' + En[0].ToString);
for i := 1 to N do begin
Sn := 0; // 分子初期値
Sd := 1; // 分母初期値
for k := 0 to i - 1 do begin
Tn := En[k] * binomialBig(i, k); // 分子 計算
Td := BigInteger.Pow(two, k); // 分母 計算
Sn := Sn * Td + Tn * Sd; // 分子 分数加算
Sd := Sd * Td; // 分母
end;
Sn := Sn * BigInteger.Pow(two, i - 1) div Sd; // 分母は1になるので分子のみ計算
En[i] := one - Sn;
if i mod 2 = 0 then begin
if En[i] < 0 then mstr := ' = '
else mstr := ' = ';
if (N = i) or (i <= 40) then begin
if (N = i) and (i > 40) then Memo1.Lines.Append('');
memo1.Lines.Append('E' + intTostr(i) + mstr + En[i].ToString);
end;
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(5);
end;
// オイラー数 BigInteger
procedure TForm1.BigInteger6BtnClick(Sender: TObject);
var
n, n2, nd2 : integer;
En : array of BigInteger;
i, k : integer;
s : BigInteger;
mstr : string;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
val(labelededit2.Text, n, i);
if i <> 0 then begin
memo1.Text := 'オイラーNo に間違いがあります。';
exit;
end;
if n > Nk then begin
memo1.Text := 'オイラーNo が大きすぎます' + intTostr(Nk) + '迄です。';
exit;
end;
memo1.Clear;
StopWatch := TStopwatch.StartNew;
nd2 := n div 2;
setlength(En, nd2 + 1); // No=0 + 偶数No用配列確保
i := 0; // Kn [0] [1] [2] [3] [4] [5] [6]・・
En[i] := 1; // En 0 2 4 6 8 10 12・・
memo1.Lines.Append('E' + intTostr(i) + ' = ' + En[i].ToString);
for i := 1 to nd2 do begin
s := 0;
n2 := i shl 1; // n2 = i * 2
for k := 0 to i - 1 do
s := s + binomialBig(n2, k shl 1) * En[k]; // Σ
En[i] := -s; // En[No/2] = -s
if (nd2 = i) or (i <= 20) then begin
if s > 0 then mstr := ' = '
else mstr := ' = ';
if (nd2 = i) and (i > 20) then Memo1.Lines.Append('');
memo1.Lines.Append('E' + intTostr(n2) + mstr + En[i].ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Memo1.text := Memo1.text + '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
loadimage(6);
end;
// 階乗テーブル作成
procedure TForm1.FormCreate(Sender: TObject);
var
i : integer;
fr : BigInteger;
begin
// w～6 用 階乗テーブル
fr := 1;
fact[0] := fr;
fact[1] := fr;
for i := 2 to Nk do begin // Nk = 1000
fr := fr * i;
fact[i] := fr;
end;
labeledEdit1.EditLabel.Caption := 'En n<=' + intTostr(Nj);
labeledEdit2.EditLabel.Caption := 'En n<=' + intTostr(Nk);
end;
end.
Euler_number_2_biginteger.zip
各種プログラム計算例に戻る