2024-03-30 Σ値の計算中のk乗の計算方法を変更し計算を速くしました。
2024-03-19
FA,MVG,PVGの配列の長さが長すぎたので修正しました。
ｘの値が複素数の場合の第１種ケルビン関数
ケルビン関数の場合、ｘの値が複素数の場合、一般的な複素数の取扱いと違い、通常の計算の実数に相当する複素数の値と、通常の計算の虚数部に相当する複素数の値となります。
左の式は、三種類ありますが、全て同じ計算です。
次数が非整数の場合、複素数平面で-45°＜φ≦45°の範囲の計算となります、計算結果はでますが、正しい値と成らない様です。
整数の値の場合は全ての値で正しい値を返します。
三番目の計算式場合、Σの計算の分子にsin、cosがあるので、収束判定に注意が必要です。
Kの値が大きくなると、Σ値が変化しなくなるのを検出して収束と判定するのですが、sin,
cosの為分子が0になる事があり、kの値が小さくても、Σ値が変化しない場合があります。
そこで、Σ値が連続して何回か同じ値が続いた場合を収束として判定とします。
左図は、プログラムの実行画面ですが、ｘの実数と、虚数部の入力が出来るようになっています。
入力を分けてあるのは、ｘ＋ｙｉ、ｘ－ｙｉの様に入力するようにしても良いのですが、入力処理が面倒なので、分けて入力するようにしました。
虚数部の計算結果は、beiv(x)のチェックボックスにチェックを入れると、計算されます。
プログラム
プログラムには、二種類の多倍長計算が組み込んでありますが、組み込み方は第１種ケルビン関数ｖ，ｘ実数を参照して下さい。
次のプログラム例は、三種類の計算式の一番上の計算例です。
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;
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);
mpf_set_int(two, 2); // 2
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_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;
// xの整数部とyの虚数部が0なら
if mpf_is0(x.re) and mpf_is0(y.im) then begin
mpf_frac(y.re, ai); // 乗数の小数部
if mpf_is0(ai) then begin // 整数なら
mpf_div(y.re, two, ay); // 偶数奇数確認
mpf_frac(ay, ai); // 小数部
if mpf_is0(ai) then mpf_set0(ans.im) // 偶数なら 虚数部0
else mpf_set0(ans.re); // 奇数なら実数部0
end;
end;
mpf_abs(x.re, ay);
mpf_abs(x.im, ai);
// xの実数部と虚数部の絶対値が等しくてyがゼロでなかったら
if mpf_is_eq(ay, ai) and mpf_is0(y.im) and not mpc_is0(y) then begin
mpf_frac(y.re, ai); // 乗数の小数部
if mpf_is0(ai) then begin // 整数なら
mpf_div(y.re, two, ay); // 偶数奇数確認
mpf_frac(ay, ai); // 小数部
if mpf_is0(ai) then begin
mpf_add(two, two, two); // 4
mpf_div(y.re, two, ay); // 4偶数奇数確認
mpf_frac(ay, ai); // 小数部
if mpf_is0(ai) then mpf_set0(ans.im) // 偶数なら 虚数部0
else mpf_set0(ans.re); // 偶数なら 虚数部0
end;
end;
end;
mpf_clear5(ay, ai, zero, harf, two);
end;
//------------------------------------------------------------------------------
// jv(x) 多倍長
// X 値 複素数
// v 次数
// F False 第一種ベッセル関数 true 変形
procedure Jvx(var v : mp_float; var x, ri: mp_complex; F: boolean);
var
k : integer;
s, x24k, tmp, tmp0, tmp1 : mp_complex;
kc, 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_init3(kc, 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) // Γ(v+k+1)
else
mpf_copy(MVG[k], tmf); // Γ(-v+k+1)
mpf_mul(FA[k], tmf, khg); // k!Γ(n+k+1) vkが0,負の整数の時±∞
// mpc_set_mpf(kc, kf, zero); // *** kc = k + 0i 複素数
// mpc_powa(xc, kc, x24k); // *** (i(x^2)/4)^k
mpc_div_mpf(x24k, khg, tmp0); // ((i(x^2)/4)^k)/(k!Γ(v+k+1))
if not F 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_clear3(kc, xc, xs2);
mpf_clear2(khg, kf);
mpf_clear3(kd, nd, vk);
mpf_clear2(tmf, tmf0);
mpc_clear(sb);
end;
// z 値 複素数
// v 次数
// F False berv(z) True beiv(z)
procedure berivz(v : mp_float; var x, ri: mp_complex; F : boolean);
var
xehpipis4, xehmipis4, ehpipiv, ehmipiv : mp_complex;
jvs, jve : mp_complex;
tmpc0, tmpc1 : mp_complex;
tmpf0, tmpf1 : mp_float;
begin
mpc_init4(xehpipis4, xehmipis4, ehpipiv, ehmipiv);
mpc_init2(jvs, jve);
mpc_init2(tmpc0, tmpc1);
mpf_init2(tmpf0, tmpf1);
mpf_div(pai, four, tmpf0); // π/4
mpc_set_mpf(tmpc0, zero, tmpf0); // 0+iπ/4
mpc_exp(tmpc0, xehpipis4); // e^(iπ/4)
mpc_mul(xehpipis4, x, xehpipis4); // e^(iπ/4)x *
mpc_chs(tmpc0, tmpc1); // 0-iπ/4
mpc_exp(tmpc1, xehmipis4); // e^(-iπ/4)
mpc_mul(xehmipis4, x, xehmipis4); // e^(-iπ/4)x *
mpf_mul(pai, v, tmpf1); // πv
mpc_set_mpf(tmpc0, zero, tmpf1); // iπv
mpc_exp(tmpc0, ehpipiv); // e^(iπv) *
mpc_chs(tmpc0, tmpc1); // -iπv
mpc_exp(tmpc1, ehmipiv); // e^(-iπv) *
if not F then begin
Jvx(v, xehmipis4, jvs, false); // jxv(v,e^(-iπ/4))
mpc_mul(jvs, ehpipiv, jvs); // e^(iπv) * jvs
Jvx(v, xehpipis4, jve, false); // jxv(v,e^(iπ/4))
mpc_mul(jve, ehmipiv, jve); // e^(-iπv) * jve
mpc_add(jvs, jve, tmpc0); // jvs+jve
mpc_div_mpf(tmpc0, two, ri); // (jvs+jve) /2
end
else begin
Jvx(v, xehpipis4, jvs, false); // jvx(v,e^(iπ/4))
mpc_mul(jvs, ehmipiv, jvs); // e^(-iπv) * jvs
Jvx(v, xehmipis4, jve, false); // jve(v,e^(-iπ/4))
mpc_mul(jve, ehpipiv, jve); // e^(iπv) * jve
mpc_sub(jvs, jve, tmpc0); // jvs-jve
mpc_div_mpf(tmpc0, two, tmpc1); // (jvs+jve) /2
mpc_set_mpf(tmpc0, zero, one); // i
mpc_mul(tmpc1, tmpc0, ri) // i(jvs+jve) /2
end;
mpc_clear4(xehpipis4, xehmipis4, ehpipiv, ehmipiv);
mpc_clear2(jvs, jve);
mpc_clear2(tmpc0, tmpc1);
mpf_clear2(tmpf0, tmpf1);
end;
// Kelvin function
// F False berv(z) True beiv(z)
// v 次数
// z 値 複素数
procedure kelvincomplexNo1(var v : mp_float; var z, ri: mp_complex; F: boolean);
var
ar, ai, vn : mp_float;
begin
mpf_init3(ar, ai, vn);
berivz(v, z, ri, F);
mpf_abs(z.re, ar);
mpf_abs(z.im, ai);
mpf_frac(v, vn);
// zの実数部と虚数部の絶対値が等しく次数が整数だったら
if mpf_is_eq(ar, ai) and mpf_is0(vn) then begin
mpf_div(v, two, ar);
mpf_frac(ar, ai);
if mpf_is0(ai) then // 偶数だったら
if not F then mpf_set0(ri.im)
else mpf_set0(ri.re);
end;
if mpf_is0(z.im) then mpf_set0(ri.im);
// 実数部がゼロだったら
if mpf_is0(z.re) and not mpf_is0(z.im) then begin
mpf_frac(v, ar);
if mpf_is0(ar) then begin // 次数整数だったら
mpf_div(v, two, ai);
mpf_frac(ai, ar);
if mpf_is0(ar) then mpf_set0(ri.im) // 次数偶数だったら
else mpf_set0(ri.re); // 次数奇数だったら
end;
end;
mpf_clear3(ar, ai, vn);
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;
F : boolean;
// 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);
F := False;
if Checkbox3.Checked then F := true;
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)));
series1.Clear;
series2.Clear;
series3.Clear;
series4.Clear;
series5.Clear;
series6.Clear;
mpf_int(vm, tmf);
mpf_sub(vm, tmf, tmf0);
if not mpf_is0(tmf0) then begin
if mpf_is_lt(xm, zero) then begin
istr := 'vの値が整数でない場合、xの実数部は正数でないと' + #13#10 +
'正しい答えが得られません計算を中断しますか';
if MessageDlg(istr , mtConfirmation, [mbYes, mbNo], 0, mbYes) = mrYes then goto EXT;
end;
if mpf_is_gt(xm, zero) then begin
mpf_div(ixm, xm, tmf);
mpf_chs(one, tmf0);
if mpf_is_gt(tmf, one) or mpf_is_le(tmf, tmf0) then begin
istr := 'vの値が整数でない場合、xの実数部が正の値で' + #13#10 +
'複素数平面で-45°<ψ≦45°の条件でないと' + #13#10 +
'正しい答えが得られません計算を中断しますか';
if MessageDlg(istr , mtConfirmation, [mbYes, mbNo], 0, mbYes) = mrYes then goto EXT;
end;
end;
end;
mpf_read_decimal(zz, PAnsiChar(ansistring(zeromg + #00)));
mpc_set_mpf(xcb, xm, ixm); // xcb 計算用 xの複素数
if ix >= 0 then istr := ' +' + floatTostr(ix) + 'i'
else istr := ' ' + floatTostr(ix) + 'i';
memo1.Lines.Append('n = ' + 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
kelvincomplexNo1(vm, xc, ri, F); // ゼロ近辺のプラス側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;
// 表示値の計算
if not mpc_is0(xcb) then
kelvincomplexNo1(vm, xcb, ri, F); // xcb 複素数
// x=0 v < 0 v = 非整数 の場合 ±∞
if mpc_is0(xcb) then begin
if mpf_is_lt(vm, zero) and mpf_is_ne(nd, zero) then begin
berxe := rk;
beixe := ik;
end
else begin
berxe := 0;
beixe := 0;
end;
if mpf_is0(vm) and not F then berxe := 1;
if checkbox3.Checked = false then
memo1.Lines.Append('berv(x) = ' + floatTostr(berxe)) // 数値表示
else
memo1.Lines.Append('beiv(x) = ' + floatTostr(berxe)); // 数値表示
end
else begin
if checkbox3.Checked = false then
memo1.Lines.Append( string('berv(x) = ' + mpf_decimal(ri.re, 50)))
else
memo1.Lines.Append( string('beiv(x) = ' + mpf_decimal(ri.re, 50)));
if not mpf_is0(ixm) then begin
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'));
beixe := mpf_todouble(ri.im);
end
else
beixe := 0;
berxe := mpf_todouble(ri.re);
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;
if mpf_is_ge(xm, ixm) and mpf_is_ge(ixm, zero) then begin
mpf_int(vm, tmf);
mpf_sub(vm, tmf, tmf0);
if not mpf_is0(tmf0) then begin
xmin := mpf_todouble(ixm);
GCF := 0;
end;
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);
kelvincomplexNo1(vm, xc, ri, F);
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);
kelvincomplexNo1(vm, xc, ri, F);
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);
kelvincomplexNo1(vm, xc, ri, F);
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.
kelvin_function_ComplexS_of_1st_kind.zip
三角関数、逆三角関数、その他関数 に戻る