picalc.zip三角関数、逆三角関数、その他関数 に戻る
2022-03-24 ガウス＝ルジャンドルのアルゴリズムのBigDecimal計算の精度が不足していたので修正しました。
BigDecimalによるπの計算
BigDecimal には関数が用意されていないので、関数のルーチンを作成するにあたっての、演算の特徴について検討をします。
演算の特徴を調べる為、πの値を求めるプログラムを作成してみました。
答え合わせとして、MPArith多倍長でのプログラムも作成しました、mp_floatにはπを求める関数が組み込まれていますが、一応プログラムを組んでみました。
演算結果が正しいかは、インターネットでπの値を検索し確認しました(10000桁)。
検討結果
１．デフォルトの有効桁数の設定はあるが、これが動作するのは、除算と平方根のみで、加算、減算、乗算では動作しない。
２．.除算
ａ．答えが有効桁数を超える場合はデフォルト有効桁数に丸められる。
ｂ．分子の、有効桁数がデフォルトの有効桁数の二倍近い、或いは二倍以上ある場合で分母との有効桁数の差がデフォルトより大きいと除算出来ない。
ｃ．分母の有効桁数がデフォルト値に対して何倍もの大きさになると演算不可になる場合があります。
ｄ．有効桁数を個別指定した場合は、ａ，ｂのデフォルトの有効桁数となります。
３．平方根
ａ． 元の有効桁数がデフォルト値より小さく二分の一より大きい場合 元の有効桁数とデフォルト値との平均値。
ｂ．
元の有効桁数がデフォルト値の二分の一より小さい場合 デフォルト値の二分の一。
ｃ． 元の有効桁数がデフォルト値と等しいか大きい場合 元の有効桁数。
ｄ． 有効桁数に指定がある場合指定値がa,b,c,のデフォルト値。
ｅ．
平方根が開ける値の場合で、元の値の末尾にゼロを含む場合は、ゼロの数が半分になる有効桁数。
４．加減算 丸め無し。
５．累乗 a^n 有効桁数は元の有効桁数×n。
６． 小数点以下のゼロも有効桁数 a := '2.0000' <- 有効桁数 5 b
:= '123.12300' <- 有効桁数 8
累乗で、小数点以下の値の末尾にゼロを含む場合は、RemoveTrailingZeros(n)で末尾のゼロの消去をして不要に有効桁数が増えるのを防止します。
加算、減算、乗算では、丸めがないので、非常に正確に計算されます。
除算の計算を行う場合は、前もって値を演算可能な有効桁数に設定する、ことにより演算エラーの防止することが必用です。
特に乗算の後の除算には注意が必用です、有効桁数30桁と50桁の値の乗算を行うと有効桁数80桁になります、この場合小数点があっても丸めは行われません。
乗算時二つの値の有効桁数がデフォルト値になっていると、乗算後の有効桁数が倍になるので、次に除算があると、演算不可となる場合があります。
演算エラーにならないエラーで値がゼロになります。
分子と、分母が有効桁数を超えていても、分子と分母の有効桁数が同じか、分母の有効桁数が大きければ除算は正常に演算されます。
元の有効桁数或いは、デフォルト値と違う有効桁数が必要な場合は有効桁数を指定した除算をおこないます。
平方根、値がデフォルトの有効桁数に達していない場合は、有効桁数が小さくなってしまうので、値をデフォルトの有効桁数かそれ以上に設定してから、平方根を開きます。
π計算プログラム
10000桁迄、計算するプログラムを作成
BigDecimaで関数を使用する場合、MPArithには殆ど揃っているので、これを利用しても良いかと思います。
その時のデーターのやり取りは、String形式で行います。
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.Buttons, Vcl.ExtCtrls;
type
TForm1 = class(TForm)
Memo1: TMemo;
BitBtn1: TBitBtn;
LabeledEdit1: TLabeledEdit;
BitBtn2: TBitBtn;
procedure BitBtn1Click(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses Velthuis.BigDecimals, mp_real, mp_types, mp_base, system.Math;
// https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html
// T.Ooura AGM アルゴリズム
// pcs: πの精度 有効桁数
// pai: πの値
// mp_arithは演算毎に設定有効桁数数に調整されます。
procedure pi_mpf(pcs: integer; var pai: mp_float);
var
SQRT_SQRT_EPSILON, c, b, e, tmp : mp_float;
npow, n : cardinal;
a : mp_float;
begin
mpf_set_default_decprec(pcs + 10); // 指定された有効桁数+αに設定
mpf_init5(SQRT_SQRT_EPSILON, c, b, e, tmp); // 変数の初期化
mpf_init(a); // 変数の初期化
mpf_set_int(SQRT_SQRT_EPSILON, 1); // 1 設定
mpf_set_int(tmp, 2); // 2 設定
mpf_expt_int(tmp, pcs shl 0, tmp); // 2^(pcs >> 0)
mpf_div(SQRT_SQRT_EPSILON, tmp, SQRT_SQRT_EPSILON); // 1 / 2^pcs
mpf_set_dbl(c, 1 / 8); // 0.125
mpf_sqrt(c, c); // c = sqrt(0.125);
mpf_set_int(tmp, 3); // 3
mpf_mul(tmp, c, a); // 3 * c
mpf_add_int(a, 1, a); // a = 1 + c * 3;
mpf_sqrt(a, b); // b = sqrt(a)
mpf_set_dbl(e, 5 / 8); // 0.625
mpf_sub(b, e, e); // e = b - 0.625
mpf_add(b, b, b); // b = 2 * b
mpf_sub(e, c, c); // c = e - c
mpf_add(a, e, a); // a = a + e
npow := 4;
n := 0;
repeat
npow := npow + npow; // npow = 2 * npow
mpf_add(a, b, tmp); // 'e = a + b'
mpf_div_int(tmp, 2, e); // e = (a + b) / 2
mpf_mul(a, b, tmp); // 'b = a * b'
mpf_sqrt(tmp, b); // b = sqrt(a * b)
mpf_sub(e, b, e); // e = e - b
mpf_add(b, b, b); // b = 2 * b
mpf_sub(c, e, c); // c = c - e
mpf_add(e, b, a); // a := e + b
inc(n);
until mpf_cmp_mag(e, SQRT_SQRT_EPSILON) <= 0; // while e > SQRT_SQRT_EPSILON
mpf_mul(e, e, e);
mpf_div_int(e, 4, e); // e = e * e / 4
mpf_add(a, b, a); // a = a + b
mpf_mul(c, a, c); // 'c = a * c
mpf_sub(c, e, c); // 'c = a * c - e
mpf_mul(a, a, a); // 'a = a * a
mpf_sub(a, e, a); // 'a = a * a - e
mpf_div_int(e, 2, e); // 'e = e / 2
mpf_sub(a, e, a); // 'a = a * a - e - e / 2
mpf_div(a, c, a); // 'a = (a * a - e - e / 2)/(a * c - e)
mpf_div_int(a, npow, a); // a = (a * a - e - e / 2)/(a * c - e) / npow
mpf_copy(a, pai);
mpf_clear(a); // 変数の解放
mpf_clear5(SQRT_SQRT_EPSILON, c, b, e, tmp); // 変数の解放
Form1.Memo1.Lines.Append('T.Ooura mpf n = ' + intTostr(n));
end;
// https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html
// T.Ooura AGM アルゴリズム
// pcs πの精度
// *Bigdecimal 使用時注意点
// *各演算での有効桁数の値に注意が必要です。
// 1 乗算 有効桁数が 加算される。
// 2 除算
// a デフォルトの有効桁数。
// b 割り切れる場合は割り切れた値の桁数が有効桁数。
// c 分子の有効桁数がデフォルト値より大きくて、分母の有効桁数がデフォルトより小さい場合演算不可あり
// その場合、分子か、分母の有効桁数をデフォルト値に設定する必要あり。
// 3 平方根
// a 元の有効桁数がデフォルト値より小さく二分の一より大きい場合 元の有効桁数とデフォルト値との平均値。
// b 元の有効桁数がデフォルト値の二分の一より小さい場合 デフォルト値の二分の一。
// c 元の有効桁数がデフォルト値と等しいか大きい場合 元の有効桁数。
// d 有効桁数に指定がある場合指定値がa,b,c,のデフォルト値。
// e 平方根が開ける値の場合で、元の値の末尾にゼロを含む場合は、ゼロの数が半分になる有効桁数。
// 4 加減算 有効桁数が一番大きい値の有効桁数。
// 5 累乗 a^n 有効桁数は元の有効桁数×n。
// 6 小数点以下のゼロも有効桁数 a := '2.0000' <- 有効桁数 5
// 7 乗算、累乗で、小数点以下の値の末尾にゼロを含む場合は、RemoveTrailingZeros(n)で末尾のゼロの消去
// 特に、累乗の場合不必要に有効桁数が大きくなるのを防止します。
// 8 収束判定が必要な場合は、丸めRoundToPrecision(n)指定が必要。
function big_pi(pcs: integer): Bigdecimal;
var
n, dpcs, pcsBack : integer;
rmback : BigDecimal.RoundingMode;
SQRT_SQRT_EPSILON, c, b, e : BigDecimal;
npow : Bigdecimal;
a, one, two, four, five, eight : Bigdecimal;
begin
pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ
rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ
dpcs := pcs + 10; // 指定精度 + α
BigDecimal.DefaultPrecision := dpcs;
BigDecimal.DefaultRoundingMode := rmNearesteven;
one := '1';
one := one.RoundToPrecision(dpcs); // oneの有効桁数をdpcsに設定 sqrtでの有効桁数制御
two := '2';
four := '4';
five := '5';
eight := '8';
SQRT_SQRT_EPSILON := one / BigDecimal.IntPower(two, pcs shl 0); // 収束判定値
c := BigDecimal.Sqrt(one / eight, dpcs * 2); // 平方根有効桁数での丸め
a := one + c + c + c;
b := BigDecimal.Sqrt(a, dpcs * 2);
e := b - five / eight;
b := b + b;
c := e - c;
a := a + e;
npow := '4';
n := 0;
while e > SQRT_SQRT_EPSILON do begin
npow := npow + npow;
e := (a + b) / two; // 平均値 除算有効桁数での丸め
b := a * b;
b := b.RoundToPrecision(dpcs); // pcs + α 丸め sqrt前の有効桁数とwhile内の有効桁数制御
b := BigDecimal.Sqrt(b, dpcs * 2); // 平方根有効桁数での丸め
e := e - b;
b := b + b;
c := c - e;
a := e + b;
inc(n);
end;
e := e * e;
e := e.RoundToPrecision(dpcs); // pcs + α 丸め 次の除算の為の有効桁数設定
e := e / four; // e = e + e / 4
a := a + b;
result := (a * a - e - e / two) / (a * c - e) / npow; // 除算の順番に注意 ()/() 有効桁数同じで除算有効桁数 /npow
result := result.RoundToPrecision(pcs); // 指定の精度に丸め
BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰
BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰
Form1.Memo1.Lines.Append('T.Ooura Big n = ' + intTostr(n));
end;
// ガウス＝ルジャンドルのアルゴリズム
// pcs πの精度
function pi_big(pcs: integer): Bigdecimal;
var
a, b, t, x, an, two, one, four, EPSILON : BigDecimal;
n, dpcs, pcsBack : integer;
rmback : BigDecimal.RoundingMode;
begin
pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ
rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ
dpcs := pcs + 10; // 指定精度 * 2
BigDecimal.DefaultPrecision := dpcs;
BigDecimal.DefaultRoundingMode := rmNearesteven;
one := '1';
one := one.RoundToPrecision(dpcs); // oneの有効桁数をdpcsに設定 a, b, t, EPCILON の有効桁数になる
two := '2';
EPSILON := one / BigDecimal.IntPower(two, pcs shl 1); // 収束判定値
four := '4';
a := one;
b := one / BigDecimal.Sqrt(two, dpcs * 2); // 平方根有効桁数での丸め
t := one / four; // 除算精度
x := '1';
n := 0;
repeat
an := (a + b) / two; // 平均値 除算有効桁数での丸め
b := a * b;
b := b.RoundToPrecision(dpcs); // pcs + α 丸め 平方根前の丸め
b := BigDecimal.Sqrt(b, dpcs * 2); // 平方根有効桁数での丸め
t := t - x * (an - a) * (an - a);
t := t.RoundToPrecision(dpcs); // pcs + α 丸め until後でもok
x := x + x;
a := an;
inc(n);
until (BigDecimal.abs(a - b) < EPSILON) or (n > 100); // n>100はデバッグ用
a := a + b;
result := a * (a / (four * t));
result := result.RoundToPrecision(pcs); // 指定の精度に丸め
BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰
BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰
Form1.Memo1.Lines.Append('ガウス＝ルジャンドル Big n = ' + intTostr(n));
end;
// ガウス＝ルジャンドルのアルゴリズム
function pi_calc: double;
var
a, b, t, x, an: double;
n : integer;
EPSILON : double;
begin
EPSILON := 1 / power(2, 17); // 収束判定値の設定 17桁
a := 1;
b := 1 / sqrt(2);
t := 1 / 4;
x := 1;
n := 0;
repeat
an := (a + b) / 2; // 平均値
b := sqrt(a * b);
t := t - x * (an - a) * (an - a);
x := x * 2;
a := an;
inc(n);
until (abs(a - b) < EPSILON) or (n > 100); // n>100はデバッグ用
result := (a + b) * (a + b) / (4 * t);
Form1.Memo1.Text := 'ガウス＝ルジャンドル Double n = ' + intTostr(n);
end;
procedure TForm1.BitBtn1Click(Sender: TObject);
const
PREC = 10000;
MPRC = 3;
var
pai: double;
pai_big : BigDecimal;
s, ch : integer;
pai_mpf : mp_float;
pcsBack : integer;
begin
val(labelededit1.Text, s, ch);
if ch <> 0 then begin
memo1.Text := '精度の値に間違いがあります';
exit;
end;
if (s < MPRC) or (s > PREC) then begin
memo1.Text := '精度の値は' + intTostr(MPRC) + '以上' + intTostr(PREC) + '以下にしてください。';
exit;
end;
memo1.Width := 820; // 一行100桁表示に設定
memo1.Font.Size := 11; // 一行100桁表示にフォントサイズ設定
memo1.Clear;
pai := pi_calc; // double ガウス＝ルジャンドルのアルゴリズム
memo1.Lines.Append(floatTostr(pai));
pai_big := pi_big(s); // BigDecimal ガウス＝ルジャンドルのアルゴリズム
pai_big := BigDecimal(pai_big).RemoveTrailingZeros(0); // 末尾のゼロ消去
memo1.Lines.Append(pai_big.ToString);
pai_big := big_pi(s); // BigDecimal T.Ooura AGM アルゴリズム
pai_big := BigDecimal(pai_big).RemoveTrailingZeros(0); // 末尾のゼロ消去
memo1.Lines.Append(pai_big.ToString);
pcsBack := mpf_get_default_prec; // 除算演算精度バックアップ
mpf_set_default_decprec(s);
mpf_init(pai_mpf);
pi_mpf(s, pai_mpf); // mp_arith T.Ooura AGM アルゴリズム
memo1.Lines.Append(string(mpf_adecimal_alt(pai_mpf, s)));
mpf_clear(pai_mpf);
mpf_set_default_prec(pcsBack); // 演算精度の復帰
end;
end.
picalc.zip
三角関数、逆三角関数、その他関数 に戻る