2023-11-04
プログラム整理と、演算精度の向上
2024-10-08
zの値が、-1から-2の時、高速に計算が出来ないため、1/(1-z)の線形接続公式を追加しました。
ガウスの超幾何関数 (Hypergeometric
function)
前の超幾何関数の修正版です。
(x)n はポッホハマー記号で表した昇冪 (x)0 = 1、(x)n = x (x+1) (x+2)…(x+n−1)
値は、プログラムを組んで、n=0、n=1、n=2、と順に計算Σ値を求めて、Σ値が変化しなくなるまで計算することで、簡単に求める事が出来ます。
オイラー積分表示
オイラー積分表示は、上記ですが、ガウスの超幾何関数の資料に必ずでてくるのでプログラムに組み込んで確認してみました。
積分出来る条件は、a>=1 c-a>=1 |z|<1になり計算で来る範囲が可成り狭くなります、又、a,bの虚数部特にbの虚数部の値が大きくなると計算誤差が大きくなるようです。
上記の線形接続公式では、Γ(x)のｘの値が0(ZERO)或いは0以下の整数であると、Γ値が±∞になり正しい値を返さないことになります。
線形接続公式には、殆どの場合ガンマ関数の分数式が入ります。
前回の超幾何関数のデバッグ時に、ガンマ関数±∞の時の計算に法則があることに気が付きました。
例えば、
Γ(0)/Γ(-1) = -1 = 1/-(1!)
Γ(0)/Γ(-2) = 0.5 = 1/(2!)
Γ(0)/Γ(-3) = -0.166667 = 1/-(3!)
Γ(-2)/Γ(-3) = -3
= {1/(2!) ]/[1/-(3!)]
F(-3) = 1/-3 * 1/-2
* 1/-1 * 1/0
F(-2) = 1/-2 * 1/-1 * 1/0
F(-2)/Γ(-3) = (1/-2 * 1/-1 * 1/0) / (1/-3 * 1/-2 * 1/-1 * 1/0)
=
1 / (1/-3)
= -3
F(-2)×Γ(-3) = 1/-3 *
(1/-2 * 1/-1 * 1/0)2 = -∞
の様にガンマ関数に与えられた0以下の整数の階乗の逆数に従う事になります。
奇数の場合と、偶数の場合で符号の扱いに注意が必要です。
分子と分母側に同じ数の0以下の整数のガンマ関数がある場合は、階乗の計算で正しい値を求めることが出来ます。
分母の方が多い場合は、1/±∞となり限りなくゼロに近づき、分子の方が多い場合は、±∞に近づきます。
限りなく無限大に近づくと言うことから、近づいて行く途中の、適切な値を値を与えることで、正解に近い値(近似値)を求めることが出来ます。
なるべく正解に近づく様に、+Δ値と-Δ値の平均値をとって、値とします。
Γ分数計算値が無限大になる場合、僅かでも、0或いは負の整数から外すと、答えが得られる事に疑問があり、検討した結果、a,b,,c,zの値は同じても、Γ分数計算値が無限大になるものとならない、線形接続公式があり、双方で、計算、+Δ値と-Δ値の平均値で近似値が計算出来る事を確認しました。
又、±∞の場合 [2/0] /
[4/0] = 0.5 の様な計算がΓ関数では成り立つ事になります。
基本的に0での除算は出来ないのですが、もし、ゼロでの除算を∞として与えるならば、分子の値を伴って∞の値を与える必要があることになります。
Γ関数の計算の上記法則を取り入れて、
新しくプログラムを組んでみました。
近似値の計算をしなくて良い場合が、格段に多くなります。
上記法則を取り入れたのは、1/Zの線形接続公式のみです。
上記以外で使用している公式は次のようなものです。
線形接続公式のどれを使用するかは、2F1(a,b,c,z) z<1の計算値を渡すとき、いかにz相当の値を小さくするかにあります。
左図の演算桁数は、演算上の有効桁数で、答えの有効桁数ではありません。
一応殆どの値に対して、答えを返しますが、近似計算時の精度は不明です。
zの値が1以下であれば、時間は掛かるが、誤差の無い計算ができるのですが、1以上の、近似計算になると、一部を除いて正解の確認が出来ません。
計算速度を上げる為、pfaff変換を使用している場合があります、その時はpfaff変換を外すと、計算は遅いですが正確な値を返す場合があります。
近似計算の場合、a=bの場合も同様で、一度チェックを外して計算をしてみて下さい。
a,a+1,c or b,b+1,c は1-1/zの計算を行うのですが、誤差が大きいです、なるべく使用を避けてください。
プログラム上チェックを外していても a,a+1,c or b,b+1,c の計算をする場合があります。
* a, b=a+2, c=a+1, z=c
の時 計算結果がゼロになります、近似計算で無くても、ゼロになるので、疑問があるのですが、原因はわかっていません。
z値がを1越えていく時、1より小さくて、1に近いときは、+の大きな値となり1を越えて1に近いとき大きな－となり 1 と c の間では-となり、c=z でゼロとなり c を越えると+になり、更に大きくするとゼロに近づくことになります。
プログラムで近似値の計算になるとき、精度に応じて、表示する桁数の変更をするようにしてあります。
実際には高精度で計算結果が出ていても、表示する桁数が少ない場合があります。
計算選択のチェックボックスがありますが、選択を変更することで、近似計算ではなく正確な値が計算できる場合があります。
問題点
z=0.5+√0.75i |z|=1 の時は、特殊な場合を除き収束しません、この値を高速に計算する方法が分かりません。
z=0.6+0.8i |z|=1ですが、これは容易に収束します。
プログラム
デバッグは、不完全です注意してください。(余りにも大変なので打ち切りました。)
プログラムには、多倍長演算が組み込んであります。
組み込み方は第１種ケルビン関数を参照して下さい。
// Gaussの超幾何関数 // 値が1より小さい場合は必ず収束しますが // 1の場合は専用の計算があります、a又はbの値が負数の整数の場合は収束します。 // 1より大きい場合は、線形接続公式を計算しますが、計算は難しくなります。 // Zの値により、収束の速い線形接続公式を利用すれば良いのですが、ここでは。 // 三種類にしています。 // 二次変換も一部使用しています。 // arctan(x)はxの値が1より大きくても必ず正しい値をかえします。 // Z > 1の場合Γ(x)関数を使用するのですが、xの値が0又は負の整数となると // Γ(x) +∞ or -∞となるので微小値を加減算して整数から外して計算する場合があります // +Δと-Δの値を計算し平均値を近似値として答えとしていますが、誤差が大きくなります。 // 特にa=b c-a=1(c-b=1)の時は難しくなります。 // ガンマ値計算の分数部分で分子と分母の∞の数が等しい場合は近似計算は必要ありません。 // ガンマ計算時分子のみにガンマ値±∞がが入る場合は推定値の計算となる為、精度が低くなります。 // a, b-a+2, c=a+1, z=c z >= 2 の時 計算結果がゼロになります正しいのか疑問があります。 // z = 0.5 /z/=1 (z.im = √0.75) 時は収束しません 、近似の時も収束が遅くなりますが // 速くする方法が分かりません。 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, Vcl.Buttons, system.UITypes, system.Math; type TForm1 = class(TForm) Memo1: TMemo; aEdit: TLabeledEdit; bEdit: TLabeledEdit; cEdit: TLabeledEdit; zEdit: TLabeledEdit; Button1: TButton; iaEdit: TLabeledEdit; ibEdit: TLabeledEdit; icEdit: TLabeledEdit; izEdit: TLabeledEdit; PrecisionEdit: TLabeledEdit; Label1: TLabel; epsEdit: TLabeledEdit; BitBtn1: TBitBtn; artanBtn: TBitBtn; ArtanedEdit: TLabeledEdit; arsinEdit: TLabeledEdit; arsinBtn: TBitBtn; CheckBox1: TCheckBox; CheckBox2: TCheckBox; CheckBox3: TCheckBox; Button2: TButton; CheckBox4: TCheckBox; CheckBox5: TCheckBox; CheckBox6: TCheckBox; Label2: TLabel; CheckBox7: TCheckBox; CheckBox8: TCheckBox; procedure Button1Click(Sender: TObject); procedure FormCreate(Sender: TObject); procedure BitBtn1Click(Sender: TObject); procedure FormClose(Sender: TObject; var Action: TCloseAction); procedure artanBtnClick(Sender: TObject); procedure arsinBtnClick(Sender: TObject); procedure Button2Click(Sender: TObject); procedure FormCloseQuery(Sender: TObject; var CanClose: Boolean); private { Private 宣言 } function inputcheck(var pre: integer): boolean; public { Public 宣言 } end; var Form1: TForm1; implementation {$R *.dfm} uses mp_real, mp_cmplx, mp_types, mp_base, Velthuis.BigIntegers; const NB = 120; // ベルヌーイ数 配列数 NB + 1 var BM : array of mp_float; // ベルヌーイ数配列 NumeratorString : array[0..NB] of string; // 分子 DenominatorString : array[0..NB] of string; // Γ関数用分母 zero, one, two, three, four : mp_float; pai, maxmpf, infs : mp_float; log_2pis2 : mp_complex; // 最大公約数 ユークリッドの互除法 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]; // logΓ関数用に分母値Bの値計算 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_complex); var v, w, cone, ctwo : mp_complex; tmp, tmp0, s : mp_complex; i : integer; begin mpc_init4(v, w, cone, ctwo); mpc_init3(tmp, tmp0, s); mpc_set1(cone); // 1 + 0i mpc_add(cone, cone, ctwo); // 2 + 0i mpc_set1(v); // 1 + 0i mpc_set0(tmp); // 0 + 0i mpf_set_int(tmp.re, NB); // NB + 0i while mpf_is_lt(x.re, tmp.re) do begin // while x < NB do mpc_mul(v, x, v); // v = v * x mpf_add(x.re, one, x.re); // x := x + 1 end; mpc_mul(x, x, tmp); // x^2 mpc_div(cone, tmp, w); // w = 1 / x^2 mpc_set0(s); // s = 0 + 0i for i := NB downto 1 do begin mpc_add_mpf(s, BM[i], tmp); // tmp = s + B[i] mpc_mul(tmp, w, s); // s = tmp * w end; mpc_add_mpf(s, BM[0], tmp); // tmp = s + B[0] mpc_div(tmp, x, s); // s = (s + B[0]) / x mpc_add(s, log_2pis2, s); // s = s + ln(2π)/2 mpc_ln(v, tmp); // ln(v) mpc_sub(s, tmp, s); // s := s - ln(v) mpc_sub(s, x, s); // s := s - x mpc_div(cone, ctwo, tmp); // tmp = 1/2 mpc_sub(x, tmp, tmp0); // tmp0 = x - 1/2 mpc_ln(x, tmp); // ln(x) mpc_mul(tmp0, tmp, tmp0); // tmp0 = (x - 1/2) * ln(x) mpc_add(s, tmp0, ans); // ans = s + (x - 1/2) * ln(x) mpc_clear4(v, w, cone, ctwo); mpc_clear3(tmp, tmp0, s); end; // 多倍長ガンマ 複素数 // xの値が 0 と負整数の時Γは∞になるので∞フラグセット値は1を返します。 function gammaMul(xx : mp_complex; var ans: mp_complex) : integer; label EXT; var tmp, sinx, logG, cone, cpai : mp_complex; x : mp_complex; begin mpc_init5(tmp, sinx, logG, cone, cpai); mpc_init(x); // xxの値を変化させるxxの元値を解放不可になるのでコピーして使用します mpf_read_decimal(x.re, PAnsiChar(ansistring(string(mpf_decimal(xx.re, 200)) + #00))); mpf_read_decimal(x.im, PAnsiChar(ansistring(string(mpf_decimal(xx.im, 200)) + #00))); // 虚数部がゼロで実数部がゼロを含み負の整数だったら∞ if mpf_is0(x.im) and mpf_is_le(x.re, zero) then begin mpf_frac(x.re, tmp.re); // 実数部の小数点以下取り出し if mpf_is0(tmp.re) then begin // 小数点以下がゼロだったら mpf_copy(maxmpf, ans.re); mpf_div(x.re, two, tmp.re); mpf_frac(tmp.re, tmp.re); if not mpf_is0(tmp.re) then mpf_chs(ans.re, ans.re); // mpf_set_dbl(ans.re, maxdouble); mpf_set0(ans.im); // mpc_set0(ans); result := 1; // ∞フラグセット goto EXT; // 終了 end; end; mpc_set_mpf(cpai, pai, zero); // π+ 0i mpc_set1(cone); // 1 + 0i if mpf_is_lt(x.re, zero) then begin // x.real < 0 mpc_mul_mpf(x, pai, tmp); // x*π mpc_sin(tmp, sinx); // sin(πx); mpc_sub(cone, x, tmp); // 1-x log_GammaMul(tmp, logG); // loggamma(1-x); mpc_exp(logG, logG); // exp(loggamma) mpc_div(cpai, sinx, tmp); // π/sin(πx) mpc_div(tmp, logG, ans); // ans = π/(sin(πx) * logG(1-x)) end else begin // x.real >= 0 log_GammaMul(x, logG); // logGamma(x) mpc_exp(logG, ans); // exp(logGamma) end; result := 0; // 成功フラグ値 EXT: mpc_clear5(tmp, sinx, logG, cone, cpai); mpc_clear(x); end; // x(-0.5, 0) y(1.5, 0) 有効桁数70前後 // x(-0.5, -0.5) y(4, 0) 有効桁数60前後 // x(-0.5, 0.5) y(2, 0) 有効桁数70前後 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 y= 4 8 else mpf_set0(ans.re); // 偶数なら 虚数部0 y= 2 6 10 end; end; end; mpf_clear5(ay, ai, zero, harf, two); end; const NC = 10000; var BR : boolean; // 計算打切りフラグ TN : integer; // 計算実行 0 通常 1 artan 2 arsin DS : integer; // 計算表示 0 通常 1 artan 2 arsin procedure comment; var str : string; begin with form1.Memo1 do begin Clear; str := ' ガウスの超幾何関数'; lines.Append(str); str := ' |z|の値が1に近づくと、収束に時間が掛かります。'; lines.Append(str); str := 'その時にPfaffの変換をすると、収束が速くなるのですが、a,b,cの値に整数であると'; str := str + '近似値計算となる場合があります。'; lines.Append(str); str := 'zの値が0.85程度であれば、Pfaffの変換を使用しなくても、結構速く収束するので'; str := str + 'Pfaffの変換のある場合とない場合で同じ答えがでるか確認して下さい。'; lines.Append(str); str := '答えが違う場合は、その時の"a,b,c"値ではPfaffの変換は使用できません。'; lines.Append(str); str := ' 無限大回避は、"|z|>1"の時Γ関数を使用するのですが、Γ(x)のxの値が"0"又は'; str := str + '負の整数の時±無限大になるのを避ける為のものです。'; lines.Append(str); str := '無限大回避を使用した時と'; str := str + 'しない場合で答えが違う場合は、無限大回避を使用してください。'; lines.Append(str); str := ' Pfaff変換をした場合、計算上"|z|"の値が1を越え、Γ(x)の計算が使用されます。'; lines.Append(str); str := '無限大回避を使用した場合としない場合で答えが違う場合は、無限大回避を使用してください。'; lines.Append(str); end; end; // オイラー積分 Z < 1 // RE(c)>RE(a)>0、 RE(a)>=1 RE(c-a)>=1 でないと正しく計算出来ません。 // RE(a-1)>=0 RE(c-a-1)>=0 // aとcの虚数部の値の差が大きくなる値誤差が大きくなります。 procedure Euler_integral(a, b, c, z : mp_complex; var ans : mp_complex); var am1, cma, cmam1, mb, onemt : mp_complex; gc, ga, gcma, tham1, onemthcmam1 : mp_complex; onemtz, onemtzhmb, s, tmc, dt : mp_complex; cone, t, st : mp_complex; n, nmax : integer; nfmax : mp_float; begin mpc_init5(am1, cma, cmam1, mb, onemt); mpc_init5(gc, ga, gcma, tham1, onemthcmam1); mpc_init5(onemtz, onemtzhmb, s, tmc, dt); mpc_init3(cone, t, st); mpf_init(nfmax); mpc_set1(cone); nmax := 2000; // nmax mpc_sub(a, cone, am1); // a-1 mpc_sub(c, a, cma); // c-a mpc_sub(cma, cone, cmam1); // c-a-1 mpc_chs(b, mb); // -b mpf_set_int(nfmax, nmax); // mpf nmax mpc_set0(s); mpc_div_mpf(cone, nfmax, dt); // dt for n := 0 to nmax do begin mpf_set_int(nfmax, n); // mpf nmax mpc_mul_mpf(dt, nfmax, t); // t if n = 0 then mpc_set0(tham1) // t = 0 時 0^(re=0,im<>0) エラー対策 else mpc_powa(t, am1, tham1); // t^(a-1) mpc_sub(cone, t, onemt); // 1-t mpc_mul(t, z, tmc); // tz mpc_sub(cone, tmc, onemtz); // 1-tz mpc_powa(onemtz, mb, onemtzhmb); // (1-tz)^(mb) if n = nmax then mpc_set0(onemthcmam1) // 1-t = 0 時 0^(re=0,im<>0) エラー対策 else mpc_powa(onemt, cmam1, onemthcmam1); // (1-t)^(c-a-1) mpc_mul(tham1, onemthcmam1, tmc); // t^(a-1) * 1-t)^(c-a-1) mpc_mul(tmc, onemtzhmb, st); // st = t^(a-1) * 1-t)^(c-a-1) * (1-tz)^(mb) mpc_add(s, st, s); // s = s + st end; gammaMul(c, gc); // Γ(c) gammaMul(a, ga); // Γ(a) gammaMul(cma, gcma); // Γ(c-a) mpc_div(gc, ga, tmc); // Γ(c)/Γ(a) mpc_div(tmc, gcma, tmc); // Γ(c)/Γ(a)/Γ(c-a) mpc_mul(s, tmc, ans); // Γ(x) * s mpc_mul(ans, dt, ans); // Γ(x) * s * dt // form1.memo1.Lines.Append('Euler integral = re' + string(' ' + mpf_decimal(ans.re, 4))); // form1.memo1.Lines.Append('Euler integral = im' + string(' ' + mpf_decimal(ans.im, 4)) + ' i'); mpc_clear5(am1, cma, cmam1, mb, onemt); mpc_clear5(gc, ga, gcma, tham1, onemthcmam1); mpc_clear5(onemtz, onemtzhmb, s, tmc, dt); mpc_clear3(cone, t, st); mpf_clear(nfmax); end; // ガウスの超幾何関数 z < 1 procedure Hypergeometric_function(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); label EXT; var ap, bp, cp, zhn : mp_complex; ah, bh, ch, ab : mp_complex; s, sb, tmp, tmpb : mp_complex; nbig, nk : mp_float; smla, smlb, ftmp : mp_float; dpcs : integer; begin mpc_init4(ap, bp, cp, zhn); mpc_init4(ah, bh, ch, ab); mpc_init4(s, sb, tmp, tmpb); mpf_init2(nbig, nk); mpf_init3(smla, smlb, ftmp); dpcs := mpf_get_default_prec; mpc_set0(s); // 0 mpf_set0(smla); // 0 mpf_set0(smlb); // 0 mpc_copy(s, sb); mpf_set1(nk); // n! = 1 mpc_copy(z, zhn); // n = 1時値セット n=0 は1なので後で加算 mpc_copy(a, ah); // a(a+1)(a+2)(a+n-1) mpc_copy(b, bh); mpc_copy(c, ch); mpc_copy(a, ap); // a+n-1 mpc_copy(b, bp); mpc_copy(c, cp); mpf_set1(nbig); // n = 1 n := 1; // loop数 =1 mpc_set0(tmp); // 計算値初期化 0 repeat // n = 1から計算 mpc_mul(ah, bh, ab); // (a)n * (b)n if mpc_is0(ab) and not mpc_is0(ch) then begin // (a)n * (b)n = 0 ch <> 0なら if f <> 8 then f := 4; // z=1 ∞符号確認でなかったら a*b がゼロフラグセット終了 break; end; if mpc_is0(ch) then begin // (c)n = 0 なら mpc_mul(ab, zhn, ans); if f <> 8 then // z=1 ∞符号確認でなかったら if mpf_is_ge(ans.re, zero) then f := 1 // ∞ 無限大フラグセットして終了 else f := -1; // -∞ 無限大フラグセットして終了 if mpc_is0(ab) then begin f := 1; mpc_set1(ab); end; break; end; mpc_copy(s, sb); mpc_copy(tmp, tmpb); mpc_set0(tmp); // Σ 1～ (a(n)b(n)Z^n)/(c(n)n^n) mpc_div_mpf(ab, nk, tmp); // (a)n * (b)n / n! mpc_mul(tmp, zhn, tmp); // (a)n * (b)n / n! * z^n mpc_div(tmp, ch, tmp); // (a)n * (b)n / n! * z^n / (c)n mpc_add(s, tmp, s); // Σ // form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(tmp.re, 20))); mpf_add(ap.re, one, ap.re); // a = a + 1 mpf_add(bp.re, one, bp.re); // b = b + 1 mpf_add(cp.re, one, cp.re); // c = c + 1 mpc_mul(ah, ap, ah); // a(a+1)(a+2)(a+3)(a+n-1) mpc_mul(bh, bp, bh); // b(b+1)(b+2)(b+3)(b+n-1) mpc_mul(ch, cp, ch); // c(a+1)(c+2)(c+3)(c+n-1) mpc_mul(zhn, z, zhn); // z^n inc(n); // ループカウンターインクリメント mpf_add(nbig, one, nbig); // n = n + 1 mpf_mul(nk, nbig, nk); // n! mpc_abs(tmp, ftmp); if n mod 100 = 0 then begin // 指定回数に途中経過表示 form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(ftmp, 20))); application.ProcessMessages; if BR then begin form1.memo1.Lines.Append('途中停止しました。'); break; end; end; // 収束判定 until (mpf_is_lt(ftmp, eps) or (n >= NC)); if f = 8 then mpc_copy(s, ans); // f= 8 z=1 ∞符号確認用 if (f = 1) or (f = -1) or (f = 8) then else mpc_add_mpf(s, one, ans); // n = 0の値1を加算 form1.memo1.Lines.Append('loop数 = ' + intTostr(n)); EXT: mpc_clear4(ap, bp, cp, zhn); mpc_clear4(ah, bh, ch, ab); mpc_clear4(s, sb, tmp, tmpb); mpf_clear2(nbig, nk); mpf_clear3(smla, smlb, ftmp); end; // z > 1 a = b, c <> a // a = b c <> a 専用ルーチン procedure Hypergeometric_function_of_first_kind_aeb(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); label EXT; var ga, gb, gc, gcmamb, gapbmc : mp_complex; ma, amcpone, apbmcpone, cma, onema : mp_complex; onemz, cmamb, chsa, amc, onemonesz : mp_complex; cone, onmzhcab, zhma, tmc : mp_complex; cmb, gcmb, cmambpone, gcma : mp_complex; apbmc, zhamc : mp_complex; fa, fb, g1, g2: mp_complex; afa, afb : mp_complex; f1, f2 : integer; begin mpc_init5(ga, gb, gc, gcmamb, gapbmc); mpc_init5(ma, amcpone, apbmcpone, cma, onema); mpc_init5(onemz, cmamb, chsa, amc, onemonesz); mpc_init5(cone, onmzhcab, zhma, tmc, cmambpone); mpc_init5(cmb, gcmb, gcma, apbmc, zhamc); mpc_init4(fa, fb, g1, g2); mpc_init2(afa, afb); mpc_set1(cone); mpc_sub(c, a, cma); // c-a mpc_sub(cma, b, cmamb); // c-a-b mpc_sub(a, c, amc); // a-c mpc_sub(c, b, cmb); // c-b mpc_add(a, b, tmc); mpc_sub(tmc, c, apbmc); // a+b-c mpc_chs(a, chsa); // -a mpc_sub(cone, z, onemz); // 1-z mpc_add(amc, cone, amcpone); // a-c+1 mpc_add(amc, b, tmc); mpc_add(tmc, cone, apbmcpone); // a+b-c+1 mpc_sub(cone, a, onema); // 1-a mpc_sub(cma, b, tmc); mpc_add(tmc, cone, cmambpone); // c-a-b+1 mpc_inv(z, tmc); mpc_sub(cone, tmc, onemonesz); // 1-1/z gammaMul(a, ga); // Γ(a) gammaMul(b, gb); // Γ(b) gammaMul(c, gc); // Γ(c) gammaMul(cmamb, gcmamb); // Γ(c-a-b) gammaMul(cma, gcma); // Γ(c-a) gammaMul(cmb, gcmb); // Γ(c-b) gammaMul(apbmc, gapbmc); // Γ(a+b-c) mpc_powa(z, chsa, zhma); // z^-a mpc_powa(onemz, cmamb, onmzhcab); // (1-z)^(c-a-b) mpc_powa(z, amc, zhamc); // z~(a-c) mpc_mul(gc, gcmamb, g1); // Γ(c)Γ(c-a-b) mpc_mul(gcma, gcmb, tmc); // Γ(c-a)Γ(c-b) mpc_div(g1, tmc, g1); // g1 =Γ(c)Γ(c-a-b)/(Γ(c-a)Γ(c-b)) mpc_mul(gc, gapbmc, g2); // Γ(c)Γ(a+b-c) mpc_mul(ga, gb, tmc); // Γ(a)Γ(b) mpc_div(g2, tmc, g2); // g2 =Γ(c)Γ(a+b-c)/(Γ(a)Γ(b)) f1 := 0; Hypergeometric_function(a, amcpone, apbmcpone, onemonesz, eps, fa, n, f1); // 2F1(a, a-c+1, a+b-c+1, 1-1/z) f2 := 0; Hypergeometric_function(cma, onema, cmambpone, onemonesz, eps, fb, n, f2); // 2F1(c-a, 1-a, c-a-b+1, 1-1/z) mpc_mul(g1, zhma, tmc); mpc_mul(tmc, fa, afa); // A1 mpc_mul(g2, onmzhcab, tmc); mpc_mul(tmc, zhamc, tmc); mpc_mul(tmc, fb, afb); // A2 mpc_add(afa, afb, ans); // ans = A1 + A2 if f1 <> 0 then f := f1; if f2 <> 0 then f := f2; mpc_clear5(ga, gb, gc, gcmamb, gapbmc); mpc_clear5(ma, amcpone, apbmcpone, cma, onema); mpc_clear5(onemz, cmamb, chsa, amc, onemonesz); mpc_clear5(cone, onmzhcab, zhma, tmc, cmambpone); mpc_clear5(cmb, gcmb, gcma, apbmc, zhamc); mpc_clear4(fa, fb, g1, g2); mpc_clear2(afa, afb); end; // x <= 0 で 整数なら result = true function zeroorounder(xx: mp_complex): boolean; label EXT; var xc : mp_float; x : mp_complex; begin mpf_init(xc); mpc_init(x); // xxの値を変化させるxxの元値を解放不可になるのでコピーして使用します mpc_copy(xx, x); result := false; if not mpf_is0(x.im) then goto EXT; // frac 桁落ち対策 mpf_read_decimal(x.re, PAnsiChar(mpf_decimal(x.re, 200) + #00)); mpf_frac(x.re, xc); if not mpf_is0(xc) then goto EXT; if mpf_is_le(x.re, zero) then result := true; EXT: mpf_clear(xc); mpc_clear(x); end; // +∞ = +1; // -∞ = -1; // 整数 <=0 奇数 -1 偶数 +1 function infpm(x: mp_float): integer; var tmp , fra: mp_float; begin mpf_init2(tmp, fra); result := 0; mpf_div(x, two, tmp); mpf_frac(tmp, fra); if mpf_is0(fra) then result := 1 else result := -1; mpf_clear2(tmp, fra); 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; // {Γ(a)Γ(b)}/{Γ(a)Γ(a)} function gammaCalc(a, b, c, d: mp_complex; var ans: mp_complex): integer; var af, bf, cf, df : integer; ac, bc, cc, dc : mp_complex; tmp : mp_complex; dlt : mp_float; xp, xm : mp_complex; procedure gmmadlt(var x, g : mp_complex); var n : integer; dn : double; begin dn := mpf_todouble(x.re); n := round(dn); n := abs(n); mpf_set0(g.im); factorialMul(n, g.re); mpc_inv(g, g); if n mod 2 <> 0 then mpc_chs(g, g); end; { procedure gmmadlt(var x, g : mp_complex); begin mpc_add_mpf(x, dlt, xp); mpc_sub_mpf(x, dlt, xm); gammaMul(xp, g); gammaMul(xm, tmp); mpc_chs(tmp, tmp); mpc_add(g, tmp, tmp); mpc_div_mpf(tmp, two, g); end; } begin mpc_init4(ac, bc, cc, dc); mpc_init(tmp); mpf_init(dlt); mpc_init2(xm, xp); result := 1; mpc_set0(ans); af := 0; bf := 0; cf := 0; df := 0; mpf_read_decimal(dlt, PAnsiChar(ansistring('1E-55' + #00))); if zeroorounder(a) then af := infpm(a.re); if zeroorounder(b) then bf := infpm(b.re); if zeroorounder(c) then cf := infpm(c.re); if zeroorounder(d) then df := infpm(d.re); if af <> 0 then gmmadlt(a, ac) // Γ(a) else gammaMul(a, ac); if bf <> 0 then gmmadlt(b, bc) // Γ(b) else gammaMul(b, bc); if cf <> 0 then gmmadlt(c, cc) // Γ(c) else gammaMul(c, cc); if df <> 0 then gmmadlt(d, dc) // Γ(d) else gammaMul(d, dc); if (abs(af) + abs(bf)) = (abs(cf) + abs(df)) then begin mpc_mul(ac, bc, ans); // Γ(a) * Γ(b) mpc_div(ans, cc, ans); // Γ(a) * Γ(b) / Γ(c) mpc_div(ans, dc, ans); // Γ(a) * Γ(b) / Γ(c) / Γ(d) result := 0; end; if (abs(af) + abs(bf)) > (abs(cf) + abs(df)) then begin // result := 0; // mpf_copy(maxmpf, ans.re); result := 1; mpc_set1(ans); end; if (abs(af) + abs(bf)) < (abs(cf) + abs(df)) then begin result := 0; mpc_set0(ans); end; mpc_clear4(ac, bc, cc, dc); mpc_clear(tmp); mpf_clear(dlt); mpc_clear2(xm, xp); end; // z > 1 Γ値のみ補正計算 // c < 0 b < 0 c < 0 は x < 1 のルーチンで計算されます procedure Hypergeometric_function_of_NO_delt(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); label EXT; var g1, g2 , mzhma, mzhmb, mz : mp_complex; onepamc, onepbmc, onepamb, onemapb : mp_complex; ma, mb, sa, sb, tmc : mp_complex; onesz, fa, fb : mp_complex; amb, bma, cma, cmb : mp_complex; inf1, inf2, f1, f2 : integer; begin mpc_init5(g1, g2, mzhma, mzhmb, mz); mpc_init4(onepamc, onepbmc, onepamb, onemapb); mpc_init5(ma, mb, sa, sb, tmc); mpc_init3(onesz, fa, fb); mpc_init4(amb, bma, cma, cmb); mpc_sub(a, b, amb); // a-b mpc_sub(b, a, bma); // b-a mpc_sub(c, a, cma); // c-a mpc_sub(c, b, cmb); // c-b mpc_chs(z, mz); // -z mpc_chs(a, ma); // -a mpc_chs(b, mb); // -b mpc_powa(mz, ma, mzhma); // (-z)^-a mpc_powa(mz, mb, mzhmb); // (-z)^-b mpc_add_mpf(a, one, tmc); mpc_sub(tmc, c, onepamc); // 1+a-c mpc_sub(tmc, b, onepamb); // 1+a-b mpc_add_mpf(b, one, tmc); mpc_sub(tmc, c, onepbmc); // 1+b-c mpc_sub(tmc, a, onemapb); // 1-a+b mpc_inv(z, onesz); // 1/z // 1+a-b <=0 の時 無限大の判定 2F1(a, b, c, z)の c<=0 による判定 if zeroorounder(onepamb) then begin // 1+a-c <=0 の時 if zeroorounder(onepamc) then begin // 1+a-c <= 1+a-b if mpf_is_le(onepamc.re, onepamb.re) then begin f := 128; // 無限大フラグセット goto EXT; end; end // 1+a-c > 0 の時 else begin f := 128; // 無限大フラグセット goto EXT; end; end; // 1-a+b <= 0 の時 無限大の判定 2F1(a, b, c, z)の c<=0 による判定 if zeroorounder(onemapb) then begin // 1+b-c <=0 の時 if zeroorounder(onepbmc) then begin // 1+b-c <= 1-a+b if mpf_is_le(onepbmc.re, onemapb.re) then begin f := 128; // 無限大フラグセット goto EXT; end; end // 1+b-c > 0 の時 else begin f := 128; // 無限大フラグセット goto EXT; end; end; f1 := 0; Hypergeometric_function(a, onepamc, onepamb, onesz, eps, sa, n, f1); // 2F1(a, 1+a-c. 1+a-b, 1/z) form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(sa.re, 20))); form1.memo1.Lines.Append(' ' + string(' ' + mpf_decimal(sa.im, 20)) + ' i'); // BR := false; // 割込みフラグ解除 f2 := 0; Hypergeometric_function(b, onepbmc, onemapb, onesz, eps, sb, n, f2); // 2F1(b, 1+b-c, 1-a+b, 1/z) form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(sb.re, 20))); form1.memo1.Lines.Append(' ' + string(' ' + mpf_decimal(sb.im, 20)) + ' i'); inf1 := gammaCalc(c, bma, b, cma, g1); // form1.memo1.Lines.Append(' Γ1 ' + string(' ' + mpf_decimal(g1.re, 20))); // form1.memo1.Lines.Append(' Γ1 ' + string(' ' + mpf_decimal(g1.im, 20)) +' i'); mpc_mul(g1, mzhma, tmc); // Γ(c) * Γ(b-a) / Γ(b) / Γ(c-a) * (-z)^-a mpc_mul(tmc, sa, fa); // Γ(c) * Γ(b-a) / Γ(b) / Γ(c-a) * (-z)^-a * 2F1(a) inf2 := gammaCalc(c, amb, a, cmb, g2); // form1.memo1.Lines.Append(' Γ2 ' + string(' ' + mpf_decimal(g2.re, 20))); // form1.memo1.Lines.Append(' Γ2 ' + string(' ' + mpf_decimal(g2.im, 20)) +' i'); mpc_mul(g2, mzhmb, tmc); // Γ(c) * Γ(a-b) / Γ(a) / Γ(c-b) * (-z)^-b mpc_mul(tmc, sb, fb); // Γ(c) * Γ(a-b) / Γ(a) / Γ(c-b) * (-z)^-b * 2F1(b) if (f1 = -1) or (f1 = 1) then f := 16; if (f2 = -1) or (f2 = 1) then f := 16; if inf1 =1 then f := 16; if inf2 =1 then f := 16; if (f1 = -2) or (f2 = -2) then f := -2; mpc_add(fa, fb, ans); if mpc_is0(amb) then // a = b 時は fa=fbとなり 解はfa 又は fbとなります mpc_div_mpf(ans, two, ans); EXT: mpc_clear5(g1, g2, mzhma, mzhmb, mz); mpc_clear4(onepamc, onepbmc, onepamb, onemapb); mpc_clear5(ma, mb, sa, sb, tmc); mpc_clear3(onesz, fa, fb); mpc_clear4(amb, bma, cma, cmb); end; // z > 1 近似計算 // c < 0 b < 0 c < 0 は x < 1 のルーチンで計算されます procedure Hypergeometric_function_of_first_kind(a, b, c, amb, bma, cma, cmb, aeb, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); label EXT; var ga, gb, gc, gamb, gbma : mp_complex; gcma, gcmb, mzhma, mzhmb, mz : mp_complex; onepamc, onepbmc, onepamb, onemapb : mp_complex; ma, mb, sa, sb, tmc : mp_complex; onesz, fa, fb : mp_complex; f1, f2 : integer; begin // if mpc_is0(a) or mpc_is0(b) then begin // mpc_set1(ans); // exit; // end; mpc_init5(ga, gb, gc, gamb, gbma); mpc_init5(gcma, gcmb, mzhma, mzhmb, mz); mpc_init4(onepamc, onepbmc, onepamb, onemapb); mpc_init5(ma, mb, sa, sb, tmc); mpc_init3(onesz, fa, fb); gammaMul(a, ga); // Γ(a) gammaMul(b, gb); // Γ(b) gammaMul(c, gc); // Γ(c) gammaMul(amb, gamb); // Γ(a-b) gammaMul(bma, gbma); // Γ(b-a) gammaMul(cma, gcma); // Γ(c-a) gammaMul(cmb, gcmb); // Γ(c-b) mpc_chs(z, mz); // -z mpc_chs(a, ma); // -a mpc_chs(b, mb); // -b mpc_powa(mz, ma, mzhma); // (-z)^-a mpc_powa(mz, mb, mzhmb); // (-z)^-b mpc_set1(tmc); // mpc_sub(tmc, cma, onepamc); // 1+a-c mpc_sub(tmc, cmb, onepbmc); // 1+b-c mpc_add(tmc, amb, onepamb); // 1+a-b mpc_add(tmc, bma, onemapb); // 1-a+b mpc_inv(z, onesz); // 1/z Hypergeometric_function(a, onepamc, onepamb, onesz, eps, sa, n, f1); // 2F1(a, 1+a-c. 1+a-b, 1/z) form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(sa.re, 20))); form1.memo1.Lines.Append(' ' + string(' ' + mpf_decimal(sa.im, 20)) + ' i'); // BR := false; // 割込みフラグ解除 Hypergeometric_function(b, onepbmc, onemapb, onesz, eps, sb, n, f2); // 2F1(b, 1+b-c, 1-a+b, 1/z) form1.memo1.Lines.Append('loop = ' + intTostr(n) + string(' ' + mpf_decimal(sb.re, 20))); form1.memo1.Lines.Append(' ' + string(' ' + mpf_decimal(sb.im, 20)) + ' i'); mpc_mul(gc, gbma, tmc); // Γ(c) * Γ(b-a) mpc_div(tmc, gb, tmc); // Γ(c) * Γ(b-a) / Γ(b) mpc_div(tmc, gcma, tmc); // Γ(c) * Γ(b-a) / Γ(b) / Γ(c-a) // form1.memo1.Lines.Append(' Γ1 ' + string(' ' + mpf_decimal(tmc.re, 20))); // form1.memo1.Lines.Append(' Γ1 ' + string(' ' + mpf_decimal(tmc.im, 20)) +' i'); mpc_mul(tmc, mzhma, tmc); // Γ(c) * Γ(b-a) / Γ(b) / Γ(c-a) * (-z)^-a mpc_mul(tmc, sa, fa); // Γ(c) * Γ(b-a) / Γ(b) / Γ(c-a) * (-z)^-a * 2F1(a) mpc_mul(gc, gamb, tmc); // Γ(c) * Γ(a-b) mpc_div(tmc, ga, tmc); // Γ(c) * Γ(a-b) / Γ(a) mpc_div(tmc, gcmb, tmc); // Γ(c) * Γ(a-b) / Γ(a) / Γ(c-b) // form1.memo1.Lines.Append(' Γ2 ' + string(' ' + mpf_decimal(tmc.re, 20))); // form1.memo1.Lines.Append(' Γ2 ' + string(' ' + mpf_decimal(tmc.im, 20)) +' i'); mpc_mul(tmc, mzhmb, tmc); // Γ(c) * Γ(a-b) / Γ(a) / Γ(c-b) * (-z)^-b mpc_mul(tmc, sb, fb); // Γ(c) * Γ(a-b) / Γ(a) / Γ(c-b) * (-z)^-b * 2F1(b) mpc_add(fa, fb, ans); // Σ // a = b 時は fa=fbとなり 解はfa 又は fbとなりますので1/2にします。 if mpc_is0(aeb) then mpc_div_mpf(ans, two, ans); if (f1 = 1) or (f1 = 2) then f := 16; if (f2 = 1) or (f2 = 2) then f := 16; EXT: mpc_clear5(ga, gb, gc, gamb, gbma); mpc_clear5(gcma, gcmb, mzhma, mzhmb, mz); mpc_clear4(onepamc, onepbmc, onepamb, onemapb); mpc_clear5(ma, mb, sa, sb, tmc); mpc_clear3(onesz, fa, fb); end; // z=1 超幾何定理計算 procedure Hypergeometric_function_zeq1(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); label EXT; var cmamb, cma, cmb, tmp: mp_complex; gc, gcma, gcmb, gcmamb : mp_complex; begin mpc_init4(cmamb, cma, cmb, tmp); mpc_init4(gc, gcma, gcmb, gcmamb); f := 0; n := 0; mpc_sub(c, a, cmamb); // c - a mpc_sub(cmamb, b, cmamb); // c - a - b if mpf_is_le(cmamb.re, zero) then begin // Re(c-a-b) <= 0 f := 8; // ∞±符号計算フラグ mpc_set0(ans); goto EXT; end; // Re(c-a-b) > 0 時の計算 mpc_sub(c, a, cma); mpc_sub(c, b, cmb); gammaMul(c, gc); // Γ(c) gammaMul(cma, gcma); // Γ(c-a) gammaMul(cmb, gcmb); // Γ(c-b) gammaMul(cmamb, gcmamb); // Γ(c-a-b) mpc_mul(gc, gcmamb, tmp); // Γ(c) * Γ(c-a-b) mpc_div(tmp, gcma, tmp); // Γ(c) * Γ(c-a-b) / Γ(c-a) mpc_div(tmp, gcmb, ans); // Γ(c) * Γ(c-a-b) / Γ(c-a) / Γ(c-b) EXT: mpc_clear4(cmamb, cma, cmb, tmp); mpc_clear4(gc, gcma, gcmb, gcmamb); end; // z < -1 -1>z>-∞ 実数部の値が-1より小さい場合の専用計算 procedure Hypergeometric_function_mz(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); var onemz, bma, cma, cmb, ambp1 : mp_complex; amb, bmap1 : mp_complex; gc, gbma, gb, gcma : mp_complex; gamb, ga, gcmb : mp_complex; onec, tmp, tmc, invz : mp_complex; g1, g2, onemzhma, onemzhmb : mp_complex; tmf : mp_float; f1, f2 : integer; begin mpc_init5(onemz, bma, cma, cmb, ambp1); mpc_init2(amb, bmap1); mpc_init4(gc, gbma, gb, gcma); mpc_init3(gamb, ga, gcmb); mpc_init4(onec, tmp, tmc, invz); mpc_init4(g1, g2, onemzhma, onemzhmb); mpf_init(tmf); mpc_set1(onec); mpc_sub(onec, z, onemz); mpc_inv(onemz, invz); mpc_sub(b, a, bma); mpc_sub(c, a, cma); mpc_sub(a, b, amb); mpc_sub(c, b, cmb); mpc_add(amb, onec, ambp1); mpc_add(bma, onec, bmap1); mpc_inv(onemz, invz); gammaMul(a, ga); // Γ(a) gammaMul(b, gb); // Γ(b) gammaMul(c, gc); // Γ(c) gammaMul(bma, gbma); // Γ(b-a) gammaMul(cma, gcma); // Γ(c-a) gammaMul(amb, gamb); // Γ(a-b) gammaMul(cmb, gcmb); // Γ(c-b) mpc_chs(a, tmp); // -a mpc_powa(onemz, tmp, onemzhma); // (1-z)^-a mpc_chs(b, tmp); // -a mpc_powa(onemz, tmp, onemzhmb); // (1-z)^-b mpc_div(gc, gb, tmp); // Γ(c)/Γ(b) mpc_mul(tmp, gbma, tmp); // Γ(c)/Γ(b) * Γ(bma) mpc_div(tmp, gcma, g1); // Γ(c)/Γ(b) * Γ(bma)/Γ(c-a) mpc_div(gc, ga, tmp); // Γ(c)/Γ(a) mpc_mul(tmp, gamb, tmp); // Γ(c)/Γ(b) * Γ(amb) mpc_div(tmp, gcmb, g2); // Γ(c)/Γ(b) * Γ(amb)/Γ(c-b) f1 := 0; Hypergeometric_function(a, cmb, ambp1, invz, eps, tmc, n, f1); mpc_mul(tmc, onemzhma, tmc); // F*(1-z)^^a mpc_mul(tmc, g1, tmc); // F*(1-z)^a * Γ/Γ f2 := 0; Hypergeometric_function(b, cma, bmap1, invz, eps, ans, n, f2); mpc_mul(ans, onemzhmb, ans); // F*(1-z)^b mpc_mul(ans, g2, ans); // F*(1-z)^b * Γ/Γ mpc_add(tmc, ans, ans); if mpc_is0(amb) then // a=bの時は、1/2 mpc_div_mpf(ans, two, ans); if f1 <> 0 then f := f1; if f2 <> 0 then f := f2; mpc_clear5(onemz, bma, cma, cmb, ambp1); mpc_clear2(amb, bmap1); mpc_clear4(gc, gbma, gb, gcma); mpc_clear3(gamb, ga, gcmb); mpc_clear4(onec, tmp, tmc, invz); mpc_clear4(g1, g2, onemzhma, onemzhmb); mpf_clear(tmf); end; // 0.5 < z < 1 procedure Hypergeometric_function_harftone(a, b, c, z : mp_complex; eps : mp_float; var ans : mp_complex; var n, f : integer); var cmamb, cma, cmb, apbmcp1 : mp_complex; apbmc, cmambp1, tmc, onemz, onemzhcab : mp_complex; ga, gb, gc, g1, g2: mp_complex; gcmamb, gapbmc, tmp, onec : mp_complex; gcma, gcmb : mp_complex; tmf : mp_float; f1, f2: integer; begin mpc_init4(cmamb, cma, cmb, apbmcp1); mpc_init5(apbmc, cmambp1, tmc, onemz, onemzhcab); mpc_init5(ga, gb, gc, g1, g2); mpc_init4(gcmamb, gapbmc, tmp, onec); mpc_init2(gcma, gcmb); mpf_init(tmf); mpc_set1(onec); mpc_sub(c, a, cma); // c-a mpc_sub(c, b, cmb); // c-b mpc_sub(cma, b, cmamb); // c-a-b mpc_add(a, b, tmp); // a+b mpc_sub(tmp, c, apbmc); // a+b-c mpc_add(apbmc, onec, apbmcp1); // a+b-c+1 mpc_add(cmamb, onec, cmambp1); // c-a-b+1 mpc_sub(onec, z, onemz); // 1-z gammaMul(a, ga); // Γ(a) gammaMul(b, gb); // Γ(b) gammaMul(c, gc); // Γ(c) gammaMul(cma, gcma); // Γ(c-a) gammaMul(cmb, gcmb); // Γ(c-b) gammaMul(cmamb, gcmamb); // Γ(c-a-b) gammaMul(apbmc, gapbmc); // Γ(a+b-c) mpc_div(gc, gcma, tmp); // Γ(c)/Γ(c-a) mpc_mul(tmp, gcmamb, tmp); // Γ(c)/Γ(c-a)*Γ(c-a-b) mpc_div(tmp, gcmb, g1); // g1 = Γ(c)/Γ(c-a)*Γ(c-a-b)/Γ(c-b) mpc_div(gc, ga, tmp); // Γ(c)/Γ(a) mpc_mul(tmp, gapbmc, tmp); // Γ(c)/Γ(a)*Γ(a+b-c) mpc_div(tmp, gb, g2); // g2 = Γ(c)/Γ(a)*Γ(a+b-c)/Γ(b) mpc_powa(onemz, cmamb, onemzhcab); // (1-z)^(c-a-b) f1 := 0; Hypergeometric_function(a, b, apbmcp1, onemz, eps, tmc, n, f1); mpc_mul(tmc, g1, tmc); form1.memo1.Lines.Append('1-z tmc.re ' + string(' ' + mpf_decimal(tmc.re, 20))); f2 := 0; Hypergeometric_function(cma, cmb, cmambp1, onemz, eps, ans, n, f2); mpc_mul(ans, g2, ans); mpc_mul(ans, onemzhcab, ans); form1.memo1.Lines.Append('1-z ans.re ' + string(' ' + mpf_decimal(ans.re, 20))); mpc_add(ans, tmc, ans); form1.memo1.Lines.Append('1-z ans+tmc ' + string(' ' + mpf_decimal(ans.re, 20))); if f1 <> 0 then f := f1; if f2 <> 0 then f := f2; mpc_clear4(cmamb, cma, cmb, apbmcp1); mpc_clear5(apbmc, cmambp1, tmc, onemz, onemzhcab); mpc_clear5(ga, gb, gc, g1, g2); mpc_clear4(gcmamb, gapbmc, tmp, onec); mpc_clear2(gcma, gcmb); mpf_clear(tmf); end; // arcsin計算 procedure TForm1.arsinBtnClick(Sender: TObject); var ch : integer; ars : double; tmp0, tmp1, tmp2 : mp_float; begin val(arsinEdit.Text, ars, ch); if ch <> 0 then begin application.MessageBox('arsin の値に間違いがあります。','注意',0); exit; end; if abs(ars) > 1 then begin application.MessageBox('arsin の値は±1迄にして下さい。','注意',0); exit; end; mpf_init3(tmp0, tmp1, tmp2); mpf_read_decimal(tmp0, PAnsiChar(ansistring(arsinEdit.Text + #00))); if (abs(ars) < 0.75) or (abs(ars) = 1) then TN := 2 else TN := 1; aedit.Text := '0.5'; if (abs(ars) < 0.75) or (abs(ars) = 1) then bedit.text := '0.5' else bedit.text := '1'; cedit.text := '1.5'; if (abs(ars) < 0.75) or (abs(ars) = 1) then zedit.text := arsinEdit.Text else begin mpf_mul(tmp0, tmp0, tmp1); // ars * ars mpf_sub(one, tmp1, tmp2); // 1 - ars * ars mpf_sqrt(tmp2, tmp1); // sqrt(1 - ars * ars) mpf_div(tmp0, tmp1, tmp2); // ars / sqrt(1 - ars * ars) // ars := ars / sqrt(1 - ars * ars); zedit.Text := string(mpf_decimal(tmp2, 50)); // zedit.text := floattostr(ars); end; iaedit.Text := '0'; ibedit.text := '0'; icedit.text := '0'; izedit.text := '0'; DS := 2; mpf_clear3(tmp0, tmp1, tmp2); Button1Click(nil); end; // arctan計算 procedure TForm1.artanBtnClick(Sender: TObject); var ch : integer; art : double; tmp0, tmp1, tmp2 : mp_float; begin val(ArtanedEdit.Text, art, ch); if ch <> 0 then begin application.MessageBox('artan の値に間違いがあります。','注意',0); exit; end; mpf_init3(tmp0, tmp1, tmp2); mpf_read_decimal(tmp0, PAnsiChar(ansistring(ArtanedEdit.Text + #00))); if (abs(art) < 0.75) or (abs(art) > 1.4) then TN := 1 else TN := 2; aedit.Text := '0.5'; if (abs(art) < 0.75) or (abs(art) > 1.4) then bedit.text := '1' else bedit.text := '0.5'; cedit.text := '1.5'; if (abs(art) < 0.75) or (abs(art) > 1.4) then zedit.text := artanededit.Text else begin mpf_mul(tmp0, tmp0, tmp1); // art * art mpf_add(tmp1, one, tmp2); // 1 + art * art mpf_sqrt(tmp2, tmp1); // sqrt(1 + art * art) mpf_div(tmp0, tmp1, tmp2); // art / sqrt(1 + art * art) zedit.Text := string(mpf_decimal(tmp2, 50)); end; iaedit.Text := '0'; ibedit.text := '0'; icedit.text := '0'; izedit.text := '0'; DS := 1; mpf_clear3(tmp0, tmp1, tmp2); Button1Click(nil); end; // 計算打切り procedure TForm1.BitBtn1Click(Sender: TObject); begin BR := True; button1.Enabled := True; end; // 数値の後ろのゼロ消去 function ZeroErase(s : string): string; const EP = 'E'; ZC = '0'; dt = '.'; var c : char; i, j, k, l : integer; begin l := length(s); j := 1; for i := 1 to l do begin c := s[i]; if c = EP then begin j := i - 1; break; end; j := i; end; result := ''; if j < l then begin for i := l downto j + 1 do result := s[i] + result; end; K := 1; for i := j downto 1 do begin c := s[i]; if c <> ZC then begin k := i; if c = DT then k := k + 1; break; end; end; for i := k downto 1 do result := s[i] + result; end; // 入力チェック function TForm1.inputcheck(var pre: integer): boolean; var chd: double; ch : integer; begin result := false; val(aedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('a の値に間違いがあります。','注意',0); exit; end; val(bedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('b の値に間違いがあります。','注意',0); exit; end; val(cedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('c の値に間違いがあります。','注意',0); exit; end; val(zedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('z の値に間違いがあります。','注意',0); exit; end; val(iaedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('a の虚数値に間違いがあります。','注意',0); exit; end; val(ibedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('b の虚数値に間違いがあります。','注意',0); exit; end; val(icedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('c の虚数値に間違いがあります。','注意',0); exit; end; val(izedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('z の虚数値に間違いがあります。','注意',0); exit; end; val(precisionedit.Text, pre, ch); if ch <> 0 then begin application.MessageBox('有効桁数 の値に間違いがあります。','注意',0); exit; end; if pre < 10 then begin application.MessageBox('有効桁数 は10以上にして下さい。','注意',0); exit; end; val(epsedit.Text, chd, ch); if ch <> 0 then begin application.MessageBox('収束判定値に間違いがあります。','注意',0); exit; end; if chd < 1e-100 then begin application.MessageBox('収束判定値は1e-100より大きくして下さい。','注意',0); exit; end; if chd > 1e-1 then begin application.MessageBox('収束判定値は1e-1より小さく下さい。','注意',0); exit; end; result := true; end; // 計算開始 procedure TForm1.Button1Click(Sender: TObject); label EXT, TN0E, ZG01, ZG02, ZG03; var a, b, c, z, ans : mp_complex; mz2, tans, dat, tmc, onec : mp_complex; n, pre : integer; f : integer; eps, az: mp_float; amb, bma, cma, cmb, aeb : mp_complex; d1, d2 : mp_float; XLE, EUF : boolean; ap, bp, cp : mp_complex; am, bm, cm : mp_complex; zback, aback, zbackp : mp_complex; euans : mp_complex; dnine, tmf, absz : mp_float; df : integer; // 1E-35 以下はゼロにセット procedure zerounde1em35(var ans: mp_complex); begin mpf_read_decimal(d1, PAnsiChar(ansistring('1E-35' + #00))); mpf_abs(ans.re, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.re); mpf_abs(ans.im, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.im); end; // 1E-42 以下はゼロにセット procedure zerounde1em42(var ans: mp_complex); begin mpf_read_decimal(d1, PAnsiChar(ansistring('1E-42' + #00))); mpf_abs(ans.re, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.re); mpf_abs(ans.im, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.im); end; // 1E-52 以下はゼロにセット procedure zerounde1em52(var ans: mp_complex); begin mpf_read_decimal(d1, PAnsiChar(ansistring('1E-52' + #00))); mpf_abs(ans.re, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.re); mpf_abs(ans.im, tmf); if mpf_is_lt(tmf, d1) then mpf_set0(ans.im); end; // 表示桁数を収束反対値に合わせる procedure epslimit(var ans: mp_complex); const mch = '-'; mcd = '.'; var str, nstr : string; lng, pc, flng : integer; begin str := epsedit.Text; lng := length(str); pc := pos(mch, str); nstr := copy(str, pc+1, lng - pc); // 1e-(xx) nstr = xx flng := strToint(nstr); str := (ZeroErase(string(mpf_decimal(ans.re, 100)))); pc := pos(mcd, str); // xx.zzzzzzzzzz の'.'の位置 lng := pc + flng - 1; // 表示桁数 nstr := (ZeroErase(string(mpf_decimal(ans.re, lng)))); // ans.re表示桁数文字変換 str := (ZeroErase(string(mpf_decimal(ans.im, lng)))); // ans.im表示桁数文字変換 mpf_read_decimal(ans.re, PAnsiChar(ansistring(nstr + #00))); // re 文字数値変換 mpf_read_decimal(ans.im, PAnsiChar(ansistring(str + #00))); // im 文字数値変換 end; begin // 入力チェック if not inputcheck(pre) then exit; button1.Enabled := False; BR := false; // 停止フラグリセット mpf_set_default_prec(pre); // 有効桁数セット mpc_init5(a, b, c, z, ans); mpf_init3(eps, az, tmf); mpc_init5(mz2, tans, dat, tmc, onec); mpc_init5(amb, bma, cma, cmb, aeb); mpf_init4(d1, d2, dnine, absz); mpc_init4(ap, bp, cp, zbackp); mpc_init5(am, bm, cm, zback, aback); mpc_init(euans); memo1.Clear; memo1.Lines.Append('計算中'); application.ProcessMessages; // 入力値 多倍長 mpf_read_decimal(a.re, PAnsiChar(ansistring(aedit.Text + #00))); mpf_read_decimal(b.re, PAnsiChar(ansistring(bedit.Text + #00))); mpf_read_decimal(c.re, PAnsiChar(ansistring(cedit.Text + #00))); mpf_read_decimal(z.re, PAnsiChar(ansistring(zedit.Text + #00))); mpf_read_decimal(a.im, PAnsiChar(ansistring(iaedit.Text + #00))); mpf_read_decimal(b.im, PAnsiChar(ansistring(ibedit.Text + #00))); mpf_read_decimal(c.im, PAnsiChar(ansistring(icedit.Text + #00))); mpf_read_decimal(z.im, PAnsiChar(ansistring(izedit.Text + #00))); mpf_read_decimal(eps, PAnsiChar(ansistring(epsedit.Text + #00))); // 値設定 mpc_set1(onec); // c1 mpc_arg(z, tmf); mpf_abs(tmf, tmf); memo1.Lines.Append(ZeroErase(string('|Arg(z)| =' + mpf_decimal(tmf, 40)))); mpc_sub(onec, z, tmc); mpc_arg(tmc, tmf); mpf_abs(tmf, tmf); memo1.Lines.Append(ZeroErase(string('|Arg(1-z)| =' + mpf_decimal(tmf, 40)))); mpc_abs(z, absz); // |z| // オイラー積分 2000等分 // R(c)>R(a)>0 |z|<1 // aとcの虚数部の値の差が大きくなる値誤差が大きくなります。 EUF := false; if checkbox8.Checked = true then begin mpc_sub(c, a, am); // RE(a) >= 1 RE(c-a) >=1 |z|<1 if mpf_is_ge(a.re, one) and mpf_is_ge(am.re, one) and mpf_is_lt(absz, one) then begin memo1.Lines.Append('オイラー積分中暫くお待ちください。'); Euler_integral(a, b, c, z, euans); EUF := true; end; end; // a>b なら入れ替え mpc_abs(a, ap.re); mpc_abs(b, bp.re); if mpf_is_gt(ap.re, bp.re) then begin mpc_copy(a, tmc); mpc_copy(b, a); mpc_copy(tmc, b); end; f := 0; n := 0; DS := 0; mpc_set0(ans); mpf_abs(z.re, az); // |z.re| mpc_copy(z, zback); // z -> zback { memo1.Lines.Append(ZeroErase(string('|z| ' + mpf_decimal(absz, 40)))); mpc_set1(tmc); mpc_sub(tmc, z, ap); // 1-z mpc_abs(ap, tmf); memo1.Lines.Append(ZeroErase(string('1-z ' + mpf_decimal(tmf, 40)))); mpc_div(tmc, z, ap); // 1/z mpc_abs(ap, tmf); memo1.Lines.Append(ZeroErase(string('1/z ' + mpf_decimal(tmf, 40)))); mpc_sub(tmc, ap, am); // 1-1/z mpc_abs(am, tmf); memo1.Lines.Append(ZeroErase(string('1-1/z ' + mpf_decimal(tmf, 40)))); mpc_sub(z, tmc, ap); // z-1 mpc_div(z, ap, am); // z/(z-1); mpc_abs(am, tmf); memo1.Lines.Append(ZeroErase(string('z/(z-1) ' + mpf_decimal(tmf, 40)))); mpc_sub(tmc, am, bp); // 1-z/(z-1) mpc_abs(bp, tmf); memo1.Lines.Append(ZeroErase(string('1-z/(z-1) ' + mpf_decimal(tmf, 40)))); goto EXT; } XLE := false; // 近似計算フラグリセット if TN = 0 then begin // 通常計算 // z = 0 if mpc_is0(z) then begin memo1.Lines.Append('zの値がゼロです'); memo1.Lines.Append(' 1.0'); memo1.Lines.Append(' +0.0i'); goto EXT; end; //-------------------------------------------------------------------------- // z = 1 超幾何定理計算 if mpc_is1(z) then begin Hypergeometric_function_zeq1(a, b, c, z, eps, ans, n, f); Memo1.Lines.Append(' z = 1 専用'); if f = 0 then begin Memo1.Lines.Append(' Re(c - a - b) > 0'); goto TN0E; end; // ∞符号設定 1から0.3 減算 0.7・・ 計算し∞±符号設定 if f = 8 then begin mpf_set_dbl(d1, 0.3); mpc_sub_mpf(z, d1, z); goto ZG03; end; end; // |Z| > 1 実数部0の時 if mpf_is_gt(absz, one) and mpf_is0(z.re) then begin goto ZG02; end; // ------------------------------------------------------------------------- // 整数 (c <= 0 or a <= 0 or b <= 0) a,b,cの何れか又は全部が0を含み負の整数の場合 // zの値に関係なく通常の超幾何関数計算 Z<1の計算使用します。 二番目に確認 // /Z/の値が1より大きい場合aかb又はcがゼロになった時点で計算計算打ち切り // cが先か同時にゼロになった場合は±∞となります。 // |z|の値が1より小さい場合は、a,bが先にゼロになるか、cがゼロになった時 // a,bがゼロになってないかで計算がかわります。 if zeroorounder(c) or (zeroorounder(a) or zeroorounder(b)) then begin memo1.Lines.Append(' (整数a<=0) or (整数b<=0) or (整数c<=0) '); Hypergeometric_function(a, b, c, z, eps, ans, n, f); goto TN0E; end; //-------------------------------------------------------------------------- // z.re < -1 -1> z > -2 -∞でも可 -2迄に限定 デバッグの為 a, b, c に非整数がある場合 // 非整数だとpfaff変換が出来ないため 1/(1-z)で計算 mpf_frac(a.re, ap.re); mpf_frac(b.re, bp.re); mpf_frac(c.re, cp.re); df := 0; if not mpf_is0(ap.re) or not mpf_is0(bp.re) or not mpf_is0(cp.re) then df := 1; if not mpf_is0(a.im) or not mpf_is0(b.im) or not mpf_is0(c.im) then df := df or 2; if (checkbox5.Checked = true) and (df <> 0) then begin // z,re < and |z| >=1 and |z|<2 if mpf_is_lt(z.re, zero) and mpf_is_ge(absz, one) and mpf_is_lt(absz, two) then begin mpc_sub(a, b, amb); // a-b mpc_sub(b, a, bma); // b-q mpc_sub(c, a, cma); // c-a mpc_add(amb, onec, ap); // a-b + 1 mpc_add(bma, onec, bp); // b-a + 1 // a,b,c<=0 整数は前で処理されます // c以外で無限大になるか確認 df := 0; // ∞になるか確認 if zeroorounder(amb) then df := df or 1; // a-b if zeroorounder(bma) then df := df or 2; // b-a if zeroorounder(ap) then df := df or 4; // a-b+1 if zeroorounder(bp) then df := df or 8; // b-a+1 // a=b=c 確認 if mpc_is0(amb) and mpc_is0(cma) then df := 16; // Γ(0)/Γ(0) になります。 if (df = 0) or (df = 16) then begin Hypergeometric_function_mz(a, b, c, z, eps, ans, n, f ); memo1.Lines.Append('Z < -1 -1>z>-2'); memo1.Lines.Append('変換 1/(1-z) Δ無し'); goto TN0E; end else begin mpf_read_decimal(d1, PAnsiChar(ansistring('1E-20' + #00))); mpc_add_mpf(c, d1, cp); mpc_add_mpf(a, d1, ap); Hypergeometric_function_mz(ap, b, cp, z, eps, tmc, n, f); mpc_sub_mpf(c, d1, cm); mpc_sub_mpf(a, d1, am); Hypergeometric_function_mz(am, b, cm, z, eps, ans, n, f); mpc_add(tmc, ans, tmc); mpc_div_mpf(tmc, two, ans); memo1.Lines.Append('Z < -1 -1>z>-2'); memo1.Lines.Append('変換 1/(1-z) a,c±Δ計算'); XLE := true; goto TN0E; end; end; end; //-------------------------------------------------------------------------- // f(a, b, b, z) z<>1 mpc_sub(c, a, am); mpc_sub(c, b, bm); if (mpc_is0(am) or mpc_is0(bm)) and not mpc_is1(z) and not mpc_is0(z) then begin if mpc_is0(am) then begin mpc_copy(a, tmc); mpc_copy(b, a); mpc_copy(tmc, b); end; mpc_chs(a, ap); mpc_sub(onec, z, tmc); // 1-z mpc_powa(tmc, ap, ans); // (1-z)^-a memo1.Lines.Append('F(a, b, b, z) z <> 1'); goto TN0E; end; //-------------------------------------------------------------------------- // a,a+1,c or b+1,b,c z>1 時の近似計算 mpf_add(a.re, one, ap.re); // 1加算 減算時桁落ち対策 mpf_add(b.re, one, bp.re); // 1加算 if mpf_is_gt(ap.re, bp.re) then begin // ap > bp なら ap <-> bp mpc_copy(ap, tmc); mpc_copy(bp, ap); mpc_copy(tmc, bp); end; mpf_sub(bp.re, ap.re, cp.re); // re(a),re(b) 差分 差分は0～+ mpf_frac(a.re, tmf); // re(a)の小数部 // im(z) = 0 なら if (checkbox4.Checked = true) then begin mpf_read_decimal(cp.re, PAnsiChar(mpf_decimal(cp.re, 100) + #00)); // 桁落ち対策 差分1確認の為 100桁で丸め // re(a)-re(b) = 1 虚数部im(a)=0 |z{>1 実数部re(z)>=0 if mpf_is_eq(cp.re, one) and mpf_is0(a.im) and mpf_is_gt(absz, one) and mpf_is_ge(z.re, zero) then begin if not mpf_is0(tmf) then begin // 小数部がゼロでなかったら 整数でなかったら if mpf_is_le(absz, four) then begin // 1<=|z|<4 だったら DS := 32; // 専用近似計算フラグセット 1-1/z goto ZG01; // pfaff 変換不可 end else begin // |z|>=4 DS := 16; // 近似計算フラグセット goto ZG03; // pfaff 変換不要 end; end else begin // a, b が整数だったら if mpf_is_gt(absz, four) then begin // /z/>=4 なら DS := 16; // 近似計算フラグセット goto ZG03; // pfaff 変換不要 /z/>1の計算へ end else begin // 1<=|z|<=4 DS := 32; // 専用近似計算フラグセット z' = 1-1/z goto ZG01; // pfaff 変換不要 end; end; end else begin // a-b=1 im(a)=0 |z|>=1 re(z)<0 if mpf_is_eq(cp.re, one) and mpf_is0(a.im) and mpf_is_ge(absz, one) and mpf_is_lt(z.re, zero) then begin DS := 0; // 近似計算なし goto ZG02; // pfaff 変換 end; end; end; ZG01: //-------------------------------------------------------------------------- // a,a+1,c or b+1,b,c 1<z<2 時の近似計算 z'= 1-1/z 変換 // a,b,c 全て整数時は±Δ 非整数時は Δ0の時あり if DS = 32 then begin mpc_add(a, b, ap); // a+b mpc_sub(ap, c, cm); // a+b-c mpc_add(cm, onec, cp); // a+b-c+1 mpc_sub(c, a, cma); // c-a mpc_sub(cma, b, am); // c-a-b mpc_add(am, onec, bp); // c-a-b+1 df := 0; if zeroorounder(c) then df := df or 1; // c if zeroorounder(am) then df := df or 2; // c-a-b if zeroorounder(cp) then df := df or 4; // a+b-c+1 if zeroorounder(cm) then df := df or 8; // a+b-c if zeroorounder(bp) then df := df or 16; // c-a-b+1 // mpc_sub(a, b, amb); if (df = 0) then begin // Γ無限大が無かったら and a=b Hypergeometric_function_of_first_kind_aeb(a, b, c, z, eps, ans, n, f); memo1.Lines.Append('変換 1-1/z Δ無し'); DS := 0; // Δ無しなので近似計算ではありません。 goto TN0E; end; mpf_read_decimal(d1, PAnsiChar(ansistring('1E-20' + #00))); // mpf_read_decimal(d1, PAnsiChar(ansistring('431E-21' + #00))); mpc_sub_mpf(c, d1, cp); Hypergeometric_function_of_first_kind_aeb(a, b, cp, z, eps, tmc, n, f); mpc_add_mpf(c, d1, cm); Hypergeometric_function_of_first_kind_aeb(a, b, cm, z, eps, ans, n, f); mpc_add(tmc, ans, tmc); mpc_div_mpf(tmc, two, ans); memo1.Lines.Append('a,a+1,c or b+1,b,c 1<z<2'); memo1.Lines.Append('変換 1-1/z c±Δ計算'); XLE := true; goto TN0E; end; //-------------------------------------------------------------------------- // 0.5 < z <= 1 if (checkbox6.Checked = true) then begin mpf_frac(a.re, ap.re); mpf_frac(b.re, bp.re); mpf_frac(c.re, cp.re); mpf_set_dbl(tmf, 0.5); if mpf_is_gt(absz, tmf) and mpf_is_lt(absz, one) and mpf_is_gt(z.re, zero) and (not mpf_Is0(ap.re) or not mpf_Is0(bp.re) or not mpf_Is0(cp.re)) then begin DS := 0; mpc_add(a, b, ap); // a+b mpc_sub(ap, c, cm); // a+b-c mpc_add(cm, onec, amb); // a+b-c+1 mpc_sub(c, a, cma); // c-a mpc_sub(cma, b, am); // c-a-b mpc_add(am, onec, bp); // c-a-b+1 df := 0; if zeroorounder(am) then df := df or 1; // c-a-b if zeroorounder(amb) then df := df or 2; // a+b-c+1 if zeroorounder(cm) then df := df or 4; // a+b-c if zeroorounder(bp) then df := df or 8; // c-a-b+1 if zeroorounder(c) then df := df or 16; // c if df = 0 then begin Hypergeometric_function_harftone(a, b, c, z, eps, ans, n, f); memo1.Lines.Append('変換 1-z Δ無し'); goto TN0E; end; mpf_read_decimal(d1, PAnsiChar(ansistring('226476e-34' + #00))); mpc_add_mpf(c, d1, cp); Hypergeometric_function_harftone(a, b, cp, z, eps, tmc, n, f); mpc_sub_mpf(c, d1, cm); Hypergeometric_function_harftone(a, b, cm, z, eps, ans, n, f); mpc_add(tmc, ans, tmc); mpc_div_mpf(tmc, two, ans); memo1.Lines.Append('変換 1-z c±Δ計算'); XLE := true; goto TN0E; end; end; //-------------------------------------------------------------------------- // F a,a+1,c, z) or ( F b+1, b, c, z) |z| > 1 mpc_sub(b, a, am); mpc_sub(a, b, bm); if (mpc_is1(am) or mpc_is1(bm)) then begin if mpf_is_gt(absz, one) then begin goto ZG02; // pfaff変換へ end; end; //-------------------------------------------------------------------------- // F(1, 1, 2, z) z <> 1 mpc_sub(c, onec, tmc); // c - (1+0i) if mpc_is1(a) and mpc_is1(b) and mpc_is1(tmc) and not mpc_is1(z) and not mpc_is0(z) then begin mpc_sub(onec, z, tmc); // 1-z mpc_ln(tmc, tmc); // ln(1-z) mpc_div(tmc, z, ans); // ln(1-z)/z mpc_chs(ans, ans); // -ln(1-z)/z memo1.Lines.Append('F(1, 1, 2, z) z <> 1'); goto TN0E; end; //-------------------------------------------------------------------------- // a = b, c / a = 2 0.6<z<2 1-1/z 変換 // cが整数 c<=0 の場合は、先に処理されるので、aのみ近似計算 // a = b, c / a = 2 の時には ±Δ近似計算以外なし // a±Δ b±Δ c±Δ どれを使用しても結果は同じです // re(z) =0 /z/=1 は不可 mpc_add_mpf(a, one, ap); // 1加算 減算時桁落ち対策 mpc_add_mpf(b, one, bp); // 1加算 mpc_sub(bp, ap, bp); // a,b 差分 差分は0～+ mpf_set_dbl(d1, 0.6); if not mpc_is0(a) then mpc_div(c, a, cp); // c/a if (checkbox2.Checked = true) then // 0.6<|z| c/a=2 a=b z<>1 if mpf_is_gt(absz, d1) and mpf_is_eq(cp.re, two) and mpc_is0(bp) and not mpc_is1(z) then begin // zの範囲 0.6 < |z| < 2 if mpf_is_lt(absz, two) and not mpf_is1(absz) and mpf_is_gt(z.re, d1) then begin mpf_read_decimal(d1, PAnsiChar(ansistring('3E-19' + #00))); mpc_add_mpf(a, d1, ap); Hypergeometric_function_of_first_kind_aeb(ap, b, c, z, eps, tmc, n, f); mpc_sub_mpf(a, d1, am); Hypergeometric_function_of_first_kind_aeb(am, b, c, z, eps, ans, n, f); mpc_add(tmc, ans, ans); mpc_div_mpf(ans, two, ans); memo1.Lines.Append('a = b, c / a = 2 0.6 < z < 2'); memo1.Lines.Append('変換 1-1/z a±Δ計算'); DS := 32; // 近似計算フラグ低精度 XLE := true; // 近似計算フラグ goto TN0E; end; end; //-------------------------------------------------------------------------- // a = b, c <> a, 1<z<= 2 a,b 非整数 // Γ計算に±∞無し 近似計算の必要なし確認 非整数時はΔ無多い if not mpc_is0(a) then mpc_div(c, a, cp); // c/a if (checkbox2.Checked = true) and not mpf_is0(ap.re) and mpf_is0(bp.re) and mpf_is_gt(absz, one) and not mpf_is_eq(cp.re, two) and mpf_is_le(absz, two) then begin mpc_add(a, b, ap); // a+b mpc_sub(ap, c, ap); // a+b-c mpc_sub(c, a, am); // c-a mpc_sub(am, b, am); // c-a-b // a+b-c c-a-b c if not zeroorounder(ap) and not zeroorounder(am) and not zeroorounder(c) then begin mpc_add(ap, onec, ap); // a+b-c+1 mpc_add(am, onec, am); // c-a-b+1 mpc_sub(a, b, amb); // a+b-c+1 c-a-b+1 a=b RE(z) > 0 if not zeroorounder(ap) and not zeroorounder(am) and mpc_is0(amb) and mpf_is_gt(z.re, zero) then begin Hypergeometric_function_of_first_kind_aeb(a, b, c, z, eps, ans, n, f); memo1.Lines.Append('a = b, 1<z<=2 '); memo1.Lines.Append('変換 1-1/z Δ無し計算'); goto TN0E; end; end; end; //-------------------------------------------------------------------------- // c/a=2 c/b=2 二次変換 (a,b 非整数 |z|>1 は前で処理) // 入力値z= 1.4 近辺で最小値-4.7程度になり入力値が大きくなるとゼロに近づく if not mpc_is0(a) then mpc_div(c, a, ap); if not mpc_is0(b) then mpc_div(c, b, bp); if (checkbox2.Checked = true) and mpf_is_eq(ap.re, two) and mpf_is_eq(bp.re, two) and not mpc_is1(z) then begin mpc_copy(z, zback); // z=zback mpf_div(one, two, tmf); // 1/2 mpc_add_mpf(b, tmf, c); // c = b + 1/2 mpc_copy(a, aback); // a=>aback a/2 mpc_div_mpf(a, two, a); // a = a/2 mpc_sub(b, a, b); // b = b - a/2 mpc_mul_mpf(z, four, tmc); // 4z mpc_sub_mpf(tmc, four, tmc); // 4z - 4 mpc_div(z, tmc, tmc); // z / (4z-4) mpc_mul(tmc, z, z); // z^2 /(4z-4) memo1.Lines.Append('c/a=2 c/b=2'); memo1.Lines.Append(ZeroErase(string('z= z^2 /(4z-4)= ' + mpf_decimal(z.re, 40)))); mpf_set_dbl(tmf, 0.9); if mpf_is_gt(zback.re, tmf) then DS := 64 // z.re > 0.9 入力値z 1.16～7.6の時 1-1/z計算 else DS := 8; // z.re <=0.9 mpc_abs(z, absz); // |z| memo1.Lines.Append(ZeroErase(string('|z| =' + mpf_decimal(absz, 40)))); // 入力値z 1.16～7.6 // 0.9 < /z/ < 2.2 mpf_set_dbl(d1, 0.9); mpf_set_dbl(d2, 2.2); if mpf_is_gt(absz, d1) and mpf_is_lt(absz, d2) and mpf_is_ge(z.re, zero) and mpf_is0(z.im) and (DS<>8) then begin Hypergeometric_function_of_first_kind_aeb(a, b, c, z, eps, ans, n, f); memo1.Lines.Append('変換 (1-1/z 0.9<|z|<2.2 Δ無し)'); goto TN0E; end else goto ZG03; // Pfafの変換無し end; ZG02: //-------------------------------------------------------------------------- // Pfaffの変換公式 // a = b c / a = 2 c / b = 2 の時はpfaff使用しない(許可設定あり // Γ計算の都合上の±Δ計算の誤差が無ければ使用しても問題なし) // zの実数部がゼロの時で/z/=1ならPfafの変換使用 // |z|>=1 IM(z)=0の場合計算結果に虚数部がある場合虚数部±が逆になる場合があるので注意 // 1近傍は1-1/z の線形接続公式を使用 // 0,75<|z|<=1 又は 1<|z|<2 a,b,c 小数部0 虚数部0 // zの虚数部の値によりPfaffを使用しない方が良い場合があります|z|=1近辺 mpf_frac(a.re, ap.re); // re(a) 小数部 mpf_frac(b.re, bp.re); // re(b) 小数部 mpf_frac(c.re, cp.re); // re(c) 小数部 mpf_read_decimal(dnine, PAnsiChar(ansistring('0.75' + #00))); // 0,75<|z|<=1 if (mpf_is_gt(absz, dnine) and mpf_is_le(absz, one)) or // 又は 1<|z|<2 re(a),re(b),re(c)の小数部0 im(a),im(b),im(c)虚数部0 (mpf_is_gt(absz, one) and mpf_is_lt(absz, two) and mpf_is0(ap.re) and mpf_is0(bp.re) and mpf_is0(cp.re) and mpf_is0(a.im) and mpf_is0(b.im) and mpf_is0(c.im)) then begin // pfaff変換 有効 なら if (checkbox3.Checked = true) then begin if not mpc_is0(a) then mpc_div(c, a, ap); // c/a if not mpc_is0(b) then mpc_div(c, b, bp); // c/b // a,とbの実数部がc/2, zの虚数部が０でないなら if (checkbox7.Checked = false) and mpf_is_eq(ap.re, two) and mpf_is_eq(bp.re, two) and mpf_is0(z.im) then else begin mpc_sub(a, b, am); // a-b mpc_sub(c, b, cp); // c-b // |z| < 0.75 {z| < 2 z<>1 なら Pfafの変換 if mpf_is_gt(absz, dnine) and mpf_is_le(absz, two) and not mpc_is1(z) then begin // a=b c=a+1 c=b+1 でないなら 例 2,2,3 3,3,4でない // c-b=1がcとbの値によって桁落ちするので cp=1 を not cp.re<1<cp.re で代用 if not (mpc_is0(am) and not (mpc_is0(am) and (mpf_is_lt(cp.re, one) or mpf_is_gt(cp.re, one)))) then begin mpc_copy(z, zbackP); DS := DS or 4; // 4 pfaff 変換フラグセット mpc_sub(z, onec, tmc); // z - 1 mpc_div(z, tmc, z); // z = z/(z-1); mpc_abs(z, absz); // |z| mpc_sub(c, b, b); // b = c - b end; // a=b c=a+1 c=b+1 の時変換後aにΔ値加減算して近似計算 // 例 2,2,3 3,3,4 なら if (mpc_is0(am) and not (mpf_is_lt(cp.re, one) or mpf_is_gt(cp.re, one))) then begin mpc_copy(z, zbackP); DS := DS or 4 or 128; // 4 pfaff 変換フラグセット 128 は近似計算フラグ mpc_sub(z, onec, tmc); // z - 1 mpc_div(z, tmc, z); // z = z/(z-1); mpc_abs(z, absz); // |z| mpc_sub(c, b, b); // b = c - b XLE := true; end; end; end; end; end else begin // 1<|z|<2 の範囲は1-1/z計算使用 // z.re < 0 時は 1/z使用 mpc_add(a, b, ap); // a+b mpc_sub(ap, c, ap); // a+b-c mpc_sub(c, a, am); // c-a mpc_sub(am, b, am); // c-a-b if (checkbox2.Checked = true) // a+b-c c-a-b c and not zeroorounder(ap) and not zeroorounder(am) and not zeroorounder(c) then begin mpc_add(ap, onec, ap); // a+b-c+1 mpc_add(am, onec, am); // c-a-b+1 // a+b-c+1 c-a-b+1 |z| > 1 if not zeroorounder(ap) and not zeroorounder(am) and mpf_is_gt(absz, one) // /z/<2 RE(z) > 0 and mpf_is_lt(absz, two) and mpf_is_gt(z.re, zero) then begin Hypergeometric_function_of_first_kind_aeb(a, b, c, z, eps, ans, n, f); memo1.Lines.Append('変換 1-1/z 1<z<2 Δ無し計算'); goto TN0E; end; end; end; ZG03: //-------------------------------------------------------------------------- // |z| > 1 // Γ値の値が±∞になる場合、微小値Δを加減算して計算平均値を答えにしています。 // a=bの計算の時は、aのみ微小値Δを加減算 a=bとa<>bでは、平均値の計算が違います。 if mpf_is_gt(absz, one) then begin // |z|が1より大きかったら if (DS = 0) or (DS = 4) then XLE := false; mpc_sub(a, b, amb); mpc_sub(b, a, bma); mpc_sub(c, a, cma); mpc_sub(c, b, cmb); mpc_sub(a, b, aeb); df := 0; // amb= 0 cma整数>0 又は 近似計算フラグ時 a±Δ近似値計算フラグセット df = 16 // aebフラグ1 (a±Δによりa=bではなくなります) if (mpc_is0(amb) and not zeroorounder(cma)) or (DS and 16 = 16) or (DS and 128 = 128) then begin // DS=16 (a, b=a+1, b+1 b) df := 16; // aのみΔ加算減算計算 mpc_set1(aeb); // df=16時はa±Δで計算する為 aeb では無くなります end; // 1/zの線接続公式の時±無限大になるΓ検出 if (df = 16) then else begin if zeroorounder(amb) then df := df or 1; if zeroorounder(bma) then df := df or 2; if zeroorounder(cma) then df := df or 4; if zeroorounder(cmb) then df := df or 8; end; // 片側のΓ計算の分子と分母に一つの|∞|になるので±Δなし if (df = 6) or (df = 9) or (df = 15) then begin Hypergeometric_function_of_NO_delt(a, b, c, z, eps, ans, n, f); if f <> 128 then begin memo1.Lines.Append(' 1/z Δ無し計算 Γ/Γ'); XLE := false; goto TN0E; end; f := 0; end; // Γ計算分母が無限大になる場合±Δなし if (df and 12<> 0) and (df and 3 = 0) then begin Hypergeometric_function_of_NO_delt(a, b, c, z, eps, ans, n, f); if f <> 128 then begin memo1.Lines.Append(' 1/z Δ無し計算 Γ/∞'); XLE := false; goto TN0E; end; f := 0; end; // 両側のΓ計算の分母に|∞|片側の分子が|∞|になるので、a±Δの計算 if (df = 13) or (df = 14) then begin df := 16; mpc_set1(aeb); // df=16時はa±Δで計算する為 aeb では無くなります end; // |z| > 1 計算で Γ値が+∞-∞になるなら近似計算 Δα値を加算減算して計算し平均値を求めます if checkbox1.Checked = true then begin // a=bでΓ計算の分子が∞になり、分母に無限大がない場合計算結果が無限大になってしまうが、 // 実際には無限大ではないので、正解に近いΔ値を与えます。 if mpc_is0(amb) then mpf_read_decimal(d1, PAnsiChar(ansistring('2.8E-19' + #00))) else mpf_read_decimal(d1, PAnsiChar(ansistring('1E-45' + #00))); if DS and 128 = 128 then mpf_read_decimal(d1, PAnsiChar(ansistring('1.8E-19' + #00))); if mpc_is0(amb) then DS:= DS or 128; // a=b の時は近似計算表示フラグセット if df <> 0 then begin // 微小値のΔ加算減算が必用なら XLE := true; // df > 0 近似計算表示選択 XLE = true // a,b,cの何処にΔ補正するか選択 if (df = 1) or (df = 6) then df := 16; // a if (df = 2) or (df = 9) then df := 32; // b if (df = 4) or (df = 8) then df := 64; // c // ＋Δ 計算 mpc_copy(a, ap); mpc_copy(b, bp); mpc_copy(c, cp); if df = 16 then mpc_add_mpf(a, d1, ap); if df = 32 then mpc_add_mpf(b, d1, bp); if df = 64 then mpc_add_mpf(c, d1, cp); mpc_sub(ap, bp, amb); mpc_sub(bp, ap, bma); mpc_sub(cp, ap, cma); mpc_sub(cp, bp, cmb); Hypergeometric_function_of_first_kind(ap, bp, cp, amb, bma, cma, cmb, aeb, z, eps, tmc, n, f); // Δ加算計算 memo1.Lines.Append(ZeroErase(string('tmc ' + mpf_decimal(tmc.re, 40)))); // -Δ 計算 mpc_copy(a, am); mpc_copy(b, bm); mpc_copy(c, cm); if df = 16 then mpc_sub_mpf(a, d1, am); if df = 32 then mpc_sub_mpf(b, d1, bm); if df = 64 then mpc_sub_mpf(c, d1, cm); mpc_sub(am, bm, amb); mpc_sub(bm, am, bma); mpc_sub(cm, am, cma); mpc_sub(cm, bm, cmb); Hypergeometric_function_of_first_kind(am, bm, cm, amb, bma, cma, cmb, aeb, z, eps, ans, n, f); // Δ減算計算 memo1.Lines.Append(ZeroErase(string('ans ' + mpf_decimal(ans.re, 42)))); mpc_add(ans, tmc, ans); mpc_div_mpf(ans, two, ans); // 平均値 if df = 16 then memo1.Lines.Append(' 1/z a ±Δ 計算'); if df = 32 then memo1.Lines.Append(' 1/z b ±Δ 計算'); if df = 64 then memo1.Lines.Append(' 1/z c ±Δ 計算'); end else begin // 微小値のΔ加算減算が必用 ないなら Hypergeometric_function_of_NO_delt(a, b, c, z, eps, ans, n, f); // Hypergeometric_function_of_first_kind(a, b, c, amb, bma, cma, cmb, aeb, z, eps, ans, n, f); // Δ無し memo1.Lines.Append(' 1/z Δ=0 計算'); XLE := false; end; end else begin // 無限大回避のチェックを外すと使用されます Hypergeometric_function_of_NO_delt(a, b, c, z, eps, ans, n, f); memo1.Lines.Append(' 1/z Δ無し計算'); XLE := false; end; end else if not mpc_is1(z) then begin // z < 1 通常の超幾何関数計算 Hypergeometric_function(a, b, c, z, eps, ans, n, f); memo1.Lines.Append(' z < 1 ±Δ無し'); XLE := false; end; end; // TN=0 end TN0E: //---------------------------------------------------------------------------- if TN = 0 then begin // 通常計算 mpc_abs(zback, tmf); memo1.Lines.Append(ZeroErase(string(' |z| = ' + mpf_decimal(tmf, 50)))); mpc_abs(z, tmf); memo1.Lines.Append(ZeroErase(string('|nz| = ' + mpf_decimal(tmf, 50)))); end; //---------------------------------------------------------------------------- if TN = 1 then begin // arctan計算 mpc_mul(z, z, mz2); mpc_chs(mz2, mz2); if mpf_is_gt(absz, one) then begin // zが1より大きかったら mpc_sub(a, b, amb); mpc_sub(b, a, bma); mpc_sub(c, a, cma); mpc_sub(c, b, cmb); Hypergeometric_function_of_first_kind(a, b, c, amb, bma, cma, cmb, amb, mz2, eps, ans, n, f); end else // zが1より小さかったら Hypergeometric_function(a, b, c, mz2, eps, ans, n, f); end; //---------------------------------------------------------------------------- if TN = 2 then begin // arcsin計算 if not mpf_is1(az) then begin mpc_mul(z, z, mz2); Hypergeometric_function(a, b, c, mz2, eps, ans, n, f); end; end; //---------------------------------------------------------------------------- if BR then memo1.Lines.Append('計算が途中で打ち切られました、計算結果は正しくありません。'); if n >= NC then begin memo1.Lines.Append('計算が' + intTostr(NC) + 'Loopで収束しませんでした。'); memo1.Lines.Append('計算結果は正しくありません。'); end else begin if TN = 0 then begin // 通常計算 memo1.Lines.Append('ガウスの超幾何関数'); end; end; // ∞表示 if (f = 1) or (f = -1) then begin if mpf_is_ge(ans.re, zero) then memo1.Lines.Append('∞') else memo1.Lines.Append('-∞'); if mpf_is_gt(ans.im, zero) then memo1.Lines.Append('+∞i'); if mpf_is_lt(ans.im, zero) then memo1.Lines.Append('-∞i'); end; // Z>1 計算時Γ計算が無限大になる場合 f=128 if f = 128 then begin memo1.Lines.Append('Γ計算が+∞あるいは-∞になりました。'); memo1.Lines.Append('無限大回避にチェックを入れてください。'); goto EXT; end; //---------------------------------------------------------------------------- if TN = 0 then begin // 通常計算 // オイラー積分があったら if EUF then begin form1.memo1.Lines.Append('Euler integral = ' + string(' ' + mpf_decimal(euans.re, 3))); form1.memo1.Lines.Append(' ' + string(' ' + mpf_decimal(euans.im, 3)) + ' i'); end; // 絶対値が指定値より小さかったら0セット if (f = 16) or (f = -2) then begin memo1.Lines.Append('演算不可'); memo1.Lines.Append('無限大回避にチェックを入れてください。'); goto EXT; end; if XLE then memo1.Lines.Append(' 近似計算'); // z=1 計算 ∞時 if f = 8 then begin memo1.Lines.Append('z = 1'); memo1.Lines.Append(' Re(c - a - b) <= 0'); if mpf_is_ge(ans.re, zero) then memo1.Lines.Append(' ∞') else memo1.Lines.Append('-∞'); if mpf_is_gt(ans.im, zero) then memo1.Lines.Append('+∞i'); if mpf_is_lt(ans.im, zero) then memo1.Lines.Append('-∞i'); goto EXT; end; // pfaff変時出力設定 if DS and 4 = 4 then begin mpc_sub(onec, zbackp, tmc); // 1-z mpc_chs(a, ap); // -a mpc_powa(tmc, ap, tmc); // (1-z)^-a mpc_mul(ans, tmc, ans); // ans = ans * ((1-z)^-a) mpc_abs(zbackp, tmf); memo1.Lines.Append(' Pfaff'); // re(z) > 1 im(z)=0 の時 虚数部の符号が反転するので修正します。 if mpf_is_gt(zbackp.re, one) and mpf_is0(zbackp.im) then mpf_chs(ans.im, ans.im); end; // 二次 z^2/(z^2/(z4-4)) 変換時出力 必ずpfaffの後に実行 if (DS and 8 = 8) or (DS and 64 = 64) then begin mpc_set1(tmc); // 1 mpc_sub(tmc, zback, tmc); // 1-z mpc_chs(aback, aback); // -a mpc_div_mpf(aback, two, aback); // -a/2 mpc_powa(tmc, aback, tmc); // (1-z)^-a/2 mpc_mul(ans, tmc, ans); // ans = ans * ((1-z)^-a/2) memo1.Lines.Append(' F(a=b c=2b, z^2/(z4-4)'); end; // 近似計算時ゼロ設定 if XLE then begin if (DS and 32 = 32) or (DS and 128 = 128) then zerounde1em35(ans) // 1E-35以下ゼロ設定 else zerounde1em42(ans); // 1E-42以下ゼロ設定 end else zerounde1em52(ans); // 1E-52以下ゼロ設定 epslimit(ans); // 収束判定値による表示制限 // 計算結果表示 近似計算時は40桁 通常は50桁表示 a=b でpfaff時は35桁表示 if (f = 0) or (f = 4) or (f = 16) then begin memo1.Lines.Append('ans'); mpf_abs(ans.re, az); if mpf_is_gt(az, infs) then begin if mpf_is_gt(ans.re, zero) then memo1.Lines.Append(' ∞') else memo1.Lines.Append(' -∞'); end else begin if XLE and ((DS and 32 = 32) or (DS and 128 = 128)) then begin memo1.Lines.Append(' 精度が低くなっています。'); end; if XLE then begin if (DS and 32 = 32) or (DS and 128 = 128) then memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.re, 33)))) else memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.re, 40)))); end else memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.re, 50)))); end; mpf_abs(ans.im, az); if mpf_is_gt(az, infs) then begin if mpf_is_gt(ans.im, zero) then memo1.Lines.Append(' +∞i') else memo1.Lines.Append(' -∞i'); end else begin if XLE then begin if (DS and 32 = 32) or (DS and 128 = 128) then memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.im, 33))) + ' i') else memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.im, 40))) + ' i'); end else memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(ans.im, 50))) + ' i'); end; end; end; //---------------------------------------------------------------------------- if TN = 1 then begin // arctan計算 mpc_mul(ans, z, tans); // z * if DS = 1 then memo1.Lines.Append('arctan(z) rad 超幾何関数計算'); if DS = 2 then memo1.Lines.Append('arcsin(z) rad 超幾何関数計算'); memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(tans.re, 50)))); mpc_arctan(z, dat); memo1.Lines.Append('rad 組込み関数計算'); memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(dat.re, 50)))); end; //---------------------------------------------------------------------------- if TN = 2 then begin // arcsin計算 if not mpf_is1(az) then mpc_mul(ans, z, tans) // z * else begin mpf_set_pi(tans.re); mpf_div(tans.re, two, tans.re); mpc_mul(tans, z, tans) // z * end; if DS = 1 then memo1.Lines.Append('arctan(z) rad 超幾何関数計算'); if DS = 2 then memo1.Lines.Append('arcsin(z) rad 超幾何関数計算'); memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(tans.re, 50)))); mpc_arcsin(z, dat); memo1.Lines.Append('rad 組込み関数計算'); memo1.Lines.Append(ZeroErase(string(' ' + mpf_decimal(dat.re, 50)))); end; EXT: button1.Enabled := True; TN := 0; // artn arsin 解除 DS := 0; // artn arsin pfaff 二次変換 近似計算 解除 mpc_clear5(a, b, c, z, ans); mpf_clear3(eps, az, tmf); mpc_clear5(mz2, tans, dat, tmc, onec); mpc_clear5(amb, bma, cma, cmb, aeb); mpf_clear4(d1, d2, dnine, absz); mpc_clear4(ap, bp, cp, zbackp); mpc_clear5(am, bm, cm, zback, aback); mpc_clear(euans); end; //------------------------------------------------------------------------------ // 初期設定 procedure TForm1.FormCreate(Sender: TObject); var i : integer; N, D, tmp : mp_float; ctmp : mp_complex; begin mpf_set_default_prec(1000); precisionedit.Text := intTostr(mpf_get_default_prec); memo1.Clear; Label1.Caption := ' 設定値によって、中々収束しません。' + #13#10 + 'その時は計算打切りボタンで停止して' + #13#10 + '下さい。'; Label2.Caption := ' RE(z)≒0.5 |z|≒1は中々収束しません、収束判定値を' + #13#10 + '大きくすれば精度は悪くなり表示桁数は少なくなりますが、' + #13#10 + '計算は早く終了します。 例 1e-10, 1e-20'; BR := false; TN := 0; // 計算実行 0 通常 1 artan 2 arsin DS := 0; // 計算表示 0 通常 1 artan 2 arsin 1 4 pfaff 8,64 二次変換 32,128 近似計算 setlength(BM, NB + 1); // Γ関数用ベルヌーイ数配列 for i := 0 to NB do mpf_init(BM[i]); mpf_init3(N, D, tmp); mpf_init5(zero, one, two, four, three); mpf_init3(pai, maxmpf, infs); mpc_init(log_2pis2); mpc_init(ctmp); // 0,1,2,3,4, π 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_read_decimal(maxmpf, PAnsiChar(ansistring('1.0e+100000000' + #00))); // ∞判別値 mpf_read_decimal(infs, PAnsiChar(ansistring('1.0e+1000000' + #00))); // ln(2π)/2 複素数 mpc_set0(ctmp); // 0 + 0i mpf_mul(pai, two, ctmp.re); // 2π+0i mpc_ln(ctmp, ctmp); // ln(2π+0i) mpc_div_mpf(ctmp, two, log_2pis2); // ln(2π+0i)/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]); // BM[] Γ関数用ベルヌーイ数配列 end; comment; mpf_clear3(N, D, tmp); mpc_clear(ctmp); end; // close時 多倍長配列解放 procedure TForm1.FormClose(Sender: TObject; var Action: TCloseAction); var i : integer; begin for i := 0 to NB do mpf_clear(BM[i]); mpf_clear5(zero, one, two, four, three); mpf_clear3(pai, maxmpf, infs); mpc_clear(log_2pis2); end; procedure TForm1.FormCloseQuery(Sender: TObject; var CanClose: Boolean); begin // 演算中 close禁止 if not button1.Enabled then CanClose := false; end; procedure TForm1.Button2Click(Sender: TObject); begin comment; end; end.
Gauss_hypergeometric_function_A.zip
三角関数、逆三角関数、その他関数 に戻る