多倍長の計算は、有効桁数74桁程度なので、ゼロがE-74になる誤差が出ています。表示の桁数は少ないですが、計算は有効桁数50桁以上の精度を求めている為、計算に可成りの時間を要します。
カールソンの楕円積分(複素数)Ｎｏ３
ＲＣ，ＲＦ，ＲＪの多倍長の計算
多倍長の計算は、Ｆｕｎｃｔｉｏｎの使用が出来ないので、計算が長くなります。
多倍長の計算は、有効桁数74桁程度なので、ゼロがE-74になる誤差が出ています。
表示の桁数は少ないですが、計算は有効桁数50桁以上の精度を求めている為、計算に可成りの時間を要します。
多倍長の組み込みは MPArithからmparith_2018-11-27.zipをダウンロードして下さい。
Zipファイルを解凍して、解凍したフォルダーにパス(ツール->オプション->Delphiオプション->ライブラリ->ライブラリパス)を追加し、Uses に mp_types, mp_real, mp_cmplx, mp_baseを追加すれば使用可能になります。
プログラム ＲＪルーチンの中に分母がゼロになった場合の回避ルーチンがありますが、有効演算桁数が多い為、プログラムミスがない限り回避ルーチンが動作することはありません。
// 1/31/2021 // プログラムは https://arxiv.org/abs/math/9409227 のレポートを元に作成しています 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; type TForm1 = class(TForm) Memo1: TMemo; RJtest1: TBitBtn; Rjdtxr: TEdit; Rjdtyr: TEdit; Rjdtzr: TEdit; Rjdtpr: TEdit; Rjansr: TEdit; Rjansi: TEdit; Rjclc0: TEdit; Rjclc1: TEdit; RjRadio: TRadioGroup; Rjdtxi: TEdit; Rjdtyi: TEdit; Rjdtzi: TEdit; Rjdtpi: TEdit; Rjtest2: TBitBtn; Rftest1: TBitBtn; Rftest2: TBitBtn; Rctest1: TBitBtn; Rctest2: TBitBtn; Rfdtxr: TEdit; Rfdtxi: TEdit; Rfdtyr: TEdit; Rfdtyi: TEdit; Rfdtzr: TEdit; Rfdtzi: TEdit; Rfansr: TEdit; Rfansi: TEdit; Rfclc0: TEdit; Rfclc1: TEdit; RfRadio: TRadioGroup; Rcdtxr: TEdit; Rcdtxi: TEdit; Rcdtyr: TEdit; Rcdtyi: TEdit; Rcclc0: TEdit; Rcclc1: TEdit; Rcansi: TEdit; Rcansr: TEdit; RcRadio: TRadioGroup; procedure RJtest1Click(Sender: TObject); procedure Rjtest2Click(Sender: TObject); procedure Rftest1Click(Sender: TObject); procedure Rftest2Click(Sender: TObject); procedure Rctest1Click(Sender: TObject); procedure Rctest2Click(Sender: TObject); private { Private 宣言 } public { Public 宣言 } end; var Form1: TForm1; implementation {$R *.dfm} uses system.Math, mp_real, mp_cmplx, mp_types, mp_base; procedure Vc(a, b: double; var ans :mp_complex); begin mpc_set_dbl(ans, a, b); end; procedure mpf_max(x, y: mp_float; var ans: mp_float); begin if mpf_is_gt(x, y) then mpf_copy(x, ans) else mpf_copy(y, ans); end; // RC procedure RCd(x, y: mp_complex; var ans: mp_complex); label Ext; const r = 1E-200; c1d = 3 / 10; c2d = 1 / 7; c3d = 3 / 8; c4d = 9 / 22; c5d = 159 / 208; c6d = 9 / 8; var xt, yt : mp_complex; lamda, Am : mp_complex; A0, s, yb : mp_complex; one, p4, w : mp_complex; c1, c2, c3 : mp_complex; c4, c5, c6 : mp_complex; d2 : mp_complex; sqrtx, sqrty : mp_complex; tmp, tmpa : mp_complex; Q, tmpf, tmpaf : mp_float; m : integer; begin mpc_init2(xt, yt); mpc_init2(lamda, Am); mpc_init3(A0, s, yb); mpc_init3(one, p4, w); mpc_init3(c1, c2, c3); mpc_init3(c4, c5, c6); mpc_init2(sqrtx, sqrty); mpc_init2(tmp, tmpa); mpc_init(d2); mpf_init3(Q, tmpf, tmpaf); mpc_abs(y, tmpf); if mpf_is0(tmpf) then begin // yがゼロだったら終了 mpc_set_dbl(ans, 0, 0); m := 0; goto Ext; end; mpc_set_dbl(one, 1, 0); mpc_set_dbl(c1, c1d, 0); mpc_set_dbl(c2, c2d, 0); mpc_set_dbl(c3, c3d, 0); mpc_set_dbl(c4, c4d, 0); mpc_set_dbl(c5, c5d, 0); mpc_set_dbl(c6, c6d, 0); mpc_set_dbl(d2, 2, 0); if s_mpf_is_gt0(y.re) or not s_mpf_is0(y.im) then begin mpc_copy(x, xt); mpc_copy(y, yt); mpc_set1(w); end else begin mpc_sub(x, y, xt); mpc_chs(y, yt); mpc_div(x, xt, w); mpc_sqrt(w, w); end; mpc_copy(yt, yb); // yb = yt mpc_add(yt, yt, tmp); // yt + yt mpc_add(tmp, xt, tmp); // 2yt + xt mpc_div_dbl(tmp, 3, A0); // (2yt + xt) / 3 mpc_copy(A0, Am); // Am = A0 mpf_set_dbl(tmpf, 3 * r); // 3 * r mpf_set_dbl(tmpaf, 1 / 8); // 1 / 8 mpf_expt(tmpf, tmpaf, tmpf); // power(3 * r, 1 / 8) mpf_inv(tmpf, tmpf); // 1 / power(3 * r, 1 / 8) mpc_sub(A0, xt, tmpa); // A0 - xt mpc_abs(tmpa, tmpaf); // abs(A0 - xt) mpf_mul(tmpf, tmpaf, Q); // Q = 1 / power(3 * r, 1 / 8) * abs(A0 - xt) m := 0; repeat mpc_sqrt(xt, sqrtx); // √x mpc_sqrt(yt, sqrty); // √y mpc_mul(sqrtx, sqrty, tmp); // sqrtx * sqrty mpc_mul(d2, tmp, tmpa); // 2 * (sqrtx * sqrty) mpc_add(tmpa, yt, lamda); // lamda = 2 * (sqrtx * sqrty) + y mpc_add(xt, lamda, tmp); // x + lamda mpc_div_dbl(tmp, 4.0, xt); // x = (x + lamda) / 4 mpc_add(yt, lamda, tmp); // y + lamda mpc_div_dbl(tmp, 4.0, yt); // y = (y + lamda) / 4 mpc_add(Am, lamda, tmp); // Am + lamda mpc_div_dbl(tmp, 4.0, Am); // Am = (Am + lamda) / 4 inc(m); mpf_set_dbl(tmpf, 4); // 4 mpf_set_dbl(tmpaf, -m); // -m mpf_expt(tmpf, tmpaf, tmpf); // 4^-m mpf_mul(tmpf, Q, tmpaf); // 4^-m * Q mpc_abs(Am, tmpf); // abs(am) until mpf_is_lt(tmpaf,tmpf); // 4^-m * Q < cabs(Am) mpf_set_int(tmpf, 4); mpf_set_int(tmpaf, m); mpf_expt(tmpf, tmpaf, tmpf); // 4^m mpf_set_int(tmpaf, 0); mpc_set_mpf(p4, tmpf, tmpaf); // p4 = 4^m + 0i mpc_sub(yb, A0, tmp); // yb - A0 mpc_div(tmp, p4, tmpa); // (yb - A0) / p4 mpc_div(tmpa, Am, s); // s = (yb - A0) / p4 / Am mpc_mul(s, c6, tmpa); mpc_add(tmpa, c5, tmp); mpc_mul(s, tmp, tmpa); mpc_add(tmpa, c4, tmp); mpc_mul(s, tmp, tmpa); mpc_add(tmpa, c3, tmp); mpc_mul(s, tmp, tmpa); mpc_add(tmpa, c2, tmp); mpc_mul(s, tmp, tmpa); mpc_add(tmpa, c1, tmp); mpc_mul(s, tmp, tmpa); mpc_mul(s, tmpa, tmp); mpc_add(one, tmp, tmpa); mpc_mul(w, tmpa, tmp); // w * ??? mpc_sqrt(Am, tmpa); mpc_div(tmp, tmpa, ans); // ans = w * ???? / sqrt(Am) Ext: Form1.Memo1.Lines.Append('RC Loop数 = ' + intTostr(m)); mpc_clear2(xt, yt); mpc_clear2(lamda, Am); mpc_clear3(A0, s, yb); mpc_clear3(one, p4, w); mpc_clear3(c1, c2, c3); mpc_clear3(c4, c5, c6); mpc_clear2(sqrtx, sqrty); mpc_clear2(tmp, tmpa); mpc_clear(d2); mpf_clear3(Q, tmpf, tmpaf); end; // RF procedure RFd(x, y, z: mp_complex; var ans : mp_complex); label Ext; const r = 1E-200; c1d = 1 / 10; c2d = 1 / 14; c3d = 1 / 24; c4d = 3 / 44; var one : mp_complex; lamda : mp_complex; xm, ym, zm : mp_complex; A0, Am : mp_complex; sqrtx, sqrty, sqrtz : mp_complex; m : integer; X0, Y0, Z0 : mp_complex; E2, E3 : mp_complex; c1, c2, c3, c4 : mp_complex; d3, d4, p4, tmp : mp_complex; tmpa, tmpb, tmpc : mp_complex; tmpf, tmpaf, tmpbf : mp_float; Q, M4 : mp_float; begin mpc_init2(one, lamda); mpc_init3(xm, ym, zm); mpc_init2(A0, Am); mpc_init3(sqrtx, sqrty, sqrtz); mpc_init3(X0, Y0, Z0); mpc_init2(E2, E3); mpc_init4(c1, c2, c3, c4); mpc_init4(d3, d4, p4, tmp); mpc_init3(tmpa, tmpb, tmpc); mpf_init3(tmpf, tmpaf, tmpbf); mpf_init2(Q, M4); m := 0; mpc_abs(x, tmpf); if mpf_is0(tmpf) then inc(m); mpc_abs(y, tmpf); if mpf_is0(tmpf) then inc(m); mpc_abs(z, tmpf); if mpf_is0(tmpf) then inc(m); if m > 1 then begin // ゼロ二個有ったら終了 mpc_set_dbl(ans, 0, 0); m := 0; goto Ext; end; mpc_set1(one); mpc_set_dbl(c1, c1d, 0); mpc_set_dbl(c2, c2d, 0); mpc_set_dbl(c3, c3d, 0); mpc_set_dbl(c4, c4d, 0); mpc_set_dbl(d3, 3, 0); mpc_set_dbl(d4, 4, 0); mpf_set_dbl(M4, 4); mpc_copy(x, xm); mpc_copy(y, ym); mpc_copy(z, zm); mpc_add(ym, zm, tmp); // x + y mpc_add(tmp, xm, tmpa); // x + y + z mpc_div(tmpa, d3, A0); // (x + y + z) / 3 mpc_copy(A0, Am); // Am = A0 mpc_sub(A0, xm, tmp); // A0 - x mpc_sub(A0, ym, tmpa); // A0 - Y mpc_sub(A0, zm, tmpb); // A0 - z mpc_abs(tmp, tmpf); // |A0 - x| mpc_abs(tmpa, tmpaf); // |A0 - y| mpc_abs(tmpb, tmpbf); // |A0 - z| mpf_max(tmpf, tmpaf, tmpf); // max(|A0-x|,|A0-y|) mpf_max(tmpf, tmpbf, tmpf); // max(|A0-x|,|A0-y|,|A0-z|) mpf_set_dbl(tmpaf, 3 * r); // (3 * r) mpf_set_dbl(tmpbf, -1 / 6); // -1/6 mpf_expt(tmpaf, tmpbf, tmpaf); // (3 * r)^-1/6 or (3 / r)^1/6 mpf_mul(tmpf, tmpaf, Q); // Q = (3 * r)^-1/6 * max((A0-x),(A0-y),(A0-z)) m := 0; repeat mpc_sqrt(xm, sqrtx); // √x mpc_sqrt(ym, sqrty); // √y mpc_sqrt(zm, sqrtz); // √z mpc_mul(sqrty, sqrtz, tmp); // √y * √z mpc_add(sqrty, sqrtz, tmpa); // √y + √z mpc_mul(tmpa, sqrtx, tmpc); // (√y + √z) * √x mpc_add(tmpc, tmp, lamda); // lamda = √x√y + √x√z + √y√z mpc_add(xm, lamda, tmp); // xm + lamda mpc_div(tmp, d4, xm); // zm = (xm + lamda) / 4 mpc_add(ym, lamda, tmp); // ym + lamda mpc_div(tmp, d4, ym); // ym = (ym + lamda) / 4 mpc_add(zm, lamda, tmp); // zm + lamda mpc_div(tmp, d4, zm); // zm = (zm + lamda) / 4 mpc_add(Am, lamda, tmp); // Am + lamda mpc_div(tmp, d4, Am); // Am = (Am + lamda) / 4 m := m + 1; // inc(m) mpf_set_int(tmpaf, -m); // -m mpf_expt(M4, tmpaf, tmpbf); // 4^-m mpf_mul(tmpbf, Q, tmpf); // 4^-m * Q mpc_abs(Am, tmpaf); // abs(Am) until mpf_is_lt(tmpf, tmpaf); // 4^-m * Q < abs(Am) mpf_set_int(tmpaf, m); // m mpf_expt(M4, tmpaf, tmpf); // 4^m mpf_set0(tmpbf); // 0 mpc_set_mpf(P4, tmpf, tmpbf); // P4 = (4^m, 0i) mpc_sub(A0, x, tmp); // A0 - x mpc_div(tmp, p4, tmpa); // (A0 - x) / 4^m mpc_div(tmpa, Am, X0); // X0 = (A0 - x) / 4^m / Am mpc_sub(A0, y, tmp); // A0 - y mpc_div(tmp, p4, tmpa); // (A0 - y) / 4^m mpc_div(tmpa, Am, Y0); // Y0 = (A0 - y) / 4^m / Am mpc_chs(X0, tmp); // -X0 mpc_sub(tmp, Y0, Z0); // Z0 = -X0 - Y0 mpc_mul(Z0, Z0, tmpa); // Z0^2 mpc_mul(X0, Y0, tmpb); // X0 * Y0 mpc_sub(tmpb, tmpa, E2); // E2 = X0 * Y0 - Z0^2 mpc_mul(X0, Y0, tmpa); // X0 * Y0 mpc_mul(tmpa, Z0, E3); // E3 = X0 * Y0 * Z0 mpc_mul(c1, E2, tmpa); // c1 * E2 mpc_sub(one, tmpa, tmp); // 1 - C1 * E2 mpc_mul(c2, E3, tmpa); // c2 * E3 mpc_add(tmp, tmpa, tmp); // 1 - C1 * E2 + c2 * E3 mpc_mul(c3, E2, tmpa); // c3 * E2 mpc_mul(tmpa, E2, tmpb); // c3 * E2^2 mpc_add(tmp, tmpb, tmp); // 1 - C1 * E2 + c2 * E3 + c3 * E2^2 mpc_mul(c4, E2, tmpa); // c4 * E2 mpc_mul(tmpa, E3, tmpb); // c4 * E2 * E3 mpc_sub(tmp, tmpb, tmp); // 1 - C1 * E2 + c2 * E3 + c3 * E2^2 - c4 * E2 * E3 mpc_sqrt(Am, tmpc); // √Am mpc_div(one, tmpc, tmpb); // 1 / √Am mpc_mul(tmpb, tmp, ans); // ans = 1 / √Am * (1 - C1 * E2 + c2 * E3 + c3 * E2^2 - c4 * E2 * E3 EXT: Form1.Memo1.Lines.Append('RF Loop数 = ' + intTostr(m)); mpc_clear2(one, lamda); mpc_clear3(xm, ym, zm); mpc_clear2(A0, Am); mpc_clear3(sqrtx, sqrty, sqrtz); mpc_clear3(X0, Y0, Z0); mpc_clear2(E2, E3); mpc_clear4(c1, c2, c3, c4); mpc_clear4(d3, d4, p4, tmp); mpc_clear3(tmpa, tmpb, tmpc); mpf_clear3(tmpf, tmpaf, tmpbf); mpf_clear2(Q, M4); end; // 複素数３個の最小値 procedure minc3d(a, b, c: mp_complex; var ans: mp_complex); var tmp : mp_complex; ma, mb, mc : mp_float; begin mpc_init(tmp); mpf_init3(ma, mb, mc); mpc_abs(a, ma); mpc_abs(b, mb); mpc_abs(c, mc); if mpf_is_lt(ma, mb) then mpc_copy(a, tmp) else mpc_copy(b, tmp); mpc_abs(tmp, ma); if mpf_is_lt(mc, ma) then mpc_copy(c, ans) else mpc_copy(tmp, ans); mpc_clear(tmp); mpf_clear3(ma, mb, mc); end; // 複素数３個の最大値 procedure maxc3d(a, b, c: mp_complex; var ans: mp_complex); var tmp : mp_complex; ma, mb, mc : mp_float; begin mpc_init(tmp); mpf_init3(ma, mb, mc); mpc_abs(a, ma); mpc_abs(b, mb); mpc_abs(c, mc); if mpf_is_lt(ma, mb) then mpc_copy(b, tmp) else mpc_copy(a, tmp); mpc_abs(tmp, ma); if mpf_is_lt(mc, ma) then mpc_copy(tmp, ans) else mpc_copy(c, ans); mpc_clear(tmp); mpf_clear3(ma, mb, mc); end; // 複素数4個の最大値 procedure maxc4d(a, b, c, d: mp_complex; var ans: mp_complex); var tmp : mp_complex; ma, mb, mc, md : mp_float; begin mpc_init(tmp); mpf_init4(ma, mb, mc, md); mpc_abs(a, ma); mpc_abs(b, mb); mpc_abs(c, mc); mpc_abs(d, md); if mpf_is_lt(mb, ma) then mpc_copy(a, tmp) else mpc_copy(b, tmp); mpc_abs(tmp, ma); if mpf_is_lt(ma, mc) then mpc_copy(c, tmp); mpc_abs(tmp, ma); if mpf_is_lt(ma, md) then mpc_copy(d, ans) else mpc_copy(tmp, ans); mpc_clear(tmp); mpf_clear4(ma, mb, mc, md); end; // RJ procedure RJd(x, y, z, p: mp_complex; var ans: mp_complex); label Ext; const r = 1E-300; var one : mp_complex; rj, xt, yt, zt, pt : mp_complex; sum, a, b, rho, A0 : mp_complex; tau, rcx, lamda, Am : mp_complex; sqrtx, sqrty, sqrtz, sqrtp : mp_complex; dm, em, delta : mp_complex; c1, c2, c3 : mp_complex; c4, c5, c6 : mp_complex; X0, Y0, Z0, P0 : mp_complex; E2, E3, E4, E5 : mp_complex; d0, d2, d3, d4, d5 : mp_complex; d6, p4, powa, def : mp_complex; tmpa, tmpb, tmpc, tmpd : mp_complex; tmpaf, tmpbf, tmpcf : mp_float; Q, pm : mp_float; test, m : integer; ERRF : integer; begin mpc_init(one); mpc_init5(rj, xt, yt, zt, pt); mpc_init5(sum, a, b, rho, A0); mpc_init4(tau, rcx, lamda, Am); mpc_init4(sqrtx, sqrty, sqrtz, sqrtp); mpc_init3(dm, em, delta); mpc_init3(c1, c2, c3); mpc_init3(c4, c5, c6); mpc_init4(X0, Y0, Z0, P0); mpc_init4(E2, E3, E4, E5); mpc_init5(d0, d2, d3, d4, d5); mpc_init4(d6, p4, powa, def); mpc_init4(tmpa, tmpb, tmpc, tmpd); mpf_init3(tmpaf, tmpbf, tmpcf); mpf_init2(Q, pm); m := 0; mpc_abs(x, tmpaf); if mpf_is0(tmpaf) then inc(m); mpc_abs(y, tmpaf); if mpf_is0(tmpaf) then inc(m); mpc_abs(z, tmpaf); if mpf_is0(tmpaf) then inc(m); mpc_abs(p, tmpaf); if mpf_is0(tmpaf) then inc(m); if m > 1 then begin // ゼロが二個有ったら終了 mpc_set_dbl(ans, 0, 0); m := 0; goto Ext; end; mpc_set_dbl(c1, 3 / 14, 0); mpc_set_dbl(c2, 1 / 6, 0); mpc_set_dbl(c3, 9 / 88, 0); mpc_set_dbl(c4, 3 / 22, 0); mpc_set_dbl(c5, 9 / 52, 0); mpc_set_dbl(c6, 3 / 26, 0); mpc_set_dbl(d0, 0, 0); mpc_set_dbl(d2, 2, 0); mpc_set_dbl(d3, 3, 0); mpc_set_dbl(d4, 4, 0); mpc_set_dbl(d5, 5, 0); mpc_set_dbl(d6, 6, 0); mpc_set_dbl(sum, 0, 0); mpc_set_dbl(one, 1, 0); ERRF := 0; if s_mpf_is_le0(p.re) and s_mpf_is_gt0(y.re) then test := -1 else test := 1; if test > 0 then begin mpc_copy(x, xt); mpc_copy(y, yt); mpc_copy(z, zt); mpc_copy(p, pt); mpc_set_dbl(a, 0, 0); mpc_set_dbl(b, 0, 0); mpc_set_dbl(rcx, 0, 0); end else begin minc3d(x, y, z, xt); // xt = min(x, y, z) maxc3d(x, y, z, zt); // zt = max(x, y, z) mpc_add(xt, zt, tmpa); // xt + zt mpc_sub(z, tmpa, tmpb); // z - (xt + zt) mpc_add(y, tmpb, tmpc); // z - (xt + zt) + y mpc_add(x, tmpc, yt); // yt := z - (xt + zt) + y + x mpc_sub(yt, p, tmpa); // yt - p mpc_abs(tmpa, tmpaf); // |yt - p| if not s_mpf_is0(tmpaf) then begin // 分母ゼロ回避 {プログラムにバグがない限り必要ありませんデバッグ用です} mpc_inv(tmpa, a); // a = 1 / (yt - p) end else begin mpc_set_dbl(a, 1E200, 0); ERRF := ERRF or 00000010; end; mpc_sub(yt, xt, tmpa); // yt - xt mpc_sub(zt, yt, tmpb); // zt - yt mpc_mul(a, tmpa, tmpc); // a * (yt - xt) mpc_mul(tmpc, tmpb, b); // b = a * (yt - xt) * (zt - yt) mpc_add(yt, b, pt); // pt = yt + b mpc_abs(yt, tmpaf); // |yt| if not s_mpf_is0(tmpaf) then begin // 分母ゼロ回避 {プログラムにバグがない限り必要ありませんデバッグ用です} mpc_mul(xt, zt, tmpa); // xt * zt; mpc_div(tmpa, yt, rho); // rho = xt * zt / yt mpc_mul(p, pt, tmpb); // p * pt mpc_div(tmpb, yt, tau); // tau := p * pt / yt end else begin mpc_set_dbl(rho, 1E200, 0); mpc_set_dbl(tau, 1E200, 0); ERRF := ERRF or 00000100; end; rcd(rho, tau, rcx); // rcx := rcd(rho, tau) end; mpc_add(pt, pt, tmpa); // 2pt mpc_add(tmpa, zt, tmpb); // zt + 2pt mpc_add(tmpb, yt, tmpa); // yt + zt + 2pt mpc_add(tmpa, xt, tmpb); // xt + yt + zt + 2pt mpc_div(tmpb, d5, A0); // A0 = (xt + yt + zt + 2pt) / 5 mpc_copy(A0, Am); // Am = A0 mpc_sub(pt, xt, tmpa); // pt - xt mpc_sub(pt, yt, tmpb); // pt - yt mpc_sub(pt, zt, tmpc); // pt - zt mpc_mul(tmpa, tmpb, tmpd); // ( pt - xt) * (pt - yt) mpc_mul(tmpd, tmpc, delta); // deta = ( pt - xt) * (pt - yt) * pt - zt mpf_set_dbl(tmpaf, r / 4); // r / 4 mpf_set_dbl(tmpbf, -1 / 6); // -1 / 6 mpf_expt(tmpaf, tmpbf, tmpcf); // (r / 4)^(-1 / 6) mpc_sub(A0, xt, tmpa); // A0 - xt mpc_sub(A0, yt, tmpb); // A0 - yt mpc_sub(A0, zt, tmpc); // A0 - yz mpc_sub(A0, pt, tmpd); // A0 - pt maxc4d(tmpa, tmpb, tmpc, tmpd, tmpa); // Max(?????) mpc_abs(tmpa, tmpaf); // |Max(?????)| mpf_mul(tmpcf, tmpaf, Q); // Q = (r / 4)^(-1 / 6) * |Max(?????)| m := 0; repeat mpc_sqrt(xt, sqrtx); // √x mpc_sqrt(yt, sqrty); // √y mpc_sqrt(zt, sqrtz); // √z mpc_sqrt(pt, sqrtp); // √p mpc_add(sqrty, sqrtz, tmpa); // √y + √z mpc_mul(tmpa, sqrtx, tmpb); // (√y + √z) * √x mpc_mul(sqrty, sqrtz, tmpc); // √y * √z mpc_add(tmpb, tmpc, lamda); // lamda = (√y + √z) * √x + √y * √z mpc_add(sqrtp, sqrtx, tmpa); // √p - √x mpc_add(sqrtp, sqrty, tmpb); // √p - √y mpc_add(sqrtp, sqrtz, tmpc); // √p - √z mpc_mul(tmpa, tmpb, tmpd); // (√p - √x)(√p - √y) mpc_mul(tmpd, tmpc, dm); // dm = (√p - √x)(√p - √y)(√p - √z) mpc_set_dbl(tmpa, -3 * m, 0); // (-3m, 0i) mpc_pow(d4, tmpa, p4); // 4^(-3m, 0i) mpc_abs(dm , tmpaf); if not s_mpf_is0(tmpaf) then begin // ゼロで除算回避 {プログラムにバグがない限り必要ありませんデバッグ用です} mpc_mul(p4, delta, tmpa); // 4^(-3m, 0i) * δ mpc_div(tmpa, dm, tmpb); // 4^(-3m, 0i) * δ / dm mpc_div(tmpb, dm, em); // em = 4^(-3m, 0i) * δ / dm / dm end else begin mpc_set_dbl(em, 2E200, 0); ERRF := ERRF or 00001000; end; mpc_add(one, em, tmpa); // 1 + em rcd(one, tmpa, tmpb); // rc(1, 1+ em) mpc_set_dbl(tmpc, -m, 0); // -m mpc_pow(d4, tmpc, p4); // 4^-m mpc_mul(p4, tmpb, tmpd); // 4^-m * rc(1 + em) if not s_mpf_is0(tmpaf) then begin // ゼロで除算回避 {プログラムにバグがない限り必要ありませんデバッグ用です} mpc_div(tmpd, dm, def); // def = 4^-m / dm * rc(1 + em) end else begin mpc_set_dbl(def, 1E200, 0); ERRF := ERRF or 00010000; end; mpc_add(sum, def, sum); // sum = sum + 4^-m / dm * rc(1 + em) mpc_add(Am, lamda, tmpd); // Am + lamda mpc_div(tmpd, d4, Am); // Am = (Am + lamda) / 4 mpc_add(xt, lamda, tmpd); // xt + lamda mpc_div(tmpd, d4, xt); // xt = (xt + lamda) / 4 mpc_add(yt, lamda, tmpd); // yt + lamda mpc_div(tmpd, d4, yt); // yt = (yt + lamda) / 4 mpc_add(zt, lamda, tmpd); // zt + lamda mpc_div(tmpd, d4, zt); // zt = (zt + lamda) / 4 mpc_add(pt, lamda, tmpd); // pt + lamda mpc_div(tmpd, d4, pt); // pt = (pt + lamda) / 4 m := m + 1; mpc_set_dbl(tmpa, -m, 0); // (-m, 0i) mpc_pow(d4, tmpa, p4); // 4^(-3m, 0i) mpc_abs(p4, tmpbf); // |4^(-3m, 0i)| mpf_mul(tmpbf, Q, tmpaf); // |4^(-3m, 0i)| * Q mpc_abs(Am, tmpcf); // |Am| until mpf_is_lt(tmpaf, tmpcf); // |4^(-3m, 0i)| * Q < |Am| mpc_set_dbl(tmpa, m, 0); // (m, 0i) mpc_pow(d4, tmpa, p4); // 4^(3m, 0i) mpc_mul(p4, Am, tmpd); // 4^(3m, 0i) * Am mpc_sub(A0, x, tmpb); // A0 - x mpc_div(tmpb, tmpd, X0); // X0 = (A0 - x) / 4^3m / Am mpc_sub(A0, y, tmpb); // A0 - y mpc_div(tmpb, tmpd, Y0); // Y0 = (A0 - y) / 4^3m / Am mpc_sub(A0, z, tmpb); // A0 - z mpc_div(tmpb, tmpd, Z0); // Z0 = (A0 - z) / 4^3m / Am mpc_chs(X0, tmpa); // -X0 mpc_sub(tmpa, Y0, tmpb); // -X0 - Y0 mpc_sub(tmpb, Z0, tmpc); // -X0 - Y0 - Z0 mpc_div(tmpc, d2, P0); // P0 = (-X0 - Y0 - Z0) / 2 mpc_sub(A0, y, tmpb); // A0 - y mpc_mul(X0, Y0, tmpa); // X0 * Y0 mpc_mul(X0, Z0, tmpb); // X0 * Z0 mpc_mul(Y0, Z0, tmpc); // Y0 * Z0 mpc_add(tmpa, tmpb, tmpd); // X0 * Y0 + X0 * Z0 mpc_add(tmpc, tmpd, tmpa); // X0 * Y0 + X0 * Z0 + Y0 * Z0 mpc_mul(P0, P0, tmpb); // P0 * p0 mpc_mul(tmpb, d3, tmpc); // 3 * P0 * p0 mpc_sub(tmpa, tmpc, E2); // E2 = X0 * Y0 + X0 * Z0 + Y0 * Z0 - 3 * P0 * p0 mpc_mul(X0, Y0, tmpa); // X0 * Y0 mpc_mul(tmpa, Z0, tmpb); // X0 * Y0 * Z0 mpc_mul(E2, P0, tmpa); // E2 * P0 mpc_mul(tmpa, d2, tmpc); // 2 * E2 * P0 mpc_mul(P0, P0, tmpa); // P0 * P0 mpc_mul(tmpa, P0, tmpd); // P0 * P0 * p0 mpc_mul(tmpd, d4, tmpa); // 4 * P0 * P0 * p0 mpc_add(tmpb, tmpc, tmpd); // X0 * Y0 * Z0 + 2 * E2 * P0 mpc_add(tmpd, tmpa, E3); // E3 = X0 * Y0 * Z0 + 2 * E2 * P0 + 4 * P0 * P0 * p0 mpc_abs(P0, Pm); // |P0| if not s_mpf_is0(Pm) then begin // |P0| <> 0 // 分母ゼロ回避 {p0がゼロにる場合誤差が十分に小さくなっている為結果には影響しません} mpc_mul(P0, P0, tmpa); // P0 * P0 mpc_mul(tmpa, P0, tmpb); // tmpb P0 * P0 * P0 mpc_mul(tmpb, d3, tmpa); // tmpa 3 * P0 * P0 * P0 mpc_mul(E2, P0, tmpb); // tmpb E2 * P0 mpc_mul(X0, Y0, tmpc); // tmpc X0 * Y0 mpc_mul(tmpc, Z0, tmpd); // tmpd X0 * Y0 * Z0 mpc_mul(tmpd, d2, tmpc); // tmpc 2 * X0 * Y0 * Z0 mpc_add(tmpc, tmpb, tmpd); // tmpd = tmpc + tmpb 2 * X0 * Y0 * Z0 + E2 * P0 mpc_add(tmpd, tmpa, tmpb); // tmpb = tmpd + tmpa 2 * X0 * Y0 * Z0 + E2 * P0 + 3 * P0 * P0 * P0 mpc_div(tmpb, P0, E4); // E4 = tmpb / po 2 * X0 * Y0 * Z0 + E2 * P0 + 3 * P0 * P0 * P0 / P0 end else begin mpc_set_dbl(E4, 0, 0); // P0がゼロの時は誤差が小さくなっているで問題なし // ERRF := ERRF or 00100000; end; mpc_mul(X0, Y0, tmpa); // X0 * Y0 mpc_mul(tmpa, Z0, tmpb); // X0 * Y0 * Z0 mpc_mul(tmpb, P0, tmpc); // X0 * Y0 * Z0 * P0 mpc_mul(tmpc, P0, E5); // E5 = X0 * Y0 * Z0 * P0 * P0 mpc_set_dbl(tmpb, -m, 0); // -m, 0i mpc_pow(d4, tmpb, P4); // p4 = d4^-m mpc_set_dbl(tmpc, -3 / 2, 0); // -3 / 2, 0i mpc_pow(Am, tmpc, powa); // Am^(-3 / 2) mpc_mul(c1, E2, tmpa); // c1 * E2 mpc_sub(one, tmpa, def); // 1 - c1 * E2 mpc_mul(c2, E3, tmpb); // c2 * E3 mpc_add(def, tmpb, tmpc); // 1 - c1 * E2 + c2 * E3 mpc_mul(E2, E2, tmpb); // E2 * E2 mpc_mul(c3, tmpb, tmpd); // c3 * E2 * E2 mpc_add(tmpd, def, tmpa); // 1 - c1 * E2 + c2 * E3 + c3 * E2 * E2 mpc_mul(c4, E4, tmpb); // c4 * E4 mpc_mul(E2, E3, tmpc); // E2 * E3 mpc_mul(c5, tmpc, tmpd); // c5 * E2 * E3 mpc_sub(tmpa, tmpb, tmpc); // 1 - c1 * E2 + c2 * E3 + c3 * E2 * E2 - c4 * E4 mpc_sub(tmpc, tmpd, tmpa); // 1 - c1 * E2 + c2 * E3 + c3 * E2 * E2 - c4 * E4 - c5 * E2 * E3 mpc_mul(c6, E5, tmpb); // c6 * E5 mpc_add(tmpa, tmpb, def); // def = 1 - c1*E2 + c2*E3 + c3*E2*E2 - c4*E4 - c5*E2*E3 + c6*E5 mpc_mul(d6, sum, tmpa); // 6 * sum mpc_mul(powa, def, tmpb); // powa * def mpc_mul(tmpb, p4, tmpc); // d4^-m * powa * def mpc_add(tmpc, tmpa, rj); // rj = d4^-m * powa * def + 6 * sum if test < 0 then begin rfd(xt, yt, zt, tmpa); // rf(xt, yt, zt) mpc_sub(rcx, tmpa, tmpb); // rcx - rf(xt, yt, zt) mpc_mul(d3, tmpb, tmpc); // 3 * (rcx - rf(xt, yt, zt)) mpc_mul(b, rj, tmpd); // b * rj mpc_add(tmpc, tmpd, tmpa); // b * rj + 3 * (rcx - rf(xt, yt, zt)) mpc_mul(a, tmpa, rj); // rj = a * (b * rj + 3 * (rcx - rf(xt, yt, zt))) end; mpc_copy(rj, ans); form1.Canvas.TextOut(250,320,' '); form1.Canvas.TextOut(250,320,'rj 判定 ' + floatTostr(test)); if ERRF <> 0 then Form1.Memo1.Lines.Append('RJ 演算エラー有り code=' + intTostr(ERRF)); Ext: Form1.Memo1.Lines.Append('RJ Loop数 = ' + intTostr(m)); mpc_clear(one); mpc_clear5(rj, xt, yt, zt, pt); mpc_clear5(sum, a, b, rho, A0); mpc_clear4(tau, rcx, lamda, Am); mpc_clear4(sqrtx, sqrty, sqrtz, sqrtp); mpc_clear3(dm, em, delta); mpc_clear3(c1, c2, c3); mpc_clear3(c4, c5, c6); mpc_clear4(X0, Y0, Z0, P0); mpc_clear4(E2, E3, E4, E5); mpc_clear5(d0, d2, d3, d4, d5); mpc_clear4(d6, p4, powa, def); mpc_clear4(tmpa, tmpb, tmpc, tmpd); mpf_clear3(tmpaf, tmpbf, tmpcf); mpf_clear2(Q, pm); end; var rcdt : array[0..5] of array[0..1] of mp_complex; rcan : array[0..5] of array[0..1] of double; procedure datasetRc; begin Vc( 0, 0, rcdt[0,0]); Vc(1/4, 0, rcdt[0, 1]); Vc(9/4, 0, rcdt[1,0] ); Vc( 2, 0, rcdt[1, 1]); Vc( 0, 0, rcdt[2,0]); Vc( 0, 1, rcdt[2, 1]); Vc( 0, -1, rcdt[3,0]); Vc( 0, 1, rcdt[3, 1]); Vc(1/4, 0, rcdt[4,0] ); Vc( -2, 0, rcdt[4, 1]); Vc( 0, 1, rcdt[5,0]); Vc( -1, 0, rcdt[5, 1]); rcan[0,0] := 3.1415926535898; rcan[0,1] := 0; rcan[1,0] := 0.69314718055995; rcan[1,1] := 0; rcan[2,0] := 1.1107207345396; rcan[2,1] := -1.1107207345396; rcan[3,0] := 1.2260849569072; rcan[3,1] := -0.34471136988768; rcan[4,0] := 0.23104906018665; rcan[4,1] := 0; rcan[5,0] := 0.77778596920447; rcan[5,1] := 0.19832484993429; end; var rfdt : array[0..5] of array[0..2] of mp_complex; rfan : array[0..5] of array[0..1] of double; procedure datasetRf; begin Vc( 1, 0, rfdt[0,0]); Vc( 2, 0, rfdt[0,1]); Vc( 0, 0, rfdt[0,2]); Vc( 0, 1, rfdt[1,0]); Vc( 0,-1, rfdt[1,1]); Vc( 0, 0, rfdt[1,2]); Vc(-1, 1, rfdt[2,0]); Vc( 0, 1, rfdt[2,1]); Vc( 0, 0, rfdt[2,2]); Vc( 2, 0, rfdt[3,0]); Vc( 3, 0, rfdt[3,1]); Vc( 4, 0, rfdt[3,1]); Vc( 0, 1, rfdt[4,0]); Vc( 0,-1, rfdt[4,1]); Vc( 2, 0, rfdt[4,2]); Vc(-1, 1, rfdt[5,0]); Vc( 0, 1, rfdt[5,1]); Vc( 1,-1, rfdt[5,2]); rfan[0,0] := 1.3110287771461; rfan[0,1] := 0; rfan[1,0] := 1.8540746773014; rfan[1,1] := 1.8540746773014; rfan[2,0] := 0.79612586584234; rfan[2,1] := -1.2138566698365; rfan[3,0] := 0.58408284167715; rfan[3,1] := 0; rfan[4,0] := 1.0441445654064; rfan[4,1] := 0; rfan[5,0] := 0.93912050218619; rfan[5,1] := -0.53296252018635; end; var rjdt : array[0..9] of array[0..3] of mp_complex; rjan : array[0..9] of array[0..1] of double; procedure datasetRj; begin Vc( 0, 0, rjdt[0,0]); Vc( 1, 0, rjdt[0,1]); Vc( 2, 0, rjdt[0,2]); Vc( 3, 0, rjdt[0,3]); Vc( 2, 0, rjdt[1,0]); Vc( 3, 0, rjdt[1,1]); Vc( 4, 0, rjdt[1,2]); Vc( 5, 0, rjdt[1,3]); Vc( 2, 0, rjdt[2,0]); Vc( 3, 0, rjdt[2,1]); Vc( 4, 0, rjdt[2,2]); Vc(-1, 1, rjdt[2,3]); Vc( 0, 1, rjdt[3,0]); Vc( 0,-1, rjdt[3,1]); Vc( 0, 0, rjdt[3,2]); Vc( 2, 0, rjdt[3,3]); Vc(-1, 1, rjdt[4,0]); Vc(-1,-1, rjdt[4,1]); Vc( 1, 0, rjdt[4,2]); Vc( 2, 0, rjdt[4,3]); Vc( 0, 1, rjdt[5,0]); Vc( 0,-1, rjdt[5,1]); Vc( 0, 0, rjdt[5,2]); Vc( 1,-1, rjdt[5,3]); Vc(-1, 1, rjdt[6,0]); Vc(-1,-1, rjdt[6,1]); Vc( 1, 0, rjdt[6,2]); Vc(-3, 1, rjdt[6,3]); Vc(-1, 1, rjdt[7,0]); Vc(-2,-1, rjdt[7,1]); Vc( 0,-1, rjdt[7,2]); Vc(-1, 1, rjdt[7,3]); Vc( 2, 0, rjdt[8,0]); Vc( 3, 0, rjdt[8,1]); Vc( 4, 0, rjdt[8,2]); Vc(-0.5, 0,rjdt[8,3]); Vc( 2, 0, rjdt[9,0]); Vc( 3, 0, rjdt[9,1]); Vc( 4, 0, rjdt[9,2]); Vc(-5, 0, rjdt[9,3]); rjan[0,0] := 0.77688623778582; rjan[0,1] := 0; rjan[1,0] := 0.14297579667157; rjan[1,1] := 0; rjan[2,0] := 0.13613945827771; rjan[2,1] := -0.38207561624427; rjan[3,0] := 1.6490011662711; rjan[3,1] := 0; rjan[4,0] := 0.94148358841220; rjan[4,1] := 0; rjan[5,0] := 1.8260115229009; rjan[5,1] := 1.2290661908643; rjan[6,0] := -0.61127970812028; rjan[6,1] := -1.0684038390007; rjan[7,0] := 1.8249027393704; rjan[7,1] := -1.2218475784827; rjan[8,0] := 0.24723819703052; rjan[8,1] := 0; rjan[9,0] := -0.12711230042964; rjan[9,1] := 0; end; procedure TForm1.Rctest1Click(Sender: TObject); var i, j : integer; x, y : mp_complex; ans0r, ans0i: double; ans : mp_complex; begin for i := 0 to 5 do for j := 0 to 1 do begin mpc_init(rcdt[i, j]); end; mpc_init3(x, y, ans); Form1.Memo1.clear; datasetRc; i := RcRadio.ItemIndex; mpc_copy(rcdt[i,0], x); mpc_copy(rcdt[i,1], y); ans0r := rcan[i,0]; ans0i := rcan[i,1]; Rcansr.Text := floatTostr(ans0r); Rcansi.Text := floatTostr(ans0i); RCd(x, y, ans); Rcdtxr.Text := floatTostr(mpf_todouble(x.re)); Rcdtxi.Text := floatTostr(mpf_todouble(x.im)); Rcdtyr.Text := floatTostr(mpf_todouble(y.re)); Rcdtyi.Text := floatTostr(mpf_todouble(y.im)); Rcclc0.Text := floatTostr(mpf_todouble(ans.re)); Rcclc1.Text := floatTostr(mpf_todouble(ans.im)); for i := 0 to 5 do for j := 0 to 1 do begin mpc_clear(rcdt[i, j]); end; mpc_clear3(x, y, ans); end; procedure TForm1.Rctest2Click(Sender: TObject); var dinxr, dinxi: double; dinyr, dinyi: double; ch: integer; x, y: mp_complex; ans : mp_complex; begin val(Rcdtxr.Text, dinxr, ch); if ch <> 0 then begin application.MessageBox('dtxrに間違いがあります。','注意',0); exit; end; val(Rcdtxi.Text, dinxi, ch); if ch <> 0 then begin application.MessageBox('dtxiに間違いがあります。','注意',0); exit; end; val(Rcdtyr.Text, dinyr, ch); if ch <> 0 then begin application.MessageBox('dtyrに間違いがあります。','注意',0); exit; end; val(Rcdtyi.Text, dinyi, ch); if ch <> 0 then begin application.MessageBox('dtyiに間違いがあります。','注意',0); exit; end; mpc_init3(x, y, ans); memo1.Clear; Vc( dinxr, dinxi, x); Vc( dinyr, dinyi, y); RCd(x, y, ans); memo1.Lines.Append(floatTostr(mpf_todouble(ans.re)) + ' ' + floatTostr(mpf_todouble(ans.im)) + 'i'); mpc_clear3(x, y, ans); end; procedure TForm1.Rftest1Click(Sender: TObject); var i, j : integer; x, y, z: mp_complex; ans0r, ans0i: double; ans : mp_complex; begin for i := 0 to 5 do for j := 0 to 2 do begin mpc_init(rfdt[i, j]); end; mpc_init4(x, y, z, ans); Form1.Memo1.clear; datasetRf; i := RfRadio.ItemIndex; mpc_copy(rfdt[i,0], x); mpc_copy(rfdt[i,1], y); mpc_copy(rfdt[i,2], z); ans0r := rfan[i,0]; ans0i := rfan[i,1]; Rfansr.Text := floatTostr(ans0r); Rfansi.Text := floatTostr(ans0i); RFd(x, y, z, ans); Rfdtxr.Text := floatTostr(mpf_todouble(x.re)); Rfdtxi.Text := floatTostr(mpf_todouble(x.im)); Rfdtyr.Text := floatTostr(mpf_todouble(y.re)); Rfdtyi.Text := floatTostr(mpf_todouble(y.im)); Rfdtzr.Text := floatTostr(mpf_todouble(z.re)); Rfdtzi.Text := floatTostr(mpf_todouble(z.im)); Rfclc0.Text := floatTostr(mpf_todouble(ans.re)); Rfclc1.Text := floatTostr(mpf_todouble(ans.im)); for i := 0 to 5 do for j := 0 to 2 do begin mpc_clear(rfdt[i, j]); end; mpc_clear4(x, y, z, ans); end; procedure TForm1.Rftest2Click(Sender: TObject); var dinxr, dinxi: double; dinyr, dinyi: double; dinzr, dinzi: double; ch: integer; x, y, z: mp_complex; ans : mp_complex; begin val(Rfdtxr.Text, dinxr, ch); if ch <> 0 then begin application.MessageBox('dtxrに間違いがあります。','注意',0); exit; end; val(Rfdtxi.Text, dinxi, ch); if ch <> 0 then begin application.MessageBox('dtxiに間違いがあります。','注意',0); exit; end; val(Rfdtyr.Text, dinyr, ch); if ch <> 0 then begin application.MessageBox('dtyrに間違いがあります。','注意',0); exit; end; val(Rfdtyi.Text, dinyi, ch); if ch <> 0 then begin application.MessageBox('dtyiに間違いがあります。','注意',0); exit; end; val(Rfdtzr.Text, dinzr, ch); if ch <> 0 then begin application.MessageBox('dtzrに間違いがあります。','注意',0); exit; end; val(Rfdtzi.Text, dinzi, ch); if ch <> 0 then begin application.MessageBox('dtziに間違いがあります。','注意',0); exit; end; mpc_init4(x, y, z, ans); memo1.Clear; Vc( dinxr, dinxi, x); Vc( dinyr, dinyi, y); Vc( dinzr, dinzi, z); Rfd(x, y, z, ans); memo1.Lines.Append(floatTostr(mpf_todouble(ans.re)) + ' ' + floatTostr(mpf_todouble(ans.im)) + 'i'); mpc_clear4(x, y, z, ans); end; procedure TForm1.RJtest1Click(Sender: TObject); var i, j : integer; x, y, z, p: mp_complex; ans0r, ans0i: double; ans : mp_complex; begin for i := 0 to 9 do for j := 0 to 3 do begin mpc_init(rjdt[i, j]); end; mpc_init5(x, y, z, p, ans); Form1.Memo1.clear; datasetRj; i := RJRadio.ItemIndex; mpc_copy(rjdt[i,0], x); mpc_copy(rjdt[i,1], y); mpc_copy(rjdt[i,2], z); mpc_copy(rjdt[i, 3], p); ans0r := rjan[i,0]; ans0i := rjan[i,1]; RJansr.Text := floatTostr(ans0r); RJansi.Text := floatTostr(ans0i); RJd(x, y, z, p, ans); RJdtxr.Text := floatTostr(mpf_todouble(x.re)); RJdtxi.Text := floatTostr(mpf_todouble(x.im)); RJdtyr.Text := floatTostr(mpf_todouble(y.re)); RJdtyi.Text := floatTostr(mpf_todouble(y.im)); RJdtzr.Text := floatTostr(mpf_todouble(z.re)); RJdtzi.Text := floatTostr(mpf_todouble(z.im)); RJdtpr.Text := floatTostr(mpf_todouble(p.re)); RJdtpi.Text := floatTostr(mpf_todouble(p.im)); RJclc0.Text := floatTostr(mpf_todouble(ans.re)); RJclc1.Text := floatTostr(mpf_todouble(ans.im)); for i := 0 to 9 do for j := 0 to 3 do begin mpc_clear(rjdt[i, j]); end; mpc_clear5(x, y, z, p, ans); end; procedure TForm1.Rjtest2Click(Sender: TObject); var dinxr, dinxi: double; dinyr, dinyi: double; dinzr, dinzi: double; dinpr, dinpi: double; ch: integer; x, y, z, p : mp_complex; ans : mp_complex; begin val(RJdtxr.Text, dinxr, ch); if ch <> 0 then begin application.MessageBox('dtxrに間違いがあります。','注意',0); exit; end; val(RJdtxi.Text, dinxi, ch); if ch <> 0 then begin application.MessageBox('dtxiに間違いがあります。','注意',0); exit; end; val(RJdtyr.Text, dinyr, ch); if ch <> 0 then begin application.MessageBox('dtyrに間違いがあります。','注意',0); exit; end; val(RJdtyi.Text, dinyi, ch); if ch <> 0 then begin application.MessageBox('dtyiに間違いがあります。','注意',0); exit; end; val(RJdtzr.Text, dinzr, ch); if ch <> 0 then begin application.MessageBox('dtzrに間違いがあります。','注意',0); exit; end; val(RJdtzi.Text, dinzi, ch); if ch <> 0 then begin application.MessageBox('dtziに間違いがあります。','注意',0); exit; end; val(RJdtpr.Text, dinpr, ch); if ch <> 0 then begin application.MessageBox('dtprに間違いがあります。','注意',0); exit; end; val(RJdtpi.Text, dinpi, ch); if ch <> 0 then begin application.MessageBox('dtpiに間違いがあります。','注意',0); exit; end; mpc_init5(x, y, z, p, ans); Vc( dinxr, dinxi, x); Vc( dinyr, dinyi, y); Vc( dinzr, dinzi, z); Vc( dinpr, dinpi, p); memo1.Clear; RJd(x, y, z, p, ans); memo1.Lines.Append(floatTostr(mpf_todouble(ans.re)) + ' ' + floatTostr(mpf_todouble(ans.im)) + 'i'); mpc_clear5(x, y, z, p, ans); end; end.
RC_RF_RJ_checkBignum_arithmetic.zip
各種プログラム計算例に戻る