2024/03/25 多倍長計算のルーチンを見直し演算を早くしました。
2024/03/20 多倍長計算のFA,MVG, PVGの配列の長さを修正しました。
2023/10/08
複素数の計算ができるプログラムを追加しました、又、多倍長演算とし、50桁の精度で答えがでる様にしてみました。
2021/09/13 次数が実数の場合の計算で、計算精度が高くなるようにした、計算式の検討をしてみました。
Xの値が大きくなった時、二けた程度、精度が高くなります。
2020/07/01 次数の値がマイナスの時の計算を考慮していませんでしたので修正しました。
第１種ベッセル関数
ベッセル関数の詳細は、インターネットで参照してください。
周波数フーリエ変換のカイザー窓関数で、第１種変形ベッセル関数使用されていたので、取り上げました。
Γ関数を使用すれば、次数nの値は、整数でなくても計算可能です。
Γ関数はΓ(n+m+1)だけ使用すれば、OKです。
整数の場合は、階乗でΓ値が計算可能です。
次数が整数で負数の時は、左図の関係となります。
次の計算式は、次数が実数の場合で、計算精度が高くなるようにしたものです。
変形の場合は、Yの符号を逆(Y=-(x/2)2)にするだけです。
次の実行画面は、次数ｎの値が整数の時のグラフです。
∑の計算で、m=0から無限大まで加算するようになっていますが、実用的には、m=0から30程度でも良いようです。
Double 32ビットの浮動小数点では、171の階乗までしか計算出来ません、171以上にするとオーバーフローします。
また、Xの値が大きくなると、演算有効桁数の不足により、誤差が大きくなり正しい結果が得られなくなります。
次のプログラムは、次数が正の整数の時のプログラムです。
// Xの値 積分計算の上限値MNによってオーバーフローを発生するので // コンパイル設定にオーバーフローチックが必須です。 unit Main; interface uses Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, Vcl.ExtCtrls, VCLTee.TeeProcs, VCLTee.Chart, system.UITypes; type TForm1 = class(TForm) Button1: TButton; Chart1: TChart; Series1: TLineSeries; CheckBox1: TCheckBox; LabeledEdit2: TLabeledEdit; Series2: TLineSeries; LabeledEdit3: TLabeledEdit; Series3: TLineSeries; CheckBox2: TCheckBox; Label1: TLabel; procedure Button1Click(Sender: TObject); procedure FormCreate(Sender: TObject); private { Private 宣言 } function gamma(z: double): double; procedure Besselgraph(N: integer; Xmax: double; flag: boolean); function Bessel(X: double; N: integer; flag: boolean): double; function BesselG(X: double; v: double; flag: boolean): double; public { Public 宣言 } end; var Form1: TForm1; implementation {$R *.dfm} uses math; var MN : integer; // 積分上限値 //----------------------------------------------------------------------------- // Gamma_Function(x) function TForm1.gamma(z: double): double; const P: array[0..7] of double = ( 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7); max_double_z: double = 171.0; var i : integer; a, t : double; begin if z > max_double_z then begin result := maxdouble; exit; end; if z < 0.5 then begin result := pi / (sin(pi * z) * gamma(1 - Z)); exit; end; z := z - 1; a := 0.99999999999980993; for i := 0 to length(p) - 1 do a := a + p[i] / (z + i + 1); t := z + length(p) - 0.5; result := sqrt(2 * pi) * power(t, z + 0.5) * exp(-t) * a; end; //***************************************************************************** //----------------------------------------------------------------------------- // 第一種ベッセル関数 Γ関数使用 // X 値 // v 次数 // flag false 通常 true 変形 function TForm1.BesselG(X: double; v: double; flag: boolean): double; var M : Integer; B, S : double; begin S := 0; for M := 0 to MN do begin B := 1 / Gamma(M + 1) / Gamma(M + v + 1); // B = 1 / (Γ(m + 1)Γ(m + n + 1)) if not flag then S := S + power(-1, M) * B * power(X / 2, 2 * M + v) else S := S + B * power(X / 2, 2 * M + v) end; result := S; end; // 第一種ベッセル関数 // X 値 // N 次数 // flag false 通常 true 変形 function TForm1.Bessel(X: double; N: integer; flag: boolean): double; var J, K: Integer; B, S, Jf: double; begin B := 1.0; if N > 1 then begin for K := 1 to N do B := B / K; end; S := B * power((X / 2), N); // j = 0 for J := 1 to MN do begin Jf := J; // オーバーフロー対策で浮動小数点を使用します B := B / (Jf * (Jf + N)); // b / m!(m + n)! if not flag then S := S + power(-1, J) * B * power(X / 2, 2 * J + N) else S := S + B * power(X / 2, 2 * J + N) end; result := S; end; // 第一種ベッセル関数 グラフ表示 // N 次数 // Xmax グラフ最大値 // flag false 通常 true 変形 procedure TForm1.Besselgraph(n: integer; Xmax: double; flag: boolean); const NMAX = 1000; var Y, X: array of double; I : integer; XMIN: double; begin setlength(X, NMAX); setlength(Y, NMAX); XMIN := 0; for I := 0 to NMAX - 1 do begin X[I] := XMIN + (XMAX - XMIN) * I / NMAX; if checkbox2.Checked = False then Y[I] := Bessel(x[I], n, flag) // ベッセル関数 else Y[I] := BesselG(x[I], n, flag); // ガンマ関数使用ベッセル関数 end; for I := 0 to NMAX - 1 do begin if n = 0 then Series1.AddXY(X[I], Y[I]); if n = 1 then Series2.AddXY(X[I], Y[I]); if n = 2 then Series3.AddXY(X[I], Y[I]); end; end; // ベッセル関数計算 procedure TForm1.Button1Click(Sender: TObject); var ch: integer; X : double; Flg: boolean; begin val(LabeledEdit2.Text, X, Ch); if ch <> 0 then begin MessageDlg('X値に間違いがあります。', mtInformation, [mbOk], 0, mbOk); end; val(LabeledEdit3.Text, MN, Ch); if ch <> 0 then begin MessageDlg('積分上限に間違いがあります。', mtInformation, [mbOk], 0, mbOk); end; Chart1.Title.Text.Clear; if checkbox1.Checked = False then begin Flg := False; Chart1.Title.Text.Add('第１種ベッセル関数'); end else begin Flg := True; Chart1.Title.Text.Add('第１種変形ベッセル関数'); end; Series1.Clear; Series2.Clear; Series3.Clear; Besselgraph(0, X, Flg); Besselgraph(1, X, Flg); Besselgraph(2, X, Flg); end; procedure TForm1.FormCreate(Sender: TObject); begin Label1.Caption := '値を大きくし過ぎると' + #13#10 + 'オーバーフローします。'; end; end.
次のプログラムは次数vが実数の場合の計算です。
通常の計算と、精度を高くした計算の両方が出来ます。
右側の実行画面は、正しい値との差を表示しています。
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, VclTee.TeeGDIPlus,
VCLTee.Series, VCLTee.TeEngine, Vcl.ExtCtrls, VCLTee.TeeProcs, VCLTee.Chart;
type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
LabeledEdit1: TLabeledEdit;
LabeledEdit2: TLabeledEdit;
Chart1: TChart;
Series1: TLineSeries;
Series2: TPointSeries;
Button2: TButton;
CheckBox1: TCheckBox;
CheckBox2: TCheckBox;
procedure Button1Click(Sender: TObject);
procedure Button2Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses system.Math;
var
// 第一種変形ベッセル関数 I0(x) X = 0 ～15 の値
I0x: array [0..15] of extended = ( Maxdouble,
1.266065877752008335598,
2.279585302336067267437,
4.880792585865024085611,
11.301921952136330449636,
27.23987182360444689454,
67.23440697647797532619,
168.5939085102896988573,
427.5641157218047851774,
1093.588354511374695846,
2815.71662846625447147,
7288.489339821248106181,
18948.92534929630886121,
49444.48958221757299925,
129418.5627006485636646,
339649.3732979138795217
);
// 第一種ベッセル関数 J0(x) X = 0 ～15 の値
J0x: array [0..15] of extended = ( 1,
0.7651976865579665514497,
0.2238907791412356680518,
-0.2600519549019334376242,
-0.3971498098638473722866,
-0.1775967713143383043474,
0.1506452572509969316623,
0.3000792705195555966503,
0.1716508071375539060909,
-0.09033361118287613433595,
-0.2459357644513483351978,
-0.1711903004071960883458,
0.04768931079683353662381,
0.2069261023770678109966,
0.1710734761104586590631,
-0.01422447282678077323386
);
//ベルヌーイ数B(2)、B(4)、B(6)、...、B(20)
B : array[1..10] 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,
43867.0 / 798,
-174611.0 / 330); // B(20)
function loggamma(x: extended): extended; // ガンマ関数の対数
var
log_2pis2 : extended;
v, w: extended;
begin
log_2pis2 := ln(2 * pi) / 2;
v := 1;
while x < 10 do begin // ベルヌーイ数の数によりx値が小さい時補正します
v := v * x;
x := x + 1;
end;
w := 1 / (x * x);
result := (((((((((
B[10] / (20 * 19) * w + B[9] / (18 * 17)) * w
+ B[8] / (16 * 15)) * w + B[7] / (14 * 13)) * w
+ B[6] / (12 * 11)) * w + B[5] / (10 * 9)) * w
+ B[4] / ( 8 * 7)) * w + B[3] / ( 6 * 5)) * w
+ B[2] / ( 4 * 3)) * w + B[1] / ( 2 * 1)) / x
+ log_2pis2 - ln(v) - x + (x - 1 / 2) * ln(x);
end;
//--------------------------------------------------------------
function gamma(x: extended): extended; // ガンマ関数
begin
if x < 0 then begin
result := pi / (sin(pi * x) * exp(logGamma(1 - x)));
end
else begin
result := exp(LogGamma(x));
end;
end;
//==============================================================================
//----------------------- new --------------------------------------------------
// X 値
// v 次数
// flag false 通常 true 変形
function Jvx(v, x: extended; flag: boolean): extended;
var
k : integer;
g, j, bkk, n: extended;
function Bk(k : integer): extended;
var
kn : integer;
Bn, Y, Zk: extended;
begin
Y := -(x / 2) * (x / 2);
if flag then Y := -Y;
Bn := 1;
for kn := 1 to k do begin
Zk := Y / (kn * (kn + v));
Bn := Zk * Bn;
end;
result := Bn
end;
begin
g := gamma(1 + v);
n := round(v);
if (n = v) and (n < 0) then v := abs(v);
bkk := 0;
for k := 0 to 100 do bkk := bkk + Bk(k);
j := power(x / 2, v) / g * bkk;
if (not flag) and (abs(n) = v) and (n < 0) then j := j * Power(-1, v);
result := j;
end;
//------------------- old -----------------------------------------------------
// X 値
// v 次数
// flag false 通常 true 変形
function Bessel1(X: extended; v: extended; flag: boolean): extended;
const
Mn = 100;
var
M : Integer;
B, S, n: extended;
begin
n := round(v);
if (n = v) and (n < 0) then v := abs(v);
S := 0;
for M := 0 to MN do begin
B := 1 / Gamma(M + 1) / Gamma(M + v + 1); // B = 1 / (Γ(m + 1)Γ(m + n + 1))
if not flag then
S := S + power(-1, M) * B * power(X / 2, 2 * M + v)
else
S := S + B * power(X / 2, 2 * M + v)
end;
if (not flag) and (abs(n) = v) and (n < 0) then S := S * Power(-1, v);
result := S;
end;
//==============================================================================
procedure TForm1.Button1Click(Sender: TObject);
var
jk : extended;
v, x, dx: extended;
i, n : integer;
flag : boolean;
begin
flag := false;
if checkbox2.Checked = true then flag := true;
memo1.Clear;
val(labelededit1.Text, v, i);
if i <> 0 then begin
memo1.Lines.Append('vの値に間違いがあります。');
exit;
end;
val(labelededit2.Text, x, i);
if i <> 0 then begin
memo1.Lines.Append('xの値に間違いがあります。');
exit;
end;
{$IFDEF CPUX64}
Memo1.Lines.Append('X64 double');
{$ELSE}
Memo1.Lines.Append('X86 extended');
{$ENDIF}
if checkbox1.Checked = true then jk := Bessel1(x, v, flag)
else jk := jvx(v, x, flag);
if checkbox2.Checked = false then memo1.Lines.Append('第一種ベッセル関数')
else memo1.Lines.Append('第一種変形ベッセル関数');
memo1.Lines.Append('');
if checkbox2.Checked = false then memo1.Lines.Append('Jv(x) = ' + floatTostr(jk))
else memo1.Lines.Append('Iv(x) = ' + floatTostr(jk));
chart1.Title.Text.Clear;
if checkbox2.Checked = false then chart1.Title.Text.Append('第一種ベッセル関数')
else chart1.Title.Text.Append('第一種変形ベッセル関数');
chart1.LeftAxis.Title.Caption := '';
Series1.Clear;
Series2.Clear;
Series2.AddXY(x, jk);
dx := 0.02;
n := 750;
if checkbox2.Checked = true then n := round((x + 1) / dx);
for i := 0 to n do begin
x := i * dx;
if checkbox1.Checked = true then jk := Bessel1(x, v, flag)
else jk := jvx(v, x, flag);
Series1.AddXY(x, jk);
end;
end;
// J0(x) X=0～15の値の誤差の計算
procedure TForm1.Button2Click(Sender: TObject);
var
flg: boolean;
i, p, l: integer;
jk: extended;
def : extended;
defs, re, ex, ads, no : string;
begin
memo1.Clear;
Series1.Clear;
Series2.Clear;
{$IFDEF CPUX64}
chart1.LeftAxis.Title.Caption := 'xxE12';
Memo1.Lines.Append('X64 double');
{$ELSE}
chart1.LeftAxis.Title.Caption := 'xxE15';
Memo1.Lines.Append('X86 extended');
{$ENDIF}
flg := False;
if checkbox2.Checked = true then flg := true;
for i := 0 to 15 do begin
if checkbox1.Checked = true then jk := Bessel1(i, 0, flg)
else jk := jvx(0, i, flg);
if flg and (i = 0) then def := 0
else begin
if flg then def := abs(I0x[i] - jk)
else def := abs(J0x[i] - jk);
end;
{$IFDEF CPUX64}
Series1.AddXY(i, def * 1E12);
{$ELSE}
Series1.AddXY(i, def * 1E15);
{$ENDIF}
if def <> 0 then
defs := floatTostr(def)
else defs := ' 0';
l := length(defs);
p := pos('E', defs, 1);
re := copy(defs, 0, 4);
ex := copy(defs, p, l - p + 1);
no := inttostr(i);
if i < 10 then no := ' ' + no;
if i mod 2 = 0 then
ads := no + ' ' + re + ex
else begin
ads := ads + ' ' + no + ' ' + re + ex;
if i = 15 then memo1.Lines.Text := memo1.Lines.Text + ads
else memo1.Lines.Append(ads);
end;
end;
end;
end.
複素数演算のプログラム
複素数でベッセル関数の計算をする為には、Delphi標準のVariantによる複素数では、有効桁数が不足するので、多倍長演算を使用します。
多倍長の組み込みは MPArithからmparith_2018-11-27.zipをダウンロードして下さい。
Zipファイルを解凍して、解凍したフォルダーにパス(ツール->オプション->Delphiオプション->ライブラリ->ライブラリパス)を追加し、Uses に mp_types, mp_real, mp_cmplx,
mp_baseを追加すれば使用可能になります。
また、今回はMPArithだけではなく、Rudy's Delphi CornerのBigintegerを使用します。
此処には、BigDecimalもあります。
特徴として、四則演算 + - * / 等がそのまま利用できることです。
但し、三角関数や対数関数はありません。
Rudy's Delphi Cornerを開いたら -> Free Coad -> Bigintegers for Delphi 又は BigDecimals for Delph -> 上の行の最後のリンク DelphiBigNumbers project on GitHub ->Coad▼-> Download ZIP でダウンロードします。
DelphiBigNumbers-master.zip がダウンロード出来たら、解凍したSorce フォルダーを、適当な場所にコピーして、そこへパス(ツール->オプション->Delphiオプション->ライブラリ->ライブラリパス)を追加し、usesに Velthuis.BigDecimals Velthuis.BigIntegers を追加すれば、BigDecimal Biginteger が使用可能となります。
使用方法は、解凍されたPDFファイルを参照してください。
Bigintegerは、ベルヌーイ数の計算にしようします、ベルヌーイ数はガンマ関数の計算に必要なのですが、固定値の配列として与えるのが大変なので、プログラム内で作成します。
Bigintegerと、mp_floatはデーター形式が違うので、テキスト形式で値の受け渡しをします。
Bigintegerは整数形式で、mp_floatの方が有効桁数が低い場合、mp_floatは、指数形式で、有効桁数に丸めて読み込むことが出来ます。
此処では、Bigintegerからmp_floatへの変換しか行いません。
* 重要
複素数の計算にべき乗があるのですが、Delphiに用意されている VerComplexPower、
Mp_complex に用意されている、mpc_powは、そのまま使用するには問題があることがわかりました。
説明には、a^b =
exp(b*ln(a))となっているので間違いないのですが、a と b
が複素数の時、Imaginary部が両方ともゼロ時で aのreal部がマイナスの時は答えのreal部はゼロにならなければなりませんが、意味不明な値が入ります。
もう一つの問題は、乗数 b の値に、**.5 の様に
．５
の値になった時は、整数部の乗数*√の計算なのですがこの時は、realかImaginaryに意味不明な値が値が入ります。
c=a^b
の複素数計算の時、aとb両方の虚数部に0が発生する時は別途正しい答えがでるルーチンを作成する必要があります。
グラフは、全ての点を演算すると時間が掛かるので、十数点を計算し、中間の値は、スプライン補間で表示します。
又、階乗、ディガンマ関数は、プログラムの起動時に、配列データーを作成、ガンマ関数は、ベッセル関数の次数が入力された後、関数の配列を作成、グラフ計算の時間を短縮しています。
有効桁数50桁を計算する為、有効桁数180桁に設定しています。
条件によって、必要な有効桁数が変化するのですが、次数±100
Xの値±100程度の範囲の有効桁数50桁を目標としています。
プログラム
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.ExtCtrls,
VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, VCLTee.TeeProcs,
VCLTee.Chart, System.Diagnostics;
type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
LabeledEdit1: TLabeledEdit;
LabeledEdit2: TLabeledEdit;
Chart1: TChart;
Series1: TLineSeries;
Series2: TLineSeries;
Series3: TPointSeries;
Series4: TPointSeries;
Series5: TLineSeries;
Series6: TLineSeries;
CheckBox1: TCheckBox;
CheckBox2: TCheckBox;
LabeledEdit3: TLabeledEdit;
CheckBox3: TCheckBox;
procedure Button1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure FormClose(Sender: TObject; var Action: TCloseAction);
private
{ Private 宣言 }
procedure chart1coment;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
uses system.Math, System.VarCmplx, Velthuis.BigIntegers,mp_real, mp_cmplx, mp_types, mp_base;
{$R *.dfm}
type
DArray = array of double;
var
xt, yt : Darray; // x,y データー akima補間用
m, t : array of mp_float; // m,t akima補間用
BM : array of mp_float; // ベルヌーイ数配列
FA : array of mp_float; // 階乗配列
PVG : array of mp_float; // +VΓ
MVG : array of mp_float; // +VΓ
zero, one, two, four : mp_float;
three, pai, log_2pis2 : mp_float;
const
KMmax = 250; // KM max 250
Vmax = 100; // v 次数 max
GL = 9; // グラフ計算基本点数
//------------------------------------------------------------------------------
NB = 120; // ベルヌーイ数 配列数 NB + 1
var
NumeratorString : array[0..NB] of string; // 分子
DenominatorString : array[0..NB] of string; // 分母
// 最大公約数 ユークリッドの互除法 BigInteger
function gcd_BigInteger(x, y: BigInteger): BigInteger;
var
t : BigInteger;
begin
while y <> 0 do begin
t := x mod y;
x := y;
y := t;
end;
result := x;
end;
// ベルヌーイ数
// Akiyama-Tanigawa algorithm
// BigInteger
procedure Bernoulli_number_BigInteger;
const
n = (NB + 1) * 2;
var
m, j, k : integer;
a : array of BigInteger; // 分子
b : array of BigInteger; // 分母
tmpN : BigInteger; // 分子
tmpD : BigInteger; // 分母
gcd : BigInteger; // 最大公約数
b0 : BigInteger;
begin
setlength(a, n + 1);
setlength(b, n + 1);
k := 0;
for m := 0 to n do begin
a[m] := 1; // 分子 初期値
b[m] := (m + 1); // 分母 初期値
for j := m - 1 downto 0 do begin
tmpN := a[j] * b[j + 1] - a[j + 1] * b[j]; // 分子
tmpN := tmpN * (j + 1); // 〃
tmpD := b[j] * b[j + 1]; // 分母
gcd := gcd_BigInteger(tmpN, tmpD); // 最大公約数
a[j] := tmpN div gcd;
b[j] := tmpD div gcd;
end;
if (m > 0) and (m mod 2 = 0) then begin
b0 := b[0];
b0 := b0 * m * (m -1); // m ベルヌーイ数No
NumeratorString[k] := a[0].tostring;
DenominatorString[k] := b0.tostring;
inc(k);
end;
end;
end;
//------------------------------------------------------------------------------
// ログガンマ多倍長
procedure log_GammaMul(var x, ans : mp_float);
var
v, w : mp_float;
tmp, tmp0, s : mp_float;
i : integer;
begin
mpf_init2(v, w);
mpf_init3(tmp, tmp0, s);
mpf_set1(v);
mpf_set_int(tmp, NB);
while mpf_is_lt(x, tmp) do begin
mpf_mul(v, x, v);
mpf_add(x, one, x);
end;
mpf_mul(x, x, tmp); // x^2
mpf_div(one, tmp, w); // w = 1 / x^2
mpf_set0(s);
for i := NB downto 1 do begin
mpf_add(s, BM[i], tmp); // tmp = s + B[i]
mpf_mul(tmp, w, s); // s = tmp * w
end;
mpf_add(s, BM[0], tmp); // tmp = s + B[0]
mpf_div(tmp, x, s); // s = (s + B[0]) / x
mpf_add(s, log_2pis2, s); // s = s + ln(2π)/2
mpf_ln(v, tmp); // ln(v)
mpf_sub(s, tmp, s); // s := s - ln(v)
mpf_sub(s, x, s); // s := s - x
mpf_div(one, two, tmp); // tmp = 1/2
mpf_sub(x, tmp, tmp0); // tmp0 = x - 1/2
mpf_ln(x, tmp); // ln(x)
mpf_mul(tmp0, tmp, tmp0); // tmp0 = (x - 1/2) * ln(x)
mpf_add(s, tmp0, ans); // ans = s + (x - 1/2) * ln(x)
mpf_clear2(v, w);
mpf_clear3(tmp, tmp0, s);
end;
// 多倍長ガンマ
// xの値が 0 と負整数の時Γは∞になるのですが処理をしていませんのでエラーになります。
// ケルビンの次数が整数の時は使用されません。
procedure gammaMul(var x, ans: mp_float);
var
tmp, tmp0, logG : mp_float;
begin
mpf_init3(tmp, tmp0, logG);
if mpf_is_lt(x, zero) then begin
mpf_mul(pai, x, tmp); // x*π
mpf_sin(tmp, tmp0); // sin(πx);
mpf_sub(one, x, tmp); // 1-x
log_GammaMul(tmp, logG); // loggamma(1-x);
mpf_exp(logG, logG); // exp(logG)
mpf_div(pai, tmp0, tmp); // π/sin(πx)
mpf_div(tmp, logG, ans); // ans = π/(sin(πx) * logG(1-x))
end
else begin
log_GammaMul(x, logG); // logG
mpf_exp(logG, ans); // exp(logG)
end;
mpf_clear3(tmp, tmp0, logG);
end;
//------------------------------------------------------------------------------
// 階乗 多倍長
procedure factorialMul(n : integer; var ans: mp_float);
label
EXT;
var
i : integer;
bi : mp_float;
begin
mpf_init(bi);
mpf_set1(ans);
mpf_copy(two, bi);
if n <= 1 then begin
goto EXT;
end;
for i := 2 to n do begin
mpf_mul(ans, bi, ans);
mpf_add(bi, one, bi);
end;
EXT:
mpf_clear(bi);
end;
//------------------------------- complex power --------------------------------
// mpc_powは**.5時正しい答えを返しません
procedure mpc_powa(x, y : mp_complex; var ans : mp_complex);
var
ay, ai, zero, harf, two : mp_float;
begin
mpf_init5(ay, ai, zero, harf, two);
mpc_pow(x, y, ans); // exp(y*ln(x))
// xの虚数部とyの虚数部0なら
if mpf_is0(x.im) and mpf_is0(y.im) then begin
mpf_set0(zero); // 0
mpf_set_int(two, 2); // 2
mpf_inv(two, harf); // 1/2 0.5
mpf_abs(y.re, ay); // yの実数部の絶対値
mpf_frac(ay, ai); // yの実数部の絶対値の小数部
mpf_sub(ai, harf, ai); // yの実数部の絶対値の小数部と0.5の差分
// 差分がゼロ(yの実数部の小数部が0.5)
if mpf_is0(ai) then
// xの実数部が0と等しいか0より大きいなら
if mpf_is_ge(x.re, zero) then mpf_set0(ans.im) // 答えの虚数部0
// 0より小さいなら
else mpf_set0(ans.re); // 答えの実数部0
mpf_frac(y.re, ai); // 乗数の小数部
if mpf_is0(ai) then mpf_set0(ans.im); // 整数なら答えの虚数部0
end;
mpf_clear5(ay, ai, zero, harf, two);
end;
//------------------------------------------------------------------------------
// jv(x) 多倍長
// X 値 複素数
// v 次数
procedure Jvx(var v : mp_float; var x, ri: mp_complex);
var
k : integer;
s, x24k, tmp, tmp0, tmp1 : mp_complex;
xc, xs2: mp_complex;
kd, nd, vk : mp_float;
khg, kf : mp_float;
tmf, tmf0 : mp_float;
sb : mp_complex;
begin
mpc_init5(s, x24k, tmp, tmp0, tmp1);
mpc_init2(xc, xs2);
mpf_init2(khg, kf);
mpf_init3(kd, nd, vk);
mpf_init2(tmf, tmf0);
mpc_init(sb);
mpc_set0(s); // Σ=0
mpc_set0(sb);
mpc_set1(x24k);
mpc_div_mpf(x, two, xs2); // x / 2
mpc_mul(xs2, xs2, xc); // (x^2)/ 4
for k := 0 to KMmax do begin
mpf_set1(nd);
mpf_set_int(kf, k); // kf = k
mpf_add(v, kf, tmf0); // k + v
mpf_add(tmf0, one, vk); // vk = k + v + 1;
if mpf_is_lt(vk, zero) then begin // vk < 0 時 nxが整数か確認
mpf_int(vk, tmf); // int(vk);
mpf_sub(vk, tmf, nd); // vk - int(vk) vkが負の整数だったら vk = 0
end;
if mpf_is_ne(vk, zero) and mpf_is_ne(nd, zero) then begin // vkが負の整数の時は計算しない
if mpf_is_ge(v, zero) then
mpf_copy(PVG[k], tmf)
else
mpf_copy(MVG[k], tmf);
// GammaMul(vk, tmf); // Γ(n+k+1)
mpf_mul(FA[k], tmf, khg); // k!Γ(n+k+1) vkが0,負の整数の時±∞
// mpc_set_mpf(kc, kf, zero); // kc = k + 0i 複素数
// mpc_pow(xc, kc, x24k); // (i(x^2)/4)^k
mpc_div_mpf(x24k, khg, tmp0); // ((i(x^2)/4)^k)/(k!Γ(v+k+1))
if Form1.CheckBox3.Checked = false then
if k mod 2 <> 0 then mpc_chs(tmp0, tmp0); // (-1)^k
mpc_add(s, tmp0, s); // Σ
mpc_sub(s, sb, tmp0);
if mpc_is0(tmp0) then break;
mpc_copy(s, sb);
end;
mpc_mul(x24k, xc, x24k);
end;
mpc_set_mpf(tmp, v, zero); // v + 0i
// x が0で次数vが負数の時power演算ゼロでの除算防止
if mpc_is0(x) and mpf_is_lt(v, zero) then // V<0 x=0 時は無限大になるので計算しない
else begin
mpc_powa(xs2, tmp, tmp0); // (x/2)^v
mpc_mul(s, tmp0, s); // (x/2)^v * Σ
end;
mpc_copy(s, ri);
mpc_clear5(s, x24k, tmp, tmp0, tmp1);
mpc_clear2(xc, xs2);
mpf_clear2(khg, kf);
mpf_clear3(kd, nd, vk);
mpf_clear2(tmf, tmf0);
mpc_clear(sb);
end;
// akima m,t テーブル作成
// x[] xデーター配列
// y[] yデーター配列
// m[] m係数テーブル
// t[] t係数テーブル
procedure akima_table;
var
ii, n : integer;
a, b, half, tmf, tmf0 : mp_float;
ytm, xtm : array of mp_float;
tmf1 : mp_float;
begin
n := high(xt) + 1;
setlength(ytm, n);
setlength(xtm, n);
for ii := 0 to n - 1 do begin
mpf_init(ytm[ii]);
mpf_init(xtm[ii]);
end;
mpf_init5(a, b, half, tmf, tmf0);
mpf_init(tmf1);
mpf_set_dbl(half, 1 / 2);
// mpf_set_int(tow, 2);
// mpf_set0(zero);
for ii := 0 to n -1 do begin
mpf_set_dbl(xtm[ii], xt[ii]);
mpf_set_dbl(ytm[ii], yt[ii]);
end;
// shift data by + 2 in the array and compute the secants
// also calculate extrapolated and end point secants
// 傾斜α両端を除く Δy/Δx
for ii := 0 to n - 2 do begin
mpf_sub(ytm[ii + 1], ytm[ii], tmf);
mpf_sub(xtm[ii + 1], xtm[ii], tmf0);
mpf_div(tmf, tmf0, m[ii + 2]);
end;
// for ii := 0 to n - 2 do
// m[ii + 2] := (yt[ii + 1] - yt[ii]) / (xt[ii + 1] - xt[ii]);
// 端点傾斜処理
mpf_mul(two, m[2], tmf);
mpf_sub(tmf, m[3], m[1]);
// m[1] := 2 * m[2] - m[3];
mpf_mul(two, m[1], tmf);
mpf_sub(tmf, m[2], m[0]);
// m[0] := 2 * m[1] - m[2];
mpf_mul(two, m[n], tmf);
mpf_sub(tmf, m[n - 1], m[n + 1]);
// m[n + 1] := 2 * m[n] - m[n - 1];
mpf_mul(two, m[n + 1], tmf);
mpf_sub(tmf, m[n], m[n + 2]);
// m[n + 2] := 2 * m[n + 1] - m[n];
// 各ポイントの傾斜係数計算
for ii := 0 to n - 1 do begin
mpf_sub(m[ii + 3],m[ii + 2],tmf0);
mpf_abs(tmf0, a);
// a := abs(m[ii + 3] - m[ii + 2]);
mpf_sub(m[ii + 1], m[ii], tmf0);
mpf_abs(tmf0, b);
// b := abs(m[ii + 1] - m[ii]);
mpf_add(a, b, tmf1);
if mpf_is_ne(tmf1, zero) then begin
mpf_mul(a, m[ii + 1], tmf);
mpf_mul(b, m[ii + 2], tmf0);
mpf_add(tmf, tmf0, tmf);
mpf_div(tmf, tmf1, t[ii]);
end
else begin
mpf_add(m[ii + 2], m[ii + 1], tmf);
mpf_mul(half, tmf, t[ii]);
end;
{
if (a + b) <> 0 then begin
t[ii] := (a * m[ii + 1] + b * m[ii + 2]) / (a + b);
end
else
t[ii] := half * (m[ii + 2] + m[ii + 1]);
}
end;
for ii := 0 to n - 1 do begin
mpf_clear(ytm[ii]);
mpf_clear(xtm[ii]);
end;
mpf_clear5(a, b, half, tmf, tmf0);
mpf_clear(tmf1);
end;
// akima 補間値計算
// xx xの値
// x[] xデーター配列
// y[] yデーター配列
// m[] m係数テーブル
// t[] t係数テーブル
// result 補間値y'
function akima_Interpolation(xx: double): double;
var
iB, iM, iT: integer;
a, b, tmf, tmf0 : mp_float;
c, d, e, f, tmf1 : mp_float;
three : mp_float;
begin
mpf_init4(a, b, tmf, tmf0);
mpf_init5(c, d, e, f, tmf1);
mpf_init(three);
mpf_set_int(three, 3);
iB := low(xt); // x[] bottom 配列no
iT := high(xt); // x[] top配列No
// xx値の上下の配列xの配列番号を探す
// XX<x[iB]の場合一番下の区間 XX>x[iT]の場合一番上の区間
while (iT - iB) > 1 do begin
iM := (iB + iT) div 2; // middle配列no
if xt[iM] > xx then
iT := iM
else
iB := iM;
end;
mpf_set_dbl(b, xt[iT] - xt[iB]);
// b := xt[iT] - xt[iB]; // 区間のxの変化量
mpf_set_dbl(a, xx - xt[iB]);
// a := xx - xt[iB]; // x[iB]からのxの値
// 3次akima spline 計算
mpf_set_dbl(c, yt[iB]); // c = yt[ib]
mpf_mul(t[iB], a, d); // d = t[ib] * a
mpf_mul(three, m[iB + 2], tmf); // 3 * m[iB + 2]
mpf_mul(two, t[ib], tmf0); // 2 * t[iB]
mpf_sub(tmf, tmf0, tmf1); // 3 * m[iB + 2] - 2 * t[iB]
mpf_sub(tmf1, t[iB + 1], tmf); // 3 * m[iB + 2] - 2 * t[iB] - t[iB + 1]
mpf_mul(tmf, a, tmf); // (3 * m[iB + 2] - 2 * t[iB] - t[iB + 1]) * a
mpf_mul(tmf, a, tmf); // (3 * m[iB + 2] - 2 * t[iB] - t[iB + 1]) * a * a
mpf_div(tmf, b, e); // (3 * m[iB + 2] - 2 * t[iB] - t[iB + 1]) * a * a / b
mpf_add(t[iB], t[iB + 1], tmf); // t[iB] + t[iB + 1]
mpf_mul(two, m[iB + 2], tmf0); // 2 * m[iB + 2]
mpf_sub(tmf, tmf0, tmf); // t[iB] + t[iB + 1] - 2 * m[iB + 2]
mpf_mul(tmf, a, tmf); // (t[iB] + t[iB + 1] - 2 * m[iB + 2]) * a
mpf_mul(tmf, a, tmf); // (t[iB] + t[iB + 1] - 2 * m[iB + 2]) * a * a
mpf_mul(tmf, a, tmf); // (t[iB] + t[iB + 1] - 2 * m[iB + 2]) * a * a * a
mpf_mul(b, b, tmf0); // b * b
mpf_div(tmf, tmf0, f); // (t[iB] + t[iB + 1] - 2 * m[iB + 2]) * a * a * a / (b * b)
mpf_add(c, d, tmf); // c + d
mpf_add(tmf, e, tmf); // c + d + e
mpf_add(tmf, f, tmf); // c + d + e + f
result := mpf_todouble(tmf);
{
result := yt[iB]
+ t[iB] * a
+ (3 * m[iB + 2] - 2 * t[iB] - t[iB + 1]) * a * a / b
+ (t[iB] + t[iB + 1] - 2 * m[iB + 2]) * a * a * a / (b * b);
}
mpf_clear4(a, b, tmf, tmf0);
mpf_clear5(c, d, e, f, tmf1);
mpf_clear(three);
end;
procedure TForm1.chart1coment;
begin
Chart1.Canvas.Font.Size := 8;
Chart1.Canvas.Font.Style := [];
chart1.Canvas.Font.Color := clred;
chart1.Canvas.TextOut(chart1.Width - 30, 1,'実数');
chart1.Canvas.Font.Color := clblue;
chart1.Canvas.TextOut(chart1.Width - 30, 13,'虚数');
chart1.Canvas.Pen.Width := 1;
chart1.Canvas.Pen.Color := clred;
chart1.Canvas.MoveTo(chart1.Width - 50, 6);
chart1.Canvas.LineTo(chart1.Width - 32, 6);
chart1.Canvas.Pen.Color := clblue;
chart1.Canvas.MoveTo(chart1.Width - 50, 18);
chart1.Canvas.LineTo(chart1.Width - 32, 18);
end;
// 計算
// xの値が大きくなると誤差が大きくなります。
// ta[]グラフ用計算点
procedure TForm1.Button1Click(Sender: TObject);
label
EXT;
const
x0m = 1e-50; // ゼロ近傍値 infinty の符号設定計算値
ta : array[0..GL - 1] of double = (0.08, 0.14, 0.23, 0.46, 0.86, 1.5, 2.5, 3.2, 4);
YPM = 1E304; // Double 演算オーバーフロー対策 最大値制限値
YMM = - YPM;
zeromg = '1E-90'; // 計算結果表示これより小さい値は0にします 演算桁数の1/2
var
ch, i, xi: integer;
berx, beix: double;
berxe, beixe: double;
xin, vin, xv, rk, ik: double;
xmin, xmax, dx, dxf : double;
ymin, ymax: double;
berl : Darray;
beil : Darray;
beru : Darray;
beiu : Darray;
xl : Darray;
xu : Darray;
xm, vm, xvm, ixm: mp_float;
nd, tmf : mp_float;
ri : mp_complex;
GCF : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
mm, ss, ms : integer;
xc, xcb : mp_complex;
GU, GS : integer; // xv = 0 がない場合 推定点数
k : integer;
vk, avm, tmf0, zz : mp_float;
ix : double;
istr : string;
// double to mpc グラフ計算用
procedure xvtoXc(xv, ia : double; var xc : mp_complex);
begin
mpf_set_dbl(xvm, xv);
mpf_set_dbl(tmf, ia);
mpc_set_mpf(xc, xvm, tmf);
end;
function maxmin(x : double): double;
begin
result := x;
if x > YPM then result := YPM;
if x < YMM then result := YMM;
end;
begin
memo1.Clear;
val(labelededit1.Text, vin, ch);
if ch <> 0 then begin
application.MessageBox('次数vの値に間違いがあります。','注意', 0);
exit;
end;
if abs(vin) > Vmax then begin
application.MessageBox('次数vの値が計算範囲外です。','注意', 0);
exit;
end;
val(labelededit2.Text, xin, ch);
if ch <> 0 then begin
application.MessageBox('xの値に間違いがあります。','注意', 0);
exit;
end;
if abs(xin) > 100 then begin
application.MessageBox('xが100を越えると条件によって誤差が大きくなります。','注意', 0);
end;
val(labelededit3.Text, ix, ch);
if ch <> 0 then begin
application.MessageBox('i の値に間違いがあります。','注意', 0);
exit;
end;
if abs(ix) > 100 then begin
application.MessageBox('iが100を越えると条件によって誤差が大きくなります。','注意', 0);
end;
mpf_init4(xm, vm, xvm, ixm);
mpf_init3(nd, tmf, zz);
mpf_init3(vk, avm, tmf0);
mpc_init3(ri, xc, xcb);
mpf_set0(zero);
mpf_set0(ixm);
mpf_read_decimal(vm, PAnsiChar(ansistring(labelededit1.Text + #00)));
mpf_read_decimal(xm, PAnsiChar(ansistring(labelededit2.Text + #00)));
mpf_read_decimal(ixm, PAnsiChar(ansistring(labelededit3.Text + #00)));
mpf_read_decimal(zz, PAnsiChar(ansistring(zeromg + #00)));
mpc_set_mpf(xcb, xm, ixm); // xcb 計算用 xの複素数
series1.Clear;
series2.Clear;
series3.Clear;
series4.Clear;
series5.Clear;
series6.Clear;
if ix >= 0 then istr := ' +' + floatTostr(ix) + 'i'
else istr := ' ' + floatTostr(ix) + 'i';
memo1.Lines.Append('v = ' + floatTostr(vin) + ' x = ' + floatTostr(xin) + istr);
application.ProcessMessages;
rk := 0;
ik := 0;
mpf_abs(vm, avm);
for k := 0 to KMmax do begin
mpf_set_int(tmf, k); // k
mpf_add(tmf, avm, tmf0); // v + k
mpf_add(tmf0, one, vk); // vk= v + k + 1
gammaMul(vk, PVG[k]); // Γ(n+k+1)
end;
mpf_chs(avm, avm); // -v
for k := 0 to KMmax do begin
mpf_set1(nd);
mpf_set_int(tmf, k); // k
mpf_add(tmf, avm, tmf0); // -v + k
mpf_add(tmf0, one, vk); // vk= -v + k + 1
if mpf_is_lt(vk, zero) then begin // vk < 0 時 nxが整数か確認
mpf_int(vk, tmf); // int(vk);
mpf_sub(vk, tmf, nd); // vk - int(vk) vkが負の整数だったら nd = 0
end;
if mpf_is_ne(vk, zero) and mpf_is_ne(nd, zero) then // vkが負の整数の時は計算しない
gammaMul(vk, MVG[k]) // Γ(n+k+1)
else
mpf_set0(MVG[k]); // Γ(n+k+1) = 0
end;
// ゼロ近傍の計算符号設定
// x=0 v < 0 v = 非整数 の場合 ±∞
if mpc_is0(xcb) and mpf_is_lt(vm, zero) and mpf_is_ne(nd, zero) then begin
xvtoxc(x0m, ix, xc); // x double to mpc
Jvx(vm, xc, ri); // ゼロ近辺のプラス側x計算
if mpf_is_gt(ri.re, zero) then rk := infinity; // x=0時の無限大±符号設定
if mpf_is_lt(ri.re, zero) then rk := -infinity;
if mpf_is_gt(ri.im, zero) then ik := infinity;
if mpf_is_lt(ri.im, zero) then ik := -infinity;
end;
mpf_int(vm, tmf); // int(vm)
mpf_sub(vm, tmf, nd); // nd = vm - int(v)
StopWatch := TStopwatch.StartNew;
// 表示値の計算
Jvx(vm, xcb, ri); // xcb 複素数
// x=0 v < 0 v = 非整数 の場合 ±∞
if mpc_is0(xcb) and mpf_is_lt(vm, zero) and mpf_is_ne(nd, zero) then begin
berxe := rk;
beixe := ik;
if checkbox3.Checked = false then
memo1.Lines.Append('Jv(x) = ' + floatTostr(berxe)) // 数値表示
else
memo1.Lines.Append('Iv(x) = ' + floatTostr(berxe)); // 数値表示
if beixe >= 0 then
memo1.Lines.Append(' +' + floatTostr(beixe) + 'i');
end
else begin
if checkbox3.Checked = false then
memo1.Lines.Append( string('Jv(x) = ' + mpf_decimal(ri.re, 50)))
else
memo1.Lines.Append( string('Iv(x) = ' + mpf_decimal(ri.re, 50)));
if mpf_is_ge(ri.im, zero) then
memo1.Lines.Append( string(' +' + mpf_decimal(ri.im, 50) + 'i'))
else
memo1.Lines.Append( string(' ' + mpf_decimal(ri.im, 50) + 'i'));
berxe := mpf_todouble(ri.re);
beixe := mpf_todouble(ri.im);
end;
StopWatch.Stop;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
ms := ElapsedMillseconds * (GL + 1) * 2 + 1000;
mm := ms div 60000;
ss := (ms mod 60000) div 1000;
memo1.Lines.Append('グラフ表示に約' + inttostr(mm) + '分' + inttostr(ss) + '秒必要です。');
if checkbox2.Checked = true then begin // グラフ無だったら終了
Chart1.Canvas.Font.Style := [fsBold];
Chart1.Canvas.Font.size := 8;
Chart1.Canvas.TextOut(170, 115,'グラフ無し');
goto EXT;
end;
// 最大値最小値の検索とグラフデーター作成
xi := round(xin);
xmin := xi - 4;
GCF := 0;
if (Xin >= -2) and (xin <= 2) then GCF := 1;
if (xin >= -3) and (xin <= -2) then GCF := 2;
if (xin >= 2) and (xin <= 3) then GCF := 3;
if (xin >= -4) and (xin <= -3) then GCF := 4;
if (xin >= 3) and (xin <= 4) then Gcf := 5;
case GCF of
0: xmin := xi - 4;
1: xmin := -4;
2: xmin := -5;
3: xmin := -3;
4: xmin := -6;
5: xmin := -2;
end;
xmax := xmin + 8;
case GCF of
0 : begin
setlength(berl, GL + GL);
setlength(beil, GL + GL);
setlength(xt, GL + GL); setlength(yt, GL + GL);
setlength(xl, GL + GL);
end;
1, 2, 3, 4, 5:
begin
setlength(berl, GL); setlength(beru, GL);
setlength(beil, GL); setlength(beiu, GL);
setlength(xt, GL); setlength(yt, GL);
setlength(xl, GL); setlength(xu, GL);
end;
end;
dxf := 0.1;
dx := 0.1;
GU := 40;
GS := 40;
case GCF of
0 :
begin
dx := (xmax - xmin) / (Gl + GL - 1);
for i := 0 to GL + GL - 1 do xl[i] := dx * i + xmin;
xl[0] := xl[0] + dxf;
xl[GL + GL - 1] := xl[GL + GL - 1] - dxf;
GS := 80;
end;
else
begin
for i := 0 to GL - 1 do begin
xu[i] := ta[i];
xl[GL - i - 1] := -ta[i];
end;
end;
end;
case GCF of
2 : begin dx := 5 / 4; dxf := 3 / 4; GS := 50; GU := 30; end;
3 : begin dx := 3 / 4; dxf := 5 / 4; GS := 30; GU := 50; end;
4 : begin dx := 6 / 4; dxf := 2 / 4; GS := 60; GU := 20; end;
5 : begin dx := 2 / 4; dxf := 6 / 4; GS := 20; GU := 60; end;
end;
if GCF >= 2 then begin
for i := 0 to Gl - 1 do xl[i] := xl[i] * dx;
for i := 0 to Gl - 1 do xu[i] := xu[i] * dxf;
end;
// グラフ用データー作成
ymin := 0;
ymax := 0;
case GCF of
0 : begin
for i := 0 to GL + GL - 1 do begin
xv := xl[i];
xvtoxc(xv, ix, xc);
Jvx(vm, xc, ri);
berl[i] := maxmin(mpf_todouble(ri.re));
beil[i] := maxmin(mpf_todouble(ri.im));
if berl[i] > ymax then ymax := berl[i];
if beil[i] > ymax then ymax := beil[i];
if berl[i] < ymin then ymin := berl[i];
if beil[i] < ymin then ymin := beil[i];
end;
end;
1, 2, 3, 4, 5:
for i := 0 to GL - 1 do begin
xv := xl[i];
xvtoxc(xv, ix, xc);
Jvx(vm, xc, ri);
berl[i] := maxmin(mpf_todouble(ri.re));
beil[i] := maxmin(mpf_todouble(ri.im));
if berl[i] > ymax then ymax := berl[i];
if beil[i] > ymax then ymax := beil[i];
if berl[i] < ymin then ymin := berl[i];
if beil[i] < ymin then ymin := beil[i];
xv := xu[i];
xvtoxc(xv, ix, xc);
Jvx(vm, xc, ri);
beru[i] := maxmin(mpf_todouble(ri.re));
beiu[i] := maxmin(mpf_todouble(ri.im));
if beru[i] > ymax then ymax := beru[i];
if beiu[i] > ymax then ymax := beiu[i];
if beru[i] < ymin then ymin := beru[i];
if beiu[i] < ymin then ymin := beiu[i];
end;
end;
// 指定値の値制御
if berxe > ymax then berxe := ymax;
if berxe < ymin then berxe := ymin;
if beixe > ymax then beixe := ymax;
if beixe < ymin then beixe := ymin;
series3.AddXY(xin, berxe);
series4.AddXY(xin, beixe);
if checkbox1.Checked = true then
case GCF of
0:
begin
for i := 0 to GL + GL - 1 do series1.AddXY(xl[i], berl[i]);
for i := 0 to GL + GL - 1 do series2.AddXY(xl[i], beil[i]);
end;
1, 2, 3, 4, 5 :
begin
for i := 0 to GL - 1 do series1.AddXY(xl[i], berl[i]);
for i := 0 to GL - 1 do series2.AddXY(xl[i], beil[i]);
for i := 0 to GL - 1 do series5.AddXY(xu[i], beru[i]);
for i := 0 to GL - 1 do series6.AddXY(xu[i], beiu[i]);
end;
end;
if checkbox1.Checked = false then
// グラフ計算
case GCF of
0:
begin
for i := 0 to GL + GL - 1 do begin
xt[i] := xl[i];
yt[i] := berl[i];
end;
akima_table;
dx := (xt[GL + GL - 1] - xt[0]) / GS;
for i := 0 to GS do begin
xv := dx * i + xt[0];
berx := maxmin(akima_Interpolation(xv));
series1.AddXY(xv, berx);
end;
for i := 0 to GL + GL - 1 do yt[i] := beil[i];
akima_table;
for i := 0 to GS do begin
xv := dx * i + xt[0];
beix := maxmin(akima_Interpolation(xv));
series2.AddXY(xv, beix);
end;
end;
1, 2, 3, 4, 5 :
begin
for i := 0 to GL- 1 do begin
xt[i] := xl[i];
yt[i] := berl[i];
end;
akima_table;
dx := xt[0] / GS;
for i := GS downto 1 do begin
xv := dx * i;
berx := maxmin(akima_Interpolation(xv));
series1.AddXY(xv, berx);
end;
for i := 0 to GL- 1 do yt[i] := beil[i];
akima_table;
for i := GS downto 1 do begin
xv := dx * i;
beix := maxmin(akima_Interpolation(xv));
series2.AddXY(xv, beix);
end;
if mpf_is0(xm) then begin
// series1.AddXY(0, herxe);
// series2.AddXY(0, heixe);
series5.AddXY(0, berxe);
series6.AddXY(0, beixe);
end;
for i := 0 to GL- 1 do begin
xt[i] := xu[i];
yt[i] := beru[i];
end;
akima_table;
dx := xt[GL - 1] / GU;
for i := 1 to GU do begin
xv := dx * i;
berx := maxmin(akima_Interpolation(xv));
series5.AddXY(xv, berx);
end;
for i := 0 to GL- 1 do yt[i] := beiu[i];
akima_table;
for i := 1 to GU do begin
xv := dx * i;
beix := maxmin(akima_Interpolation(xv));
series6.AddXY(xv, beix);
end;
end;
end;
application.ProcessMessages;
chart1coment;
EXT:
mpf_clear4(xm, vm, xvm, ixm);
mpf_clear3(nd, tmf, zz);
mpc_clear3(ri, xc, xcb);
mpf_clear3(vk, avm, tmf0);
end;
procedure TForm1.FormCreate(Sender: TObject);
var
i : integer;
N, D, tmp : mp_float;
begin
mpf_set_default_decprec(180); // 有効桁数180桁 50桁の精度に必要です。
setlength(BM, NB + 1); // ベルヌーイ数配列
setlength(FA, KMmax + 1); // K1配列
setlength(PVG, KMmax + 1); // +vΓ
setlength(MVG, KMmax + 1); // -vΓ
setlength(t, GL + GL); // akima 補間値計算 配列 t
setlength(m, GL + GL + 3); // akima 補間値計算 配列 m
for i := 0 to NB do mpf_init(BM[i]);
for i := 0 to KMmax do mpf_init(FA[i]);
for i := 0 to KMmax do mpf_init(PVG[i]);
for i := 0 to KMmax do mpf_init(MVG[i]);
for i := 0 to GL + GL - 1 do mpf_init(t[i]);
for i := 0 to GL + GL + 2 do mpf_init(m[i]);
mpf_init3(N, D, tmp);
mpf_init4(zero, one, two, four);
mpf_init3(three, pai, log_2pis2);
mpf_set0(zero);
mpf_set1(one);
mpf_set_int(two, 2);
mpf_set_int(three, 3);
mpf_set_int(four, 4);
mpf_set_pi(pai);
mpf_mul(pai, two, tmp); // 2π
mpf_ln(tmp, tmp); // ln(2π)
mpf_div(tmp, two, log_2pis2); // ln(2π)/2
Bernoulli_number_BigInteger; // ベルヌーイ数作成
for i := 0 to NB do begin
mpf_read_decimal(N, PAnsiChar(ansistring(NumeratorString[i] + #00)));
mpf_read_decimal(D, PAnsiChar(ansistring(DenominatorString[i] + #00)));
mpf_div(N, D, BM[i]);
end;
for i := 0 to KMmax do begin
factorialMul(i, N);
mpf_copy(N, FA[i]);
end;
memo1.Clear;
mpf_clear3(N, D, tmp);
end;
procedure TForm1.FormClose(Sender: TObject; var Action: TCloseAction);
var
i : integer;
begin
for i := 0 to NB do mpf_clear(BM[i]);
for i := 0 to KMmax do mpf_clear(FA[i]);
for i := 0 to KMmax do mpf_clear(PVG[i]);
for i := 0 to KMmax do mpf_clear(MVG[i]);
for i := 0 to GL + GL - 1 do mpf_clear(t[i]);
for i := 0 to GL + GL + 2 do mpf_clear(m[i]);
mpf_clear4(zero, one, two, four);
mpf_clear3(three, pai, log_2pis2);
end;
end.
Bessel1_function.zip