ディリクレラムダ関数
ディリクレラムダ関数λ(s)は、s=1に単純な極を持つ複素変数sの有理型複素数値関数です。この関数は、実数 s＞1 の λ(s) を次の様に定義することによって構築されます。
λ(s) = Σk=0∞1/（2k + 1）s
次に、s = 1で保存して、複素平面の残りの部分に分析的に続行します、ディリクレラムダ関数は、リーマンゼータ関数で次のように表すことができます。
λ（s）=（1-2-s）ζ（s）
ディリクレラムダ関数は、インターネットで探してもあまり資料が見つかりません。
最初のプログラムはMathematics Source Library C & ASMにある Dirichlet Lambda function をDelphi に変換したものですが、リーマンのζ 関数が計算出来れば、簡単にディリクレλ 関数が計算出来ます。
プログラム1
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;
LDBL_EPSILON : extended;
implementation
{$R *.dfm}
uses dirichlet_eta, Riemann_Zeta;
//==================================================================================================
// function Lambda_Star_Negative_Arg（s: extended): extended
//
// このルーチンは、s <0 のディリクレラムダスター関数を計算します。
//
// lambda *(s)=[(2^s -1)/（2^s] zeta*(s)- 1/2^s、
// zeta*(s) はzeta*関数です。
//
//==================================================================================================
function Lambda_Star_Negative_Arg(s: extended): extended;
var
zeta_star: extended;
lambda_star: extended;
two_minus_s: extended;
two_minus_s1: extended;
begin
zeta_star := xRiemann_Zeta_Star_Function(s);
if abs(zeta_star) = MAXExtended then begin
result := -zeta_star;
exit;
end;
two_minus_s := power(2.0, -s);
two_minus_s1 := two_minus_s - 1.0;
if two_minus_s1 < 1 then
two_minus_s1 := 1; // two_minus_s1が1と成る時1の筈だが最終ビット当たりに誤差が出ます
if abs(zeta_star) > MAXExtended / two_minus_s1 then begin
if zeta_star < 0.0 then result := MAXExtended else result := -MAXExtended;
exit;
end;
lambda_star := (1.0 - two_minus_s) * zeta_star;
if lambda_star > 0.0 then begin
result := lambda_star - two_minus_s;
exit;
end;
if lambda_star < two_minus_s - MAXExtended then begin
result := -MAXExtended;
exit;
end;
result := lambda_star - two_minus_s;
end;
//==================================================================================================
// function Lambda_Star_Positive_Arg(s: extended): extended;
// このルーチンは、s >= 0 のディリクレラムダスター関数を計算します。
//==================================================================================================
function Lambda_Star_Positive_Arg(s: extended): extended;
const
LDBL_MAX_EXP = 16384; // 64ビット(double)の時は1024
var
zeta_star: extended;
two_s: extended;
begin
if s = 1.0 then begin
result := MAXExtended;
exit;
end;
zeta_star := xRiemann_Zeta_Star_Function(s);
if s > LDBL_MAX_EXP then begin
result := zeta_star;
exit;
end;
two_s := power(2.0, s);
result := ((two_s - 1.0) * zeta_star - 1.0) / two_s;
end;
//==================================================================================================
// function xDirichlet_Lambda_Star_Function(s: extended): extended;
//
// このルーチンは、次のように定義されたディリクレラムダスター関数を計算します。
// λ*(s) = λ-1、
// lambda()はディリクレラムダ関数です。
//
// lambda*(s)= [zeta*(s）+ eta*(s)] / 2 であることに注意してください。
// zeta*はリーマンゼータスター関数と eta*はディリクレのイータスター関数です。
// sの値によって計算を選択します。
//==================================================================================================
function xDirichlet_Lambda_Star_Function(s: extended): extended;
begin
if s = 0.0 then begin
result := -1.0;
exit;
end;
if s > 0.0 then
result := Lambda_Star_Positive_Arg(s)
else
result := Lambda_Star_Negative_Arg(s);
end;
//==================================================================================================
// function Dirichlet_Lambda_Star_Function(s: double): edouble;
//
// このルーチンは、次のように定義されたディリクレラムダスター関数を計算します。
// λ*(s) = λ-1、
// lambda()はディリクレラムダ関数です。
//
// lambda*(s)= [zeta*(s）+ eta*(s)] / 2 であることに注意してください。
// zeta*はリーマンゼータスター関数と eta*はディリクレのイータスター関数です。
//==================================================================================================
function Dirichlet_Lambda_Star_Function(s: double): double;
var
lambda_star: extended;
begin
if s = 1 then begin
result := MAXDouble;
exit;
end;
lambda_star := xDirichlet_Lambda_Star_Function(s);
if abs(lambda_star) >= MAXDouble then begin
if lambda_star < 0.0 then result := -MAXDouble else result := MAXDouble;
exit;
end;
result := lambda_star;
end;
//==================================================================================================
// function xDirichlet_Lambda_Function(s: extended): extended;
//
// ディリクレラムダ関数は次のように定義されます
// lambda(s)= Sum(1 /(2k + 1)^s）、 k = 0、...で合計されます。
// sはRe(s)> 1の複素数です。
// 次に、解析的に複素平面の残りの部分に進みます。
// s = 1の場合、LDBL_MAXが返されます。//
// lambda（s）> LDBL_MAXの場合、LDBL_MAXが返され、
// lambda（s）<-LDBL_MAXの場合、-LDBL_MAXが返されます。
// lim lambda(s)= -inf、左からs-> 1、lim lambda(s)= inf 右からs-> 1。
//
// lambda(s)= [zeta(s)+ eta(s)] / 2であることに注意してください。
// ここで、zetaはリーマンゼータ関数とetaはディリクレのイータ関数です。
//
// lambda(s)= 1 + lambda_star(s)
//==================================================================================================
function xDirichlet_Lambda_Function(s : extended): extended;
var
lambda_star: extended;
begin
if s = 1 then begin
result := MAXExtended;
exit;
end;
lambda_star := xDirichlet_Lambda_Star_Function(s);
if abs(lambda_star) = MAXExtended then begin
result := lambda_star;
exit;
end;
result := 1.0 + lambda_star;
end;
//==================================================================================================
// function Dirichlet_Lambda_Star_Function(s: double): double;
//
// このルーチンは、次のように定義されたディリクレラムダスター関数を計算します。
// λ*(s) = λ(s)-1、
// lambda()はディリクレラムダ関数です。//
//
// lambda*(s)= [zeta*(s)+ eta*(s)] / 2 であることに注意してください。
// zeta*は リーマンゼータ*関数とeta*は ディリクレのイータスター関数です。
//==================================================================================================
function Dirichlet_Lambda_Function(s : double): double;
var
lambda: extended;
begin
if s= 1.0 then begin
result := MAXDouble;
exit;
end;
lambda := xDirichlet_Lambda_Function(s);
if abs(lambda) >= MAXDouble then begin
if lambda < 0.0 then result := -MAXDouble else result := MaxDouble;
exit;
end;
result := lambda;
end;
//--------------------------------------------------------------------------------------------------
// 計算の実行とグラフ作成
//--------------------------------------------------------------------------------------------------
procedure TForm1.BitBtn1Click(Sender: TObject);
const
n = 2000;
var
Lambda, Lambdab: double;
star: double;
s : double;
i : integer;
min, max, dx: double;
begin
val(LabeledEdit1.Text, s, i);
if i <> 0 then begin
MessageDlg('xの値に間違いがあります。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
Lambda := Dirichlet_Lambda_Function(s);
star := Dirichlet_Lambda_Star_Function(s);
Memo1.Clear;
Memo1.Lines.Add('λ関数');
// s = 1の時はMaxDoubleになるので∞を表示
if s <> 1 then begin
Memo1.Lines.Add('λ(x) = ' + floatTostr(Lambda));
Memo1.Lines.Text := Memo1.Lines.Text + 'λ(x) - 1 = ' + floatTostr(star);
// Memo1.Lines.Add('λ(x) - 1 = ' + floatTostr(star));
end
else begin
Memo1.Lines.Add('λ(x) = ∞');
Memo1.Lines.Text := Memo1.Lines.Text + 'λ(x) - 1 = ∞';
end;
Series1.Clear;
Series3.Clear;
if Lambda > 3 then Lambda := 3;
if Lambda < -2 then Lambda := -2;
Series3.AddXY(s, Lambda);
min := - 10;
if s < min then min := s;
max := 5;
if s > max then max := s;
dx := (max - min) / n;
Lambdab := 1;
for i := 0 to n do begin
s := i * dx + min;
Lambda := Dirichlet_Lambda_Function(s);
if Lambda > 3 then Lambda := 3;
if (Lambdab = -2) and (Lambda = 3) then begin
Series1.AddXY(s, Lambdab, '', clBlue);
Series1.AddXY(s, Lambda, '', clBlue);
end;
if Lambda < -2 then Lambda := -2;
if (Lambda < 3) and (Lambda > -2) then Series1.AddXY(s, Lambda, '', clBlue)
// グラフのベースの色と同じにして見えなくします。
else Series1.AddXY(s, Lambda, '', clBtnFace);
Lambdab := Lambda;
end;
end;
//-------------------------------------------
// extended の Epsilon 設定
//-------------------------------------------
procedure TForm1.FormCreate(Sender: TObject);
var
d, m: extended;
begin
Memo1.Clear;
LDBL_EPSILON := 1;
d := 1;
repeat
LDBL_EPSILON := LDBL_EPSILON / 2;
m := d + LDBL_EPSILON;
until m = d;
LDBL_EPSILON := LDBL_EPSILON * 2;
end;
end.
unit Riemann_Zeta; interface function xRiemann_Zeta_Star_Function(s: extended): extended; implementation uses system.math, dirichlet_eta; //------------------------------------------------------------------------------ // function xRiemann_Zeta_Star_Function(s : extended): extended; // // Description: // This routine calculates the Riemannt zeta star function defined as // zeta*(s) = zeta(s) - 1, // where zeta() is the Riemann zeta function for real s != 1. //------------------------------------------------------------------------------ function xRiemann_Zeta_Star_Function(s: extended): extended; var two_s: extended; eta: extended; temp: extended; begin two_s := power(2.0, s - 1.0); eta := xDirichlet_Eta_Star_Function(s); if s = 1.0 then begin result := MAXDouble; exit; end; if s < 0.0 then begin two_s := 1.0 / two_s; result := (eta + two_s) / (1.0 - two_s); exit; end; temp := two_s - 1.0; result := (two_s * eta + 1.0) / temp; end; end.
unit Gamma_Function;
interface
function xGamma_Function(x: extended): extended;
implementation
uses system.math;
const
LONG_MAX = 2147483647;
max_long_double_arg: extended = 1755.5;
g = 9.65657815377331589457187;
e = 2.71828182845904523536028747;
exp_g_o_sqrt_2pi = 6.23316569877722552586386e+3;
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
);
function xGamma(x: extended): extended; forward;
//=============================================================================
//
// この関数はLanczosの式を使用して実際のGamma(x)を計算します。
// xの範囲は -(max_long_double_arg-1) < x < max_long_double_arg。
// ガンマ関数は複素平面で有理型であり、正でない整数の値
// xaの正の整数または半分の正の整数をテストすると、
// 最大絶対相対誤差は約3.5e-16です。
//
// x > max_long_double_argの場合、lnGamma(x)を使用する必要があります。
// x < 0 の場合、ln(Gamma(x)) は複素数になる可能性があることに注意してください。
//
// 戻り値
// xが正で、max_long_double_argより小さい場合、Gamma(x)
// が返され、x > max_long_double_argの場合、maxextendedが返されます。
// xが正でない整数の場合、maxextendedが
// 返されます（Gamma(x)は値の両側で符号を変更することに注意してください）。
// xが非正の非整数の場合でx>-(max_long_double_arg + 1）の場合
// Gamma(x)が返されます。それ以外の場合は0.0が返されます
//
//=============================================================================
function 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;
//=============================================================================
//
// この関数は倍数公式を使用してGamma(two_x)を返します。
//
//=============================================================================
function 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;
//=============================================================================
//
// この関数はLanczosの式を使用して実際のGamma(x)を計算します。
// x、ここで0<x<=900。900<x<1755.5の場合、倍数公式を使用します。
// 関数power()。結果の相対誤差は約10^-16です。x = 0付近を除きます。
// x>1755.5の場合、lnGamma(x)を計算する必要があります。
// xが正で1755.5未満の場合、Gamma(x)が返され、
// x>1755.5の場合、maxextendedが返されます。
//
//=============================================================================
function 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;
end.
プログラム2
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 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;
// ========================================================================
// x が -5 より小さくなりと誤差が大きくなります
//=========================================================================
function Riemann_Zeta_Function(x: extended): extended;
const
N = 8;
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
);
var
i, N2 : 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;
N2 := N * N;
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) / (N2);
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 // -s 時誤差が大きくなるのでζを+で計算します
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_Lambda_Function(x: extended): extended;
begin
result := (1 - power(2, -x)) * RiemannFunction(x);
// result := power(2, -x) * (power(2, x) - 1) * Riemann_Zeta_Function(x);
end;
//------------------------------------------------------------------------------
// 計算の実行とグラフ作成
//------------------------------------------------------------------------------
procedure TForm1.BitBtn1Click(Sender: TObject);
const
n = 2000;
var
Lambda, Lambdab: 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;
Lambda := Dirichlet_Lambda_Function(s);
star := Lambda - 1;
Memo1.Clear;
Memo1.Lines.Add('λ関数');
Memo1.Lines.Add('λ(s) = ' + floatTostr(Lambda));
Memo1.Lines.Text := Memo1.Lines.Text + 'λ(s) - 1 = ' + floatTostr(star);
Series1.Clear;
Series3.Clear;
if Lambda > 3 then Lambda := 3;
if Lambda < -2 then Lambda := -2;
Series3.AddXY(s, Lambda);
min := -10;
if s < min then min := s;
max := 5;
if s > max then max := s;
dx := (max - min) / n;
Lambdab := 1;
for i := 0 to n do begin
s := i * dx + min;
Lambda := Dirichlet_Lambda_Function(s);
if Lambda > 3 then Lambda := 3;
if (Lambdab = -2) and (Lambda = 3) then begin
Series1.AddXY(s, Lambdab, '', clBlue);
Series1.AddXY(s, Lambda, '', clBlue);
end;
if Lambda < -2 then Lambda := -2;
if (Lambda < 3) and (Lambda > -2) then Series1.AddXY(s, Lambda, '', clBlue)
// グラフのベースの色と同じにして見えなくします。
else Series1.AddXY(s, Lambda, '', clBtnFace);
Lambdab := Lambda;
end;
end;
end.
dirichlet_Lambda_function.zip