ディリクレのイータ関数
イータ関数は、自由電子気体の電子比熱の計算に使用されています。
プログラムは、Mathematics Source Library C & ASM のディレクリ・イータ関数と、C言語辞典のリーマン・ゼータ関数を元にしたディリクレ・イータ関数です。
プログラム
次のリストは Mathematics Source Library C & ASM のイータ関数をDelphiに変換したものですが、Mathematics Source Library C & ASM のゼータ関数はイータ関数を元に計算しています。
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, system.Math, Vcl.StdCtrls, Vcl.Buttons,
Vcl.ExtCtrls, VclTee.TeeGDIPlus, VCLTee.Series, VCLTee.TeEngine,
VCLTee.TeeProcs, VCLTee.Chart;
type
TForm1 = class(TForm)
BitBtn1: TBitBtn;
LabeledEdit1: TLabeledEdit;
Memo1: TMemo;
Chart1: TChart;
Series1: TLineSeries;
Series2: TPointSeries;
procedure FormCreate(Sender: TObject);
procedure BitBtn1Click(Sender: TObject);
private
{ Private 宣言 }
function Dirichlet_Eta_Function(s: double): double;
function xDirichlet_Eta_Function(s: extended): extended;
function Dirichlet_Eta_Star_Function(s: double): double;
function xDirichlet_Eta_Star_Function(s : extended): extended;
function Reflection_Coefficient(s: extended): extended;
function Sum_Reverse_Order(s: extended): extended;
function Alternating_Series_Convergence_Acceleration(s: extended): extended;
function Duplication_Formula(two_x: extended): extended;
function xGamma(x: extended): extended;
function xGamma_Function(x: extended): extended;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
const
LONG_MAX = 2147483647;
max_long_double_arg: extended = 1755.5;
exp_g_o_sqrt_2pi = 6.23316569877722552586386e+3;
e = 2.71828182845904523536028747;
g = 9.65657815377331589457187;
a: array[0..8] of extended = (
+1.14400529453851095667309e+4,
-3.23988020152318335053598e+4,
+3.50514523505571666566083e+4,
-1.81641309541260702610647e+4,
+4.63232990536666818409138e+3,
-5.36976777703356780555748e+2,
+2.28754473395181007645155e+1,
-2.17925748738865115560082e-1,
+1.08314836272589368860689e-4
);
var
LDBL_EPSILON : extended;
//==============================================================================
// Γ関数の計算
// Γ関数の計算の説明は、ガンマ関数のページを参照してください。
//==============================================================================
function TForm1.Duplication_Formula(two_x: extended): extended;
var
x: extended;
g: extended;
n: integer;
begin
x := 0.5 * two_x;
n := round(two_x) - 1;
g := power(2.0, two_x - 1.0 - n);
g := g * power(2, n);
g := g / sqrt(pi);
g := g * xGamma_Function(x);
g := g * xGamma_Function(x + 0.5);
result := g;
end;
function TForm1.xGamma(x: extended): extended;
var
xx: extended;
temp: extended;
n, i: integer;
xx2 : extended;
begin
if x < 1 then xx := x + 1 else xx := x;
n := sizeof(a) div sizeof(extended);
if x > max_long_double_arg then begin
result := maxextended;
exit;
end;
if x > 900 then begin
result := Duplication_Formula(x);
exit;
end;
temp := 0;
for i := n - 1 downto 0 do temp := temp + a[i] / (xx + i);
temp := temp + 1;
xx2 := (xx - 0.5) / 2;
temp := temp * (power((g + xx - 0.5) / e, xx2) / exp_g_o_sqrt_2pi );
temp := temp * power((g + xx - 0.5) / e, xx2); // X64 オーバーフロー対策
if x < 1 then result := temp / x
else result := temp;
end;
function TForm1.xGamma_Function(x: extended): extended;
var
sin_x: extended;
rg: extended;
ix: int64;
begin
if x > 0 then
if x <= max_long_double_arg then begin
result := xGamma(x);
exit;
end
else begin
result := maxextended;
exit;
end;
if x > -LONG_MAX then begin
ix := round(x);
if x = ix then begin
result := maxextended;
exit;
end;
end;
sin_x := sin(pi * x);
if sin_x = 0 then begin
result := maxextended;
exit;
end;
if x < - max_long_double_arg - 1 then begin
result := 0;
exit;
end;
rg := xGamma(1 - x) * sin_x / pi;
if rg <> 0 then result := 1 / rg
else result := maxextended;
end;
//==============================================================================
// function Dirichlet_Eta_Function（s: double): double;
//
// このルーチンは、実数のEta関数を計算します
// eta(s)=（1-2^(1-s))zeta(s)
// ここで、zeta()はリーマンゼータ関数です。
// s >= 0の場合
// eta(s)= Sum(-1)^(k-1) (1 / k^s)
// ここで、合計はk = 1,...で加算 eta(1)= ln（2）
// s < 0の場合
// 反射式が使用されることに注意してください。
// s << 0の場合、
// eta（s）の傾きは、大きさが大きく、摂動が小さくなります。
// 引数は大きな絶対誤差を作成する可能性があります。
//==============================================================================
function TForm1.Dirichlet_Eta_Function(s: double): double;
var
x : double;
begin
x := xDirichlet_Eta_Function(s);
if abs(x) < MAXDouble then result := x
else
if x > 0 then result := Maxdouble
else result := -Maxdouble;
end;
function TForm1.xDirichlet_Eta_Function(s: extended): extended;
begin
if s >= 64.0 then begin
result := 1;
exit;
end;
result := 1.0 + xDirichlet_Eta_Star_Function(s);
end;
function TForm1.Dirichlet_Eta_Star_Function(s: double): double;
var
x : double;
begin
x := xDirichlet_Eta_Star_Function(s);
if abs(x) < MAXDouble then result := x
else
if x > 0.0 then result := MAXdouble
else result := -MAXdouble;
end;
function TForm1.xDirichlet_Eta_Star_Function(s : extended): extended;
var
reflection_coefficients : extended;
begin
if s >= 0.0 then
if s < 18.0 then begin
result := Alternating_Series_Convergence_Acceleration(s);
exit;
end
else begin
result := Sum_Reverse_Order(s);
exit;
end;
// Use the reflection formula s < 0
reflection_coefficients := Reflection_Coefficient(s);
result := reflection_coefficients * xDirichlet_Eta_Star_Function(1.0 - s)
+ (reflection_coefficients - 1.0);
end;
//==============================================================================
// function Reflection_Coefficient(s: extended): extended;
//
// ディリクレのイータ関数の反射公式は次のとおりです。
// eta(s)= {[2*(1-2^(1-s))/(1-2^s)]*Γ(1-s）
// /(2*pi)^(1-s)*cos((1-s*pi/2)}*eta(1-s)
// ここで、s < 0 は負の実数です。
//
// このルーチンは係数:
// [2*(1-2^(1-s))/(1-2^s)]*Γ(1-s)/(2*pi)^(1-s)*cos((1-s)*pi/2) を返します。
// ここで、s < 0は負の実数です。
// 引数
// s : extended
// ここで、s < 0です。
//
// 戻り値：
// eta（s）= c * eta（1-s） となる係数 c
//==============================================================================
function TForm1.Reflection_Coefficient(s: extended): extended;
var
one_s : extended;
k : integer;
v : extended;
c : extended;
x, temp: extended;
begin
one_s := 1.0 - s;
k := trunc(one_s / 4.0);
v := one_s - 4.0 * k;
c := cos(v * pi / 2);
if abs(c) < 1.8 * LDBL_EPSILON then begin
result := 0.0;
exit;
end;
temp := power(2.0, one_s);
x := (1.0 - temp) / (temp - 2.0);
x := x + x;
x := x * c;
x := x * xGamma_Function(one_s);
x := x / power(pi, one_s);
result := x;
end;
//==============================================================================
// function Alternating_Series_Convergence_Acceleration(s: extended): extended;
//
// 説明：
// このルーチンは、s <= 18 のEta関数を計算します。
// ゆっくりと収束します。
//
// テクニック：
// 交代級数が与えられた場合：合計は Sum(-1)^ka[k]、K = 0 ...で合計されます。
// ここで、すべてのkに対してa [k]> 0です。正の関数が存在する場合 w(x)
// [0,1]で、a[k] = I[0,1] x^kw(x)dx
// x^kw(x)のxに関して 0 から 1 まで積分
// sum(-1)^ka[k] = sum(-1)^kI[0,1] x^kw(x) dx
// = I[0,1] Sum(-1)^kx^kw(x) dx
// = I[0,1] w(x)/(1 + x) dx
//
// S = I [0,1] w(x)/（1 + x）dx とします。
// P(-1)！= 0となるような次数 n の多項式 P(x)が与えられた場合
// P(x)= p[0] - p[1] x + p[2] x^2 - ...（-1^np [n] x^n
// S[n] = I[0,1] (P(-1) - P(x)) w(x) / (1+x）dx / P(-1) を定義します。
// 次に、 S[n] = S-I [0,1] P（x）w（x）/（1 + x）dx / P（-1）
// つまり |S[n] -S| <= I[0,1] |P(x)| w(x) / (1 + x) dx / |P(x)|
// |S[n] -S| <= (M / |P(-1)|) S、ここで M = max |P(x)|、x in [0,1]
//
// 代数の後
// S[n] = (1/ P(-1)) Sum Sum p[j][(-1)^k a[k]]、
// 最初の合計 k = 0、...、n-1、2番目の合計は j = k + 1、...、n
//
// d[k] = Sum p[j] sumed for j = k + 1、...、n for k = n-1、...、0、
// S[n] = (1/c[0]）Sum d[k] (-1)^ka [k] summed k = 0、...、n-1
//
// ここでプログラムされた P(x) の選択は (-1)^n * T*[n](x) です。、
// T*[n]は 第1種の n次 の シフトされたチェビシェフ多項式です。
// P(x) をこのように選択すると、
// P(x) = Sum {(n / (n + j)) C(n+j、2j) 4^j (-x)^j}、
// 合計は j = 0、...、n and C(n+j,2j) = (n+j)!/[(2j)!(nj)!]
// and M = 1、P(-1) >= (1/2) (3 + sqrt(8))^n。
//
// 引数：
// s : extended
// s < 18
//
// 戻り値：
// eta(s)
//==============================================================================
function TForm1.Alternating_Series_Convergence_Acceleration(s: extended): extended;
const
d: array[0..28] of extended = (
1.362725501650887306817e+21, 1.362725501650887306816e+21,
1.362725501650887305248e+21, 1.362725501650886896000e+21,
1.362725501650844334208e+21, 1.362725501648488235008e+21,
1.362725501568066715648e+21, 1.362725499718371770368e+21,
1.362725469310199922688e+21, 1.362725096810094788608e+21,
1.362721590926752350208e+21, 1.362695647390018306048e+21,
1.362542007743905005568e+21, 1.361803869444099801088e+21,
1.358896740140251611136e+21, 1.349437033675348770816e+21,
1.323863206542645919744e+21, 1.266218975223368122368e+21,
1.157712186857668739072e+21, 9.872015194258554224640e+20,
7.640581139674368573440e+20, 5.220333434317674905600e+20,
3.061506212814840135680e+20, 1.496014168469232680960e+20,
5.884825485587356057600e+19, 1.781624012587768217600e+19,
3.882102878793367552000e+18, 5.404319552844595200000e+17,
3.602879701896396800000e+16
);
var
term: array[0..27] of extended;
sum: extended;
k: integer;
begin
k := 1;
while k <= 28 do begin
term[k-1] := d[k] * power(k + 1, -s);
inc(K);
term[k-1] := -d[k] * power(k + 1,-s);
inc(k);
end;
sum := term[27];
for k := 26 downto 0 do sum := sum + term[k];
sum := sum / d[0];
result := -sum;
end;
//==============================================================================
// function Sum_Reverse_Order（s: extended): extended;
//
// このルーチンは、
// s >= 18 の Eta関数を計算します(power()への呼び出しは28未満)が必要です)
//
// eta(s) = Sum (-1^(k-1) (1/ k^s
// 合計は k = 1、...で合計されます。
//
// 系列が交互になっているため、絶対誤差は最初に無視された項の絶対値未満です。
//
// 合計されるのは、連続する下限と上限を計算することによって決定されます
// 2つの境界が等しくなるまで（切り捨てエラー内で）計算され
// 次に、最終結果が最後の項から合計され、最初の項で終了します。
// （これは少し正確です）。
//
// 引数：
// s: extended;
//
// 戻り値：
// eta(s)
//==============================================================================
function TForm1.Sum_Reverse_Order(s: extended): extended;
var
term: array[0..29] of extended;
lower_bound : extended;
upper_bound : extended;
sum : extended;
k : integer;
begin
lower_bound := -power(2.0, -s);
term[0] := lower_bound;
term[1] := power(3.0, -s);
upper_bound := term[1] + lower_bound;
if lower_bound = upper_bound then begin
result := upper_bound;
exit;
end;
k := 4;
while k < 32 do begin
term[k - 2] := -power(k, -s);
lower_bound := upper_bound + term[k - 2];
if lower_bound = upper_bound then break;
inc (k);
term[k - 2] := power(k, -s);
upper_bound := lower_bound + term[ k- 2];
if lower_bound = upper_bound then break;
inc(k);
end;
k := k - 2;
sum := term[k];
for K := k - 1 downto 0 do sum := sum + term[k];
result := sum;
end;
//==============================================================================
// ディリクレのイータ関数計算とグラフ表示
//==============================================================================
procedure TForm1.BitBtn1Click(Sender: TObject);
const
n = 200;
var
ans, anss, s: double;
min, max : double;
x, dx : double;
ch : integer;
begin
val(LabeledEdit1.Text, s, ch);
if ch <> 0 then begin
application.MessageBox('注意','Sの値に間違いがあります。',0);
exit;
end;
Memo1.Clear;
ans := Dirichlet_Eta_Function(s);
Memo1.Lines.Append('s = ' + floatTostr(s));
Memo1.Lines.Append(' η(s) = ' + floatTostr(ans));
anss := Dirichlet_Eta_Star_Function(s);
Memo1.Lines.Append('η*(s) = η(s) - 1');;
Memo1.Lines.text := Memo1.Lines.text + ' η*(s) = ' + floatTostr(anss);
Series1.Clear;
Series2.Clear;
min := s - 10;
if min < - 1 then begin
min := s - 0.5;
if s >= 0 then
if s - min < 10 then min := -0.5;
end;
max := s + 10;
dx := (max - min) / n;
for ch := 0 to n do begin
x := dx * ch + min;
anss := Dirichlet_Eta_Function(x);
Series1.AddXY(x, anss);
end;
Series2.AddXY(s, ans);
end;
//==============================================================================
// LDBL_EPSILONの設定
//==============================================================================
procedure TForm1.FormCreate(Sender: TObject);
var
one, tmp : extended;
begin
one := 1;
tmp := 1;
repeat
tmp := tmp / 2;
until one + tmp = 1;
LDBL_EPSILON := tmp * 2;
memo1.Clear;
end;
end.
次のリストは、C言語辞典のゼータ関数をDelphiに変換し、更にガンマ関数を追加して、負数の値にたいして精度を向上したゼータ関数からイータ関数を計算するプログラムです。
unit main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, system.Math, Vcl.StdCtrls, Vcl.Buttons,
Vcl.ExtCtrls, VclTee.TeeGDIPlus, VCLTee.Series, VCLTee.TeEngine,
VCLTee.TeeProcs, VCLTee.Chart, system.UITypes;
type
TForm1 = class(TForm)
Memo1: TMemo;
BitBtn1: TBitBtn;
LabeledEdit1: TLabeledEdit;
Chart1: TChart;
Series1: TLineSeries;
Series3: TPointSeries;
procedure FormCreate(Sender: TObject);
procedure BitBtn1Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
function Gamma(x: Extended): Extended; // ガンマ関数
const
// ベルヌーイ数B(2)、B(4)、B(6)、...、B(16)
B : array[1..8] of Extended = ( 1.0 / 6, // B(2)
-1.0 / 30,
1.0 / 42,
-1.0 / 30,
5.0 / 66,
-691.0 / 2730,
7.0 / 6,
-3617.0 / 510); // B(16)
m: integer = sizeof(B) div sizeof(Extended);
var
i : integer;
sum : Extended;
xx, v : Extended;
xj : Extended;
lng : Extended;
ln_sqrt_2pi: Extended;
begin
ln_sqrt_2pi := ln(sqrt(2 * pi)); // = ln(2 * pi) / 2;
sum := 0;
v := 1;
while x < m do begin
v := v * x;
x := x + 1;
end;
xx := x * x;
xj := x;
lng := ln_sqrt_2pi - x + (x - 0.5) * ln(x) - ln(v);
for i := 1 to m do begin
sum := sum + B[i] / (i * 2 * (2 * i - 1)) / xj;
xj := xj * xx;
end;
result := exp(sum + lng);
end;
function Riemann_Zeta_Function(x: extended): extended;
const
coef: array[0..19] of extended = (
8.333333333333333333333333333e-2, // 1/12
-1.388888888888888888888888889e-3, // -1/720
3.306878306878306878306878307e-5, // 1/30240
-8.267195767195767195767195767e-7, // -1/1209600
2.087675698786809897921009032e-8, // 1/47900160
-5.284190138687493184847682202e-10,
1.338253653068467883282698098e-11,
-3.389680296322582866830195391e-13,
8.586062056277844564135905450e-15,
-2.174868698558061873041516424e-16,
5.509002828360229515202652609e-18,
-1.395446468581252334070768626e-19,
3.534707039629467471693229977e-21,
-8.953517427037546850402611251e-23,
2.267952452337683060310950058e-24,
-5.744790668872202445263829503e-26,
1.455172475614864901866244572e-27,
-3.685994940665310178130050728e-29,
9.336734257095044668660153106e-31,
-2.365022415700629886484029550e-32
);
N = 8;
var
i: integer;
powNx, w, z, zprev: extended;
begin
if x = 1 then begin
result := infinity;
exit;
end;
z := 1;
for i := 2 to N - 1 do begin
zprev := z;
z := z + power(i, -x);
if z = zprev then begin
result := z;
exit;
end;
end;
powNx := power(N, x);
w := x / (N * powNx);
z := z + 0.5 / powNx + N / ((x - 1) * powNx) + coef[0] * w;
i := 1;
while (i < 20) and (z <> zprev) do begin
w := w * (x + 2 * i - 1) * (x + 2 * i) / (N * N);
zprev := z;
z := z + coef[i] * w;
inc(i);
end;
result := z;
end;
function RiemannFunction(s: Extended): Extended;
var
term : Extended;
ints : integer;
def : Extended;
begin
ints := trunc(s);
def := s - ints;
if (s <= -2) and (def = 0) and (ints mod 2 = 0) then
result := 0
else
if s < -1 then begin
term := power(2, s) * power(pi, s - 1) * sin(pi * s / 2);
result := term * Riemann_Zeta_Function(1 - s) * Gamma(1 - s);
end
else
result := Riemann_Zeta_Function(s);
end;
function Dirichlet_Eta_Function(x: extended): extended;
begin
if x = 1 then begin
result := ln(2);
exit;
end;
result := (1 - power(2, 1- x)) * RiemannFunction(x);
end;
//------------------------------------------------------------------------------
// 計算の実行とグラフ作成
//------------------------------------------------------------------------------
procedure TForm1.BitBtn1Click(Sender: TObject);
const
n = 400;
var
eta: extended;
star: extended;
s : extended;
i : integer;
min, max, dx: extended;
begin
val(LabeledEdit1.Text, s, i);
if i <> 0 then begin
MessageDlg('xの値に間違いがあります。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
eta := Dirichlet_Eta_Function(s);
star := eta - 1;
Memo1.Clear;
Memo1.Lines.Add('η関数');
Memo1.Lines.Add('η(s) = ' + floatTostr(eta));
Memo1.Lines.Text := Memo1.Lines.Text + 'η(s) - 1 = ' + floatTostr(star);
Series1.Clear;
Series3.Clear;
Series3.AddXY(s, eta);
min := s - 10;
if min < - 1 then begin
min := s - 0.5;
if s >= 0 then
if s - min < 10 then min := -0.5;
end;
max := s + 9;
dx := (max - min) / n;
for i := 0 to n do begin
s := i * dx + min;
eta := Dirichlet_Eta_Function(s);
Series1.AddXY(s, eta);
end;
end;
//-------------------------------------------
procedure TForm1.FormCreate(Sender: TObject);
begin
Memo1.Clear;
end;
end.
dirichlet_eta_function.zip