2024/03/26 計算ルーチンを見直して演算を早くしました。
ｘの値を実数とする第２, ３種 ケルビン関数
第２種ケルビン関数 ここで扱うのはｘの値が実数の場合です。
第２種ケルビン関数は第２種変形ベッセル関数を元に計算されます。
第２種変形ベッセル関数は、次数が非整数と、整数の場合で違う計算を行います。
非整数の計算
整数の計算
第３種ケルビン関数
第３種ケルビン関数は、第２種に 2/(iπ) を乗じているだけなので、第２種のプログラムに追加です。
次数ｖ変数ｘは両方とも実数で、答えは実数と虚数になります。
次数ｖの最大値は±100です。
グラフは、
赤が実数で青が虚数です。。。。。
グラフを作成するのに、全ての点を計算すると、非常に時間が掛かる為、十数点を計算し、ａｋｉｍａスプライン補間を使用して作図していますので、プロットされた点と、ラインが多少ずれる場合があります。
3rd kindのチェックボックスにチェックを入れれば第３種の計算になります。
プログラム
プログラムには、計算精度向上の為、多倍長演算が組み込んであります。
組み込み方は第１種ケルビン関数を参照して下さい。
プログラム中の Σ の計算に収束判定が入っていませんので、追加すれば Σ の演算時間を短くすることが出来ます。
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, system.UITypes;
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;
Check3rd: TCheckBox;
procedure Button1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure FormDestroy(Sender: TObject);
private
{ Private 宣言 }
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; // ベルヌーイ数配列
DG : array of mp_float; // ディガンマ配列
FA : array of mp_float; // 階乗配列
PVG : array of mp_float; // +VΓ
MVG : array of mp_float; // +VΓ
KMV : integer; // Σ計算最大値
gamma : mp_float;
pai, zero, one, two, four : mp_float;
log_2pis2 : mp_float;
const
KMmax = 250; // KM max
Vmax = 100; // v 次数 max
GL = 8; // グラフ分割数
// オイラー定数 有効桁分用意
// 0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504
// 01448654283622417399764492353625350033374293733773767394279259525824709491600873520394816567085323315177661152862119950150798479374
gstr0 = '0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504';
gstr1 = '0144865428362241739976449235362535003337429373377376739427925952582470949160087352039481656708532331517766115286211995015079847';
gstr = gstr0 + gstr1;
//------------------------------------------------------------------------------
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); // 1
mpf_set_int(tmp, NB); // tmp=NB
while mpf_is_lt(x, tmp) do begin
mpf_mul(v, x, v); // v = v * x
mpf_add(x, one, x); // x = x + 1
end;
mpf_mul(x, x, tmp); // x^2
mpf_div(one, tmp, w); // w = 1 / x^2
mpf_set0(s); // s = 0
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 + 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 digammaMUL(n: integer; var ans : mp_float);
var
s, tmp : mp_float;
k : integer;
begin
mpf_init2(s, tmp);
n := n - 1;
mpf_set0(s); // Σ=0
mpf_set0(tmp);
if n >= 0 then begin
for k := 1 to n do begin
mpf_set_int(tmp, k);
mpf_div(one, tmp, tmp); // 1 / k
mpf_add(s, tmp, s);
end;
end;
mpf_sub(s, gamma, ans); // Σ + γ
mpf_clear2(s, tmp);
end;
//------------------------------------------------------------------------------
// 階乗 多倍長
procedure factorialMul(n : integer; var ans: mp_float);
label
EXT;
var
i : integer;
bi: mp_float;
begin
mpf_init(bi);
mpf_set1(ans); // 1
mpf_add(one, one, bi); // bi = 2
if n <= 1 then begin // n <= 1 ans = 1
goto EXT;
end;
for i := 2 to n do begin // n!
mpf_mul(ans, bi, ans); // ans = ans * bi
mpf_add(bi, one, bi); // bi + 1
end;
EXT:
mpf_clear(bi);
end;
//--------------------------------- 非整数計算 -------------------------------
var
rk, ik : double; // x=0 で、無限大時の値
// 第一種変形ベッセル関数
// 非整数多倍長で計算
procedure sigumaaIaxMul(var v : mp_float; var x, ri : mp_complex);
var
k : integer;
khg : mp_float;
s, x24k, tmp, tmp0, xei: mp_complex;
kc, xc : mp_complex;
tmf0 : mp_float;
begin
mpc_init5(s, x24k, tmp, tmp0, xei);
mpc_init2(kc, xc);
mpf_init(khg);
mpf_init(tmf0);
mpc_set0(s); // Σ=0
// mpc_copy(x, xei);
mpc_mul(x, x, xei);
mpc_set1(x24k);
for k := 0 to KMV do begin
if mpf_is_ge(v, zero) then
mpf_copy(PVG[k], tmf0)
else
mpf_copy(MVG[k], tmf0);
// gammaMul(vk, tmf0); // Γ(n+k+1)
mpf_mul(FA[k], tmf0, khg); // k!Γ(n+k+1)
// mpf_set_int(tmf, K + K); // 2k
// mpc_set_mpf(kc, tmf, zero); // kc = 2k + 0i 複素数
// mpc_pow(xei, kc, x24k); // (i(x^2)/4)^k
mpc_div_mpf(x24k, khg, tmp0); // ((i(x^2)/4)^k)/(k!Γ(v+k+1)) vkが0,負の整数の時±∞
mpc_add(s, tmp0, s); // Σ
mpc_mul(x24k, xei, x24k);
end;
if mpf_is_ge(v, zero) or not mpc_is0(x) then begin // v がゼロより大きいか xがゼロでな時
mpc_set_mpf(tmp, v, zero); // v + 0i
mpc_pow(x, tmp, tmp0); // ((xe^(iπ/4)) /2)^V
mpc_mul(s, tmp0, s); // Σ
end;
mpc_copy(s, ri);
mpc_clear5(s, x24k, tmp, tmp0, xei);
mpc_clear2(kc, xc);
mpf_clear(khg);
mpf_clear(tmf0);
end;
// 非整数 多倍長
// 第一種変形ベッセル関数の差分をsinで除算するため演算精度向上には大きな有効桁数が必要になります。
procedure ekazMul(var v : mp_float; var x, ri: mp_complex);
var
sv, sinvp, tmp : mp_float;
ipPv, imPv, sipv, tmpc : mp_complex;
tmp0, tmp1, tmp2 : mp_complex;
begin
mpf_init3(sv, sinvp, tmp);
mpc_init4(ipPv, imPv, sipv, tmpc);
mpc_init3(tmp0, tmp1, tmp2);
sigumaaIaxMul(v, x, ipPv); // ipv
mpf_chs(v, sv);
sigumaaIaxMul(sv, x, imPv); // imv
if not mpc_is0(x) then begin
mpc_sub(imPv, ipPv, sipv); // imv - ipv
mpf_mul(v, pai, tmp); // v * π
mpc_set_mpf(tmpc, tmp, zero);
mpc_sin(tmpc, tmpc); // sin(v * pi)
mpc_div(sipv, tmpc, sipv); // (imv - ipv)/sin(v * pi)
mpf_div(pai, two, tmp); // π /2
mpc_mul_mpf(sipv, tmp, sipv); // (Imv - Ipv) / sin(v * pi) * pi / 2
mpf_mul(sv, tmp, tmp); // -v * π/2
mpc_set_mpf(tmpc, zero, tmp); // 0 + -v * pi / 2
mpc_exp(tmpc, tmpc); // e^(-ivπ/2)
mpc_mul(sipv, tmpc, ri); // ((Imv - Ipv) / sin(v * pi) * pi / 2) * e^(-ivπ/2)
end;
mpf_clear3(sv, sinvp, tmp);
mpc_clear4(ipPv, imPv, sipv, tmpc);
mpc_clear3(tmp0, tmp1, tmp2);
end;
//---------------------------整数計算-------------------------------------------
// 整数 多倍長計算
// 計算精度の向上には、Σの最大値のUPと演算有効桁数のUPとなります。
procedure knzMul(var v : mp_float; var z, ri : mp_complex);
var
a, b, c, inz, zc : mp_complex;
z2, knza, z24, mz24: mp_complex;
tmp, tmp0, tmp1 : mp_complex;
mone : mp_float;
mf, nnf, kf, tmf, tmf0 : mp_float;
tmf1, nb, n : mp_float;
k, m, nn : integer;
begin
if mpc_is0(z) then begin // z = 0 の時は±∞となるので計算しない
mpc_set0(ri);
exit;
end;
mpc_init5(a, b, c, inz, zc);
mpc_init4(z2, knza, z24, mz24);
mpc_init3(tmp, tmp0, tmp1);
mpf_init(mone);
mpf_init4(nnf, kf, tmf, tmf0);
mpf_init3(tmf1, nb, n);
mpf_copy(v, nb);
mpf_abs(v, n);
nn := abs(round(mpf_todouble(n)));
mpf_chs(one, mone);
mpc_copy(z, z2); // z -> z2
mpc_mul(z2, z2, z24); // (ze^(iπ/4) / 2)^2
mpc_chs(z24, mz24); // -(ze^(iπ/4) / 2)^2
// form1.memo1.Lines.Append( string(mpf_decimal_alt(mz24.re, 50)));
mpc_set0(inz);
mpc_set1(tmp1);
// Inz
for m := 0 to KMV do begin
mpf_mul(FA[m], FA[m + nn], tmf);
mpf_inv(tmf, tmf0);
mpc_set_mpf(tmp, tmf0, zero);
// mpf_set_int(mf, m);
// mpc_set_mpf(tmp0, mf, zero);
// mpc_pow(z24, tmp0, tmp1);
mpc_mul(tmp, tmp1, tmp0);
mpc_add(inz, tmp0, inz);
mpc_mul(tmp1, z24, tmp1);
end;
mpf_set_int(nnf, nn);
mpc_set_mpf(tmp0, nnf, zero);
mpc_pow(z2, tmp0, tmp);
mpc_mul(inz, tmp, inz);
// form1.memo1.Lines.Append( string(mpf_decimal_alt(inz.re, 50)));
// a
mpc_set0(a);
mpc_set1(tmp1);
for k := 0 to nn - 1 do begin
mpf_div(FA[nn - k - 1], FA[k], tmf);
mpc_set_mpf(tmp0, tmf, zero);
// mpf_set_int(kf, k);
// mpc_set_mpf(tmp1, kf, zero);
// mpc_pow(mz24, tmp1, tmp1);
mpc_mul(tmp0, tmp1, tmp);
mpc_add(a, tmp, a);
mpc_mul(tmp1, mz24, tmp1);
end;
mpf_set_int(nnf, -nn);
mpc_set_mpf(tmp1, nnf, zero);
mpc_pow(z2, tmp1, tmp);
mpc_mul(a, tmp, tmp0);
mpc_div_mpf(tmp0, two, a);
// form1.memo1.Lines.Append( string(mpf_decimal_alt(a.re, 50)));
// b
if (nn + 1) mod 2 = 0 then
mpc_set_mpf(tmp1, one, zero)
else
mpc_set_mpf(tmp1, mone, zero);
mpc_ln(z2, tmp0);
mpc_mul(tmp1, tmp0, tmp);
mpc_mul(tmp, inz, b);
// form1.memo1.Lines.Append( string(mpf_decimal_alt(b.re, 50)));
// c
mpc_set0(c);
mpc_set1(mz24);
for k:= 0 to KMV do begin
mpf_add(DG[k + 1], DG[nn + k + 1], tmf);
mpc_set_mpf(tmp1, tmf, zero);
mpf_mul(FA[k], FA[nn + k], tmf0);
mpc_set_mpf(tmp0, tmf0, zero);
// mpf_set_int(kf, k);
// mpc_set_mpf(tmp, kf, zero);
// mpc_pow(z24, tmp, tmp);
mpc_mul(tmp, tmp1, tmp);
mpc_div(tmp, tmp0, tmp);
mpc_add(c, tmp, c);
mpc_mul(mz24, z24, mz24);
end;
// form1.memo1.Lines.Append('c '+ string(mpf_decimal_alt(c.re, 50)));
// form1.memo1.Lines.Append('fa '+ string(mpf_decimal_alt(FA[3], 50)));
// form1.memo1.Lines.Append('dg '+ string(mpf_decimal_alt(DG[3], 50)));
// c := c * power(-1, nn) / 2 * varcomplexpower(z2, nn);
mpf_set_int(nnf, nn);
mpc_set_mpf(tmp0, nnf, zero);
mpc_pow(z2, tmp0, tmp1); // z2^n
mpc_mul(c, tmp1, tmp); // c * z2^n
if nn mod 2 = 0 then // -1^nn
mpc_set_mpf(tmp1, one, zero)
else
mpc_set_mpf(tmp1, mone, zero);
mpc_div_mpf(tmp1, two, tmp0); // (-1~n)/2
mpc_mul(tmp, tmp0, c); // c := z2^n * c * (-1~n)/2
// knza := a + b + c;
mpc_add(a, b, tmp);
mpc_add(tmp, c, knza);
// knza * e^ -ivπ/2
mpf_set_int(tmf, -nn); // -nn
mpf_div(pai, two, tmf0); // π/2
mpf_mul(tmf, tmf0, tmf); // -nn*π
mpc_set_mpf(tmp0, zero, tmf); // -ivπ/2
mpc_exp(tmp0, tmp); // e^ -ivπ/2
mpc_mul(knza, tmp, knza);
mpf_div(n, two, tmf); // n / 2
mpf_int(tmf, tmf0); // int(n / 2)
mpf_sub(tmf, tmf0, tmf); // n / 2 - int(n / 2)
// nがゼロ以下で偶数だったら
if mpf_is_lt(nb, zero) and mpf_is_ne(tmf, zero) then // -1^n
mpc_chs(knza, ri) // 奇数時
else
mpc_copy(knza, ri); // 偶数時
mpc_clear5(a, b, c, inz, zc);
mpc_clear4(z2, knza, z24, mz24);
mpc_clear3(tmp, tmp0, tmp1);
mpf_clear(mone);
mpf_clear4(nnf, kf, tmf, tmf0);
mpf_clear3(tmf1, nb, n);
end;
// integerF true 整数 false 非整数
procedure second_or_therd_kind(vm: mp_float; xc : mp_complex; var ri : mp_complex; intgerF: Boolean);
var
tmp0, tmp1 : mp_complex;
begin
mpc_init2(tmp0, tmp1);
// 2nd kind 計算
if not intgerF then ekazMul(vm, xc, ri) // 非整数
else knzMul(vm, xc, ri); // 整数
// 3rd kind 計算
if form1.Check3rd.Checked = True then begin
mpc_set_mpf(tmp0, zero, pai); // iπ
mpc_mul_mpf(ri, two, tmp1); // ri * 2
mpc_div(tmp1, tmp0, ri); // ri * 2 / iπ
end;
mpc_clear2(tmp0, tmp1);
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;
//------------------------------------------------------------------------------
// kelvin-3rd function(x)のxの値は実数
// ta[]はグラフ計算点の配列
procedure TForm1.Button1Click(Sender: TObject);
label
EXT;
const
ta : array[0..GL - 1] of double = (0.08, 0.13, 0.23, 0.46, 0.86, 1.68, 2.7, 4);
YPM = 1E304; // doubleの最大値の制限 演算のオーバーフロー対策
YMM = -YPM;
var
ch, i, xi: integer;
herx, heix: double;
herxe, heixe: double;
x, v, xv: double;
xmin, xmax, dx, dxf: double;
ymin, ymax: double;
xm, vm, ndm, dxm: mp_float;
avm, andm: mp_float;
tmp : mp_float;
xc, epc, ri, tmc : mp_complex;
berl : Darray;
beil : Darray;
beru : Darray;
beiu : Darray;
xl : Darray;
xu : Darray;
vmaxmes : string;
GCF : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
mm, ss, ms : integer;
k : integer;
vk, tmf, tmf0 : mp_float;
GU, GS : integer;
intF : boolean; // 非整数 整数フラグ true 整数
// xの値をkelvin用の複素数に変換
procedure xtoxc(x, ia : mp_float; var xc : mp_complex);
begin
mpc_set_mpf(xc, x, zero); // xc
mpf_div(pai, four, tmp); // π/4
mpc_set_mpf(tmc, zero, tmp); // iπ/4
mpc_exp(tmc, epc); // e^(iπ/4)
mpc_mul(xc, epc, xc); // xe^(iπ/4)
mpc_div_mpf(xc, two, xc); // (xe^(iπ/4))/2
end;
// Doubleの最大値制限
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, v, ch);
if ch <> 0 then begin
application.MessageBox('次数vの値に間違いがあります。','注意', 0);
exit;
end;
if abs(v) > vmax then begin
vmaxmes := 'Vの値は ±' + intTostr(vmax) + '以下にして下さい。';
application.MessageBox(pwideChar(vmaxmes),'注意', 0);
exit;
end;
val(labelededit2.Text, x, ch);
if ch <> 0 then begin
application.MessageBox('xの値に間違いがあります。','注意', 0);
exit;
end;
if abs(x) > 100 then begin
vmaxmes := 'xの値は±100迄にしてください';
application.MessageBox(pchar(vmaxmes),'注意', 0);
exit;
end;
KMV := 50 + round(abs(x) * 10);
if KMV > KMmax then KMV := KMmax; // Σ 0～kMAX
mpf_init4(vm, xm, ndm, dxm);
mpf_init2(avm, andm);
mpf_init(tmp);
mpf_init3(vk, tmf, tmf0);
mpc_init4(xc, epc, ri, tmc);
mpf_read_decimal(vm, PAnsiChar(ansistring(labelededit1.Text + #00)));
mpf_read_decimal(xm, PAnsiChar(ansistring(labelededit2.Text + #00)));
rk := 0;
ik := 0;
memo1.Clear;
mpf_int(vm, ndm);
mpf_sub(vm, ndm, ndm);
mpf_abs(ndm, andm);
intF := True; // v 整数
if mpf_is_ne(ndm, zero) then intF := False; // v 非整数
// Γ値配列計算
mpf_abs(vm, avm);
if mpf_is_ne(ndm, zero) then begin
for k := 0 to KMV 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 KMV 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, MVG[k]); // Γ(n+k+1)
end;
end;
mpf_read_decimal(dxm, PAnsiChar(ansistring('1E-10' + #00)));
if mpf_is_lt(andm, dxm) and mpf_is_ne(ndm, zero) then begin
application.MessageBox('次数vの値の小数部は 1E-10 より大きくしてください。','注意', 0);
goto EXT;
end;
if mpf_is0(xm) then begin // xの値が0だったら
// ゼロ+側近辺
mpf_read_decimal(dxm, PAnsiChar(ansistring('1E-50' + #00)));
xToxc(dxm, zero, xc);
second_or_therd_kind(vm, xc, ri, intF);
if mpf_is_ge(ri.re, zero) then rk := infinity
else rk := -infinity;
if mpf_is_ge(ri.im, zero) then ik := infinity
else ik := -infinity;
end;
memo1.Lines.Append('v = ' + floatTostr(v) + ' x = ' + floatTostr(x));
StopWatch := TStopwatch.StartNew;
// 表示値の計算
herxe := 0;
heixe := 0;
// x <> 0
xToxc(xm, zero, xc);
if not mpc_is0(xc) then begin // xc <> 0
second_or_therd_kind(vm, xc, ri, intF);
memo1.Lines.Append('kerv(x) = ' + string(mpf_decimal(ri.re, 50)));
memo1.Lines.Append('keiv(x) = ' + string(mpf_decimal(ri.im, 50)));
herxe := mpf_todouble(ri.re);
heixe := mpf_todouble(ri.im);
end
// x = 0
else begin
// x = 0 abs(v) = 2
mpf_abs(vm, avm);
if mpf_is_eq(avm, two) then begin
mpf_inv(pai, tmp);
mpf_chs(tmp, tmp);
memo1.Lines.Append('herv(x) = ' + 'INF');
memo1.Lines.Append('heiv(x) = ' + string(mpf_decimal_alt(tmp, 50)));
herxe := infinity;
heixe := mpf_todouble(tmp);
end
else begin
// x = 0 v= 0
mpf_inv(two, tmp);
mpf_chs(tmp,tmp);
if mpf_is0(vm) then begin
memo1.Lines.Append('herv(x) = ' + string(mpf_decimal_alt(tmp, 50)));
memo1.Lines.Append('heiv(x) = ' + '-INF');
herxe := mpf_todouble(tmp);
heixe := -infinity;
end
else
// x = 0 v<> 0
if mpf_is_ne(vm, zero) then begin
memo1.Lines.Append('herv(x) = ' + floatTostr(rk));
memo1.Lines.Append('heiv(x) = ' + floatTostr(ik));
herxe := rk;
heixe := ik;
end;
end;
end;
StopWatch.Stop;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
// goto EXT;
ms := ElapsedMillseconds * (GL + 1) * 2 + 1000;
mm := ms div 60000;
ss := (ms mod 60000) div 1000;
vmaxmes := 'グラフ表示に約' + inttostr(mm) + '分' + inttostr(ss) + '秒必要です。';
memo1.Lines.Append(vmaxmes);
series1.Clear;
series2.Clear;
series3.Clear;
series4.Clear;
series5.Clear;
series6.Clear;
application.ProcessMessages;
if form1.Check3rd.Checked = True then
Chart1.Title.Caption := 'Kelvin functions of 3rd kind'
else
Chart1.Title.Caption := 'Kelvin functions of 2nd kind';
if (checkbox2.Checked = true) then begin // グラフ無だったら終了
Chart1.Canvas.Font.Style := [fsBold];
Chart1.Canvas.Font.size := 8;
Chart1.Canvas.TextOut(170, 115,'グラフ無し');
goto EXT;
end;
if ElapsedMillseconds > 500 then begin
vmaxmes := vmaxmes + #13#10 + 'グラフの表示をしますか?';
if MessageDlg(vmaxmes , mtConfirmation, [mbYes, mbNo], 0, mbYes) = mrNo then begin
Chart1.Canvas.Font.Style := [fsBold];
Chart1.Canvas.Font.size := 8;
Chart1.Canvas.TextOut(170, 115,'グラフ無し');
goto EXT;
end;
end;
// 最大値最小値の検索とグラフデーター作成
xi := round(x);
xmin := xi - 4;
GCF := 0;
if (x >= -2) and (x <= 2) then GCF := 1;
if (x >= -3) and (x <= -2) then GCF := 2;
if (x >= 2) and (x <= 3) then GCF := 3;
if (x >= -4) and (x <= -3) then GCF := 4;
if (x >= 3) and (x <= 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 := dxf;
GS := 40;
GU := 40;
ch := GCF;
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;
1 :
begin
for i := 0 to GL - 1 do begin
xu[i] := ta[i];
xl[GL - i - 1] := -ta[i];
end;
end;
else // 2,3,4,5
begin
for i := 0 to GL - 1 do begin
xu[i] := ta[i];
xl[GL - i - 1] := -ta[i];
end;
case ch 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;
for i := 0 to Gl - 1 do begin
xl[i] := xl[i] * dx;
xu[i] := xu[i] * dxf;
end;
end;
end;
// グラフ用データー作成
ymin := 0;
ymax := 0;
case GCF of
0 : begin
for i := 0 to GL + GL - 1 do begin
xv := xl[i];
mpf_set_dbl(tmp, xv);
xToxc(tmp, zero, xc);
second_or_therd_kind(vm, xc, ri, intF);
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];
mpf_set_dbl(tmp, xv);
xToxc(tmp, zero, xc);
second_or_therd_kind(vm, xc, ri, intF);
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];
mpf_set_dbl(tmp, xv);
xToxc(tmp, zero, xc);
second_or_therd_kind(vm, xc, ri, intF);
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 herxe > ymax then herxe := ymax;
if herxe < ymin then herxe := ymin;
if heixe > ymax then heixe := ymax;
if heixe < ymin then heixe := ymin;
series3.AddXY(x, herxe);
series4.AddXY(x, heixe);
if checkbox1.Checked = true then
case GCF of
0:
begin
for i := 0 to GL + GL - 1 do begin
series1.AddXY(xl[i], berl[i]);
series2.AddXY(xl[i], beil[i]);
end;
end;
1, 2, 3, 4, 5 :
begin
for i := 0 to GL - 1 do begin
series1.AddXY(xl[i], berl[i]);
series2.AddXY(xl[i], beil[i]);
series5.AddXY(xu[i], beru[i]);
series6.AddXY(xu[i], beiu[i]);
end;
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];
herx := maxmin(akima_Interpolation(xv));
series1.AddXY(xv, herx);
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];
heix := maxmin(akima_Interpolation(xv));
series2.AddXY(xv, heix);
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;
herx := maxmin(akima_Interpolation(xv));
series1.AddXY(xv, herx);
end;
for i := 0 to GL- 1 do yt[i] := beil[i];
akima_table;
for i := GS downto 1 do begin
xv := dx * i;
heix := maxmin(akima_Interpolation(xv));
series2.AddXY(xv, heix);
end;
// memo1.Lines.Append('(x) = ' + string(mpf_decimal_alt(xm, 50)));
if mpf_is0(xm) then begin
// series1.AddXY(0, herxe);
// series2.AddXY(0, heixe);
series5.AddXY(0, herxe);
series6.AddXY(0, heixe);
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;
herx := maxmin(akima_Interpolation(xv));
series5.AddXY(xv, herx);
end;
for i := 0 to GL- 1 do yt[i] := beiu[i];
akima_table;
for i := 1 to GU do begin
xv := dx * i;
heix := maxmin(akima_Interpolation(xv));
series6.AddXY(xv, heix);
end;
end;
end;
EXT:
mpf_clear4(vm, xm, ndm, dxm);
mpf_clear2(avm, andm);
mpf_clear(tmp);
mpf_clear3(vk, tmf, tmf0);
mpc_clear4(xc, epc, ri, tmc);
end;
procedure TForm1.FormCreate(Sender: TObject);
var
i : integer;
N, D, tmp : mp_float;
begin
mpf_set_default_decprec(200); // 有効桁数200桁 v=0 x= 100 時 50桁の精度に必要です。
setlength(BM, NB + 1); // ベルヌーイ数 配列
setlength(DG, KMmax + Vmax + 2); // ディガンマ
setlength(FA, KMmax + Vmax + 1); // 階乗
setlength(t, GL + GL); // akima 補間値計算 配列 t
setlength(m, GL + GL + 3); // akima 補間値計算 配列 m
setlength(PVG, KMmax + Vmax + 2); // +vΓ
setlength(MVG, KMmax + Vmax + 2); // -vΓ
for i := 0 to NB do mpf_init(BM[i]);
for i := 0 to KMmax + Vmax + 1 do mpf_init(DG[i]);
for i := 0 to KMmax + Vmax do mpf_init(FA[i]);
for i := 0 to KMmax + Vmax + 1 do mpf_init(PVG[i]);
for i := 0 to KMmax + Vmax + 1 do mpf_init(MVG[i]);
mpf_init3(N, D, tmp);
mpf_init(gamma);
mpf_init5(pai, zero, one, two, four);
mpf_init(log_2pis2);
mpf_set_pi(pai); // π
mpf_set0(zero); // 0
mpf_set1(one); // 1
mpf_set_int(two, 2); // 2
mpf_set_int(four, 4); // 4
mpf_mul(pai, two, tmp); // 2π
mpf_ln(tmp, tmp); // ln(2π)
mpf_div(tmp, two, log_2pis2); // ln(2π)/2
mpf_read_decimal(gamma, PAnsiChar(ansistring(gstr + #00))); // γ
for i := 0 to GL + GL - 1 do mpf_init(t[i]);
for i := 0 to GL + GL + 2 do mpf_init(m[i]);
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 + Vmax + 1 do begin
digammaMUL(i, D);
mpf_copy(D, DG[i]);
end;
for i := 0 to KMmax + Vmax do begin
factorialMul(i, N);
mpf_copy(N, FA[i]);
end;
memo1.Clear;
mpf_clear3(N, D, tmp);
end;
procedure TForm1.FormDestroy(Sender: TObject);
var
i : integer;
begin
for i := 0 to NB do mpf_clear(BM[i]);
for i := 0 to KMmax + Vmax + 1 do mpf_clear(DG[i]);
for i := 0 to KMmax + Vmax do mpf_clear(FA[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]);
for i := 0 to KMmax + Vmax + 1 do mpf_clear(PVG[i]);
for i := 0 to KMmax + Vmax + 1 do mpf_clear(MVG[i]);
mpf_clear(gamma);
mpf_clear5(pai, zero, one, two, four);
mpf_clear(log_2pis2);
end;
end.
kelvin_function_Mult_of_2nd_3rd_kind.zip
三角関数、逆三角関数、その他関数 に戻る