左は
Math_big.pas のテストプログラムの実行画面です、値の確認、デバックに使用しています。BigDecimalのデフォルトの有効桁数は、64桁です。
有効桁数が大きくなると、値が収束するのに時間が掛かります、一応2500桁迄は値の確認はしています。
2500桁を超える場合は、プログラムの中の収束計算のループ数上限の値を変更して下さい。
2022-12-13 BigdecimalToDouble の変換を追加しました。
2022-07-06 e^x のxの値が大きくなると誤差が増えるのを修正しました。
2022-06-11 arsinh, arcosh,
artanhの計算計算方法変更、arsinhとarcoshの値が大きいとき(1e30以上)の時正しく計算されないのを修正。
2022-05-22 Log(1)の0値の有効桁数が正しくなるように修正しました。
2022-05-08 sin,cos,tanの無駄な角度計算を省略しました。
2022-05-02 degTorad radTodeg 関数を追加して、変換の精度を上げました。
2022-05-01 数値入力チェック方法変更
2022-04-23 tan(z) 計算方法変更
2022-04-22 tan(z)、 RCによる arcsin, arccos計算の修正
BigDecimalによるMath
BigdecimalにはMathが無い為、基本的なものを作成してみました。
Bigdecimalの組み込みはベルヌーイ数 その4を参照して下さい。
Big_Math.pasをパスの通ったフォルダに置き、uses にBig_mathを追加すれば使用できます。
Velthuis.Bigdecimals.pasと同じフォルダーに置くのが良いと思います。
使用出来る関数は
pi_big, sin_big, cos_big, tan_big, arcsin_big,
arccos_big,, arctan_big, log_big, exp_big, arctan2_big, sinh_big, cosh_big,
tanh_big, arsinh_big, arcosh_big, artanh_big, power_big
逆数にしたりして簡単に変換出来るものは用意していません。
degTorad radTodeg
角度rad、deg変換用
bigdecimalTodouble 多倍長からDouble の変換
Bigdecimal は 当然 double より大きな値、又、0に近い値が扱えるので、maxdouble より大きな値は maxdouble、
mindoubleより小さな値はmindoubleに変換するようにしました。
doubleの場合maxdoubleより大きい値は、infinity(無限大)となるのですが、 演算エラーを防ぐ為です。
doubleの値を超える場合は注意して使用して下さい。
プログラムは64ビットモードでコンパイルしないと、演算が遅くなります。
演算精度のextended(80ビット)は無くなり、Double(64ビット)となりますが、BigInteger,
Bigdecimal演算を使用すれば全く問題ないと思います。
有効桁数は、全ての関数計算に反映されるので、50桁の実質的な有効桁数が必要な場合は、64桁程度で計算する必要があります。
デフォルト値が、64桁となっているので、50桁程度が保証される値となります。
関数は通常の関数名に"_big"を追加しています。
arcsinh, arccosh, arctanhはそれぞれ正式名の arsinh, arcosh,
artanhです。
log_bigは自然対数です。
左は
Math_big.pas のテストプログラムの実行画面です、値の確認、デバックに使用しています。
BigDecimalのデフォルトの有効桁数は、64桁です。
有効桁数が大きくなると、値が収束するのに時間が掛かります、一応2500桁迄は値の確認はしています。
2500桁を超える場合は、プログラムの中の収束計算のループ数上限の値を変更して下さい。
プログラム
Main.pas
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, system.Character,
system.Diagnostics;
type
TForm1 = class(TForm)
Memo1: TMemo;
GroupBox1: TGroupBox;
BitBtn1: TBitBtn;
LabeledEdit1: TLabeledEdit;
RadioGroup1: TRadioGroup;
ClearBtn: TBitBtn;
LabeledEdit2: TLabeledEdit;
GroupBox2: TGroupBox;
RadioGroup2: TRadioGroup;
BitBtn2: TBitBtn;
LabeledEdit3: TLabeledEdit;
GroupBox3: TGroupBox;
RadioGroup3: TRadioGroup;
BitBtn3: TBitBtn;
LabeledEdit4: TLabeledEdit;
LabeledEdit5: TLabeledEdit;
GroupBox4: TGroupBox;
RadioGroup4: TRadioGroup;
BitBtn4: TBitBtn;
LabeledEdit6: TLabeledEdit;
GroupBox5: TGroupBox;
RadioGroup5: TRadioGroup;
BitBtn5: TBitBtn;
LabeledEdit7: TLabeledEdit;
CheckBox1: TCheckBox;
Edit1: TEdit;
procedure BitBtn1Click(Sender: TObject);
procedure ClearBtnClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
procedure BitBtn3Click(Sender: TObject);
procedure BitBtn4Click(Sender: TObject);
procedure BitBtn5Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses Big_Math, Velthuis.Bigdecimals, Velthuis.bigintegers;
// 有効桁数入力
function InPrecision(s: string; m: integer): integer;
var
ch : integer;
begin
val(s, Result, ch);
if ch <> 0 then begin
application.MessageBox('有効桁数に間違いがあります。','注意',0);
result := -1;
exit;
end;
if Result < 1 then begin
application.MessageBox('有効桁数は１以上にしてください。','注意',0);
end;
if m = 0 then begin
if Result >2500 then begin
application.MessageBox('有効桁数は2500迄にしてください。','注意',0);
result := -1;
end;
end
else
if Result > 50000 then begin
application.MessageBox('有効桁数は50000迄にしてください。','注意',0);
result := -1;
end;
end;
// 指数の無い場合の数値確認
// x: 'e','E'のない数字文字列
// ca: 0 '.' がある場合 1 ない場合
// 問題が無い場合 0を返します
function numeralfloatcheck(x: string; ca: byte): integer;
var
len, dn, n, d : integer;
num, numb : string;
flg : boolean;
snum : integer;
begin
len := length(x); // 文字長さ
snum := 1;
num := x[1]; // １文字目
dn := 0;
flg := true;
if num = '+' then inc(snum); // '+' パス
if num = '-' then inc(snum); // '-' パス
if num = '.' then begin // '.'
inc(dn); // '.'カウントインクリメント
inc(snum); // パス
end;
for n := snum to len do begin // snum からlen 迄数値確認
num := x[n];
if num = '.' then begin // '.'
inc(dn); // '.'カウントインクリメント
continue; // パス
end;
flg := false;
for d := 0 to 9 do begin // 0～9の値か確認
numb := inttostr(d);
if numb = num then begin
flg := true;
break;
end;
end;
if flg = false then break;
end;
if ca = 0 then begin // '.'あり条件だったら
if dn > 1 then flg := false; // '.'1個迄なら0k
end
else // '.'無条件だったら
if dn > 0 then flg := false; // '.'があったらMG
if flg then result := 0
else result := 1;
end;
// 入力確認
function incheck(s, m: string): boolean;
const
ec = 'e';
eo = 'E';
dt = '.';
var
fs, es : string;
ch, p, po, leng : integer;
begin
result := false;
p := pos(ec, s); // 'e' の位置
po := pos(eo, s); // 'E' の位置
p := p + po;
if p = 0 then begin // 'e' 'E' が無かったら
ch := numeralfloatcheck(s, 0); // '.'がある条件で確認無くても良い
if ch <> 0 then begin
application.MessageBox(' 入力値に間違いがあります。','注意',0);
exit;
end;
result := true;
exit;
end;
leng := length(s);
fs := '';
fs := copy(s, 1, p - 1); // 'e'の前の文字取り出し
es := '';
es := copy(s, p + 1, leng - p); // 'e'の後ろの文字取り出し
ch := numeralfloatcheck(fs, 0);; // 'e'の前の文字の確認 '.'がある条件で確認無くても良い
if ch <> 0 then begin
application.MessageBox(' 入力値に間違いがあります。','注意',0);
exit;
end;
if es <> '' then begin // 'e'の後ろの文字確認
ch := numeralfloatcheck(es, 1);; // '.'が無い条件で確認
if ch <> 0 then begin
application.MessageBox(' 入力値に間違いがあります。','注意',0);
exit;
end;
end;
result := true;
end;
// sin(x),cos(x),tan(x)
procedure TForm1.BitBtn1Click(Sender: TObject);
var
Precision : integer;
x, y, deg : bigdecimal;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not incheck(LabeledEdit1.Text, '角度の') then exit;
Precision := InPrecision(LabeledEdit2.Text, 0); // 有効桁数入力
if Precision < 1 then exit;
BigDecimal.DefaultPrecision := Precision;
deg := LabeledEdit1.Text;
if deg > 720 then begin
application.MessageBox(' 角度は720°以下にして下さい。','注意',0);
exit;
end;
StopWatch := TStopwatch.StartNew;
x := degTorad(deg);
case Radiogroup1.ItemIndex of
0 : y := sin_big(x);
1 : y := cos_big(x);
2 : y := tan_big(x);
end;
memo1.Lines.Append(y.ToPlainString);
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Edit1.text := '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
end;
// arcsin(x),arccos(x),arctan(x),arctan2(y,x);
procedure TForm1.BitBtn2Click(Sender: TObject);
var
xi, yi, x, one, ans, y : bigdecimal;
zero : bigDecimal;
Precision : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not incheck(LabeledEdit3.Text, 'xの') then exit;
Precision := InPrecision(LabeledEdit2.Text, 0); // 有効桁数入力
if Precision < 1 then exit;
StopWatch := TStopwatch.StartNew;
BigDecimal.DefaultPrecision := Precision;
zero := '0';
one := '1';
D180 := '180';
xi := LabeledEdit3.Text;
x := BigDecimal.Abs(xi);
case Radiogroup2.ItemIndex of
0: begin
if x > one then begin
application.MessageBox('xの入力値が大きすぎます。','注意',0);
exit;
end;
if checkbox1.Checked = false then
y := arcsin_big(xi)
else
y := arcsina_big(xi);
end;
1: begin
if x > one then begin
application.MessageBox('xの入力値が大きすぎます。','注意',0);
exit;
end;
if checkbox1.Checked = false then
y := arccos_big(xi)
else
y := arccosa_big(xi);
end;
2: begin
if checkbox1.Checked = false then
y := arctan_big(xi)
else
y := arctana_big(xi);
end;
3: begin
if not incheck(LabeledEdit4.Text, 'yの') then exit;
yi := LabeledEdit4.Text;
if checkbox1.Checked = false then
y := arctan2_big(yi, xi)
else
y := arctan2a_big(yi, xi);
end;
4: begin
if xi < zero then begin
application.MessageBox('xゼロより大きくして下さい。','注意',0);
exit;
end;
if not incheck(LabeledEdit4.Text, 'yの') then exit;
yi := LabeledEdit4.Text;
y := power_big(xi, yi);
memo1.Lines.Append(y.ToString);
exit;
end;
end;
ans := radTodeg(y);
memo1.Lines.Append(ans.ToPlainString);
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Edit1.text := '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
end;
// log(x),e^(x),π
procedure TForm1.BitBtn3Click(Sender: TObject);
var
x, y: bigdecimal;
Precision : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
case Radiogroup3.ItemIndex of
0, 1: begin
if not incheck(LabeledEdit5.Text, 'xの') then exit;
x := LabeledEdit5.Text;
end;
end;
if Radiogroup3.ItemIndex <> 2 then
Precision := InPrecision(LabeledEdit2.Text, 0) // 有効桁数入力
else
Precision := InPrecision(LabeledEdit2.Text, 1); // 有効桁数入力
if Precision < 1 then exit;
StopWatch := TStopwatch.StartNew;
BigDecimal.DefaultPrecision := Precision;
case Radiogroup3.ItemIndex of
0: begin
y := log_big(x);
memo1.Lines.Append(y.ToPlainString);
end;
1: begin
y := exp_big(x);
memo1.Lines.Append(y.ToString);
end;
2: begin
y := pi_big;
y := y.RoundToPrecision(Precision);
memo1.Lines.Append(y.ToString);
end;
end;
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Edit1.text := '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
end;
// sinh(x),cosh(x),tanh(x)
procedure TForm1.BitBtn4Click(Sender: TObject);
var
x, y : bigdecimal;
Precision : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not incheck(LabeledEdit6.Text, 'xの') then exit;
Precision := InPrecision(LabeledEdit2.Text, 0); // 有効桁数入力
if Precision < 1 then exit;
StopWatch := TStopwatch.StartNew;
x := LabeledEdit6.Text;
BigDecimal.DefaultPrecision := Precision;
case Radiogroup4.ItemIndex of
0: y := sinh_big(x);
1: y := cosh_big(x);
2: y := tanh_big(x);
end;
memo1.Lines.Append(y.ToPlainString);
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Edit1.text := '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
end;
// arsinh(x),arcosh(x),artanh(x)
procedure TForm1.BitBtn5Click(Sender: TObject);
var
x, y, one : bigdecimal;
Precision : integer;
StopWatch : TStopwatch;
ElapsedMillseconds : Int64;
begin
if not incheck(LabeledEdit7.Text, 'xの') then exit;
Precision := InPrecision(LabeledEdit2.Text, 0); // 有効桁数入力
if Precision < 1 then exit;
StopWatch := TStopwatch.StartNew;
BigDecimal.DefaultPrecision := Precision;
one := '1';
x := LabeledEdit7.Text;
case Radiogroup5.ItemIndex of
0: y := arsinh_big(x);
1: begin
if x < one then begin
application.MessageBox('xの入力値を１以上にして下さい。','注意',0);
exit;
end;
y := arcosh_big(x);
end;
2: begin
if x >= one then begin
application.MessageBox('xの入力値を１以下にして下さい。','注意',0);
exit;
end;
y := artanh_big(x);
end;
end;
memo1.Lines.Append(y.ToPlainString);
ElapsedMillseconds := StopWatch.ElapsedMilliseconds;
Edit1.text := '演算時間 = ' + intTostr(ElapsedMillseconds) + 'msec';
end;
procedure TForm1.ClearBtnClick(Sender: TObject);
begin
Memo1.Clear;
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
memo1.Clear;
end;
end.
Big_Math.pas
// pi_bigは一万桁でもok その他は1000桁程度迄でそれを超えると、演算が非常に遅くなります。 // 有効桁数100桁程度であれば問題ありません。 // pi_big以外は2500桁位が上限です、それ以上必要な場合は収束ループの上限設定を外して下さい。 // 64ビットでコンパイルしないと演算が遅くなります。 unit Big_Math; interface uses Velthuis.Bigdecimals, System.SysUtils, Vcl.Forms, Velthuis.BigIntegers; //uses Velthuis.Bigdecimals, System.SysUtils, Vcl.Forms, Velthuis.BigIntegers, main; function pi_big: Bigdecimal; function sin_big(x: bigdecimal): bigDecimal; function cos_big(x: bigdecimal): bigdecimal; function tan_big(x: bigdecimal): bigdecimal; function arcsin_big(z: bigdecimal): bigdecimal; function arcsina_big(xi: bigdecimal): bigdecimal; function arccos_big(z: bigdecimal): bigdecimal; function arccosa_big(xi: bigdecimal): bigdecimal; function arctan_big(z: bigdecimal): bigdecimal; function arctana_big(xi: bigdecimal): bigdecimal; function log_big(xi: bigdecimal): bigdecimal; function exp_big(x: bigDecimal): bigdecimal; function arctan2_big(yi, xi: bigdecimal): bigdecimal; function arctan2a_big(yi, xi: bigdecimal): bigdecimal; function sinh_Big(x: bigdecimal): bigdecimal; function cosh_big(x: bigdecimal): bigdecimal; function tanh_big(x: bigdecimal): bigdecimal; function arsinh_big(z: bigdecimal): bigdecimal; function arcosh_big(z: bigdecimal): bigdecimal; function artanh_big(z: bigdecimal): bigdecimal; function power_big(xi, yi: bigdecimal): bigdecimal; function degTorad(deg: bigdecimal): bigdecimal; function radTodeg(rad: bigdecimal): bigdecimal; function BigdecimalToDouble(b: bigdecimal): double; implementation // https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html // T.Ooura AGM アルゴリズム function pi_big: 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 := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; one := '1'; one := one.RoundToPrecision(dpcs); // oneの有効桁数をdpcsに設定 two := '2'; four := '4'; five := '5'; eight := '8'; SQRT_SQRT_EPSILON := one / BigDecimal.IntPower(two, dpcs 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) and (n < 100) do begin npow := npow + npow; e := (a + b) / two; // 平均値 除算有効桁数での丸め b := a * b; b := b.RoundToPrecision(dpcs); // pcs + α 丸め 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; // 除算の順番に注意 result := result.RoundToPrecision(dpcs); // 指定の精度に丸め BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // deg単位をradに変換 function degTorad(deg: bigdecimal): bigdecimal; var pib, tmp, d180 : bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; D180 := '180'; tmp := deg * pib; tmp := tmp.RoundToPrecision(dpcs); result := tmp / D180; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // rad単位をdegに変換 function radTodeg(rad: bigdecimal): bigdecimal; var pib, tmp, d180 : bigdecimal; dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; dpcs := BigDecimal.DefaultPrecision + 5; // 演算精度 rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ BigDecimal.DefaultRoundingMode := rmNearesteven; D180 := '180'; tmp := rad * D180; tmp := tmp.RoundToPrecision(dpcs); result := tmp / pib; BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // sin(x) 実行部 // マクローリン展開 // x 角度 rad dpcs 有効桁数 function sin_bigf(x : BigDecimal; dpcs : integer): BigDecimal; var v, d, sq, u, k, two, one : BigDecimal; c : integer; begin x := x.RoundToPrecision(dpcs); one := '1'; Two := '2'; d := x; sq := x * x; // sqは小数点以下 dpcs * 2 桁 sq := sq.RoundToPrecision(dpcs); // 有効桁数dpcs u := d; k := 0; c := 0; repeat k := k + Two; v := u; d := -d * sq; // 小数点以下 dpcs×2 桁 d := d.RoundToPrecision(dpcs); // 有効桁数dpcs に丸めないと次の除算不可分母は小数点以下無し整数 d := d / (k * (k + one)); // 除算で有効桁数dpcs桁 除算Default値 u := v + d; u := u.RoundToPrecision(dpcs); // 有効桁数dpcs 丸め 収束判定の為 inc(c); until (v = u) or (c > 10000); v := v.RoundToPrecision(dpcs); // 有効桁数dpcs result := v; end; // sin(x)計算 BigDecimal // マクローリン展開 function sin_big(x : BigDecimal): BigDecimal; var pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; result := sin_bigf(x, dpcs); result := result.RoundToScale(pcsBack, rmNearesteven); // 有効桁数pcsBack // result := result.RoundToPrecision(pcsBack); // 有効桁数pcsBack BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // cos // マクローリン展開 // x (rad) dpcs 有効桁数 function cos_bigf(x: bigdecimal; dpcs: integer): bigdecimal; var d, s, e, i, two, one, im: bigdecimal; k : integer; begin x := x.RoundToPrecision(dpcs); one := '1'; two := '2'; e := '1'; s := '1'; i := '1'; x := x * x; x := x.RoundToPrecision(dpcs); // 有効桁数dpcs k := 0; repeat d := s; e := e.Negate(e); e := e * x; im := i + one; im := im * i; im := im.RoundToPrecision(dpcs); e := e / im; s := s + e; s := s.RoundToPrecision(dpcs); // 有効桁数dpcs 収束判定の為 i := i + two; inc(k); until (s = d) or (k > 10000); s := s.RoundToPrecision(dpcs); // 有効桁数dpcs result := s; end; // cos // x (rad) function cos_big(x: bigdecimal): bigdecimal; var pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmdown; result := cos_bigf(x, dpcs); result := result.RoundToScale(pcsBack, rmdown); // 有効桁数pcsBack // result := Result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // tan(x) = sin(x) / cos(x) function tan_big(x: bigdecimal): bigdecimal; label EXT; var pib, pis2, one, two, s, c, t, zero, xb: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; zero := '0'; one := '1'; two := '2'; pis2 := pib / Two; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; x := x.RoundToPrecision(dpcs); xb := x; // 分母0判定 repeat // 0～πの範囲に修正 if x < 0 then x := x + pib; if x >= pib then x := x - pib; until (x < pib) and (x >= 0); s := pis2 - x; if pcsback > 3 then s := s.RoundToScale(pcsback - 2, rmdown); if s = zero then begin application.MessageBox('tan_big(x) 分母が(0)無効な値です。','注意'); result := '0'; goto EXT; end; s := sin_bigf(xb, dpcs); c := cos_bigf(xb, dpcs); t := s / c; t := t.RoundToPrecision(pcsBack); // 有効桁数pcsBack result := t; EXT: BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // sinh // x function sinh_big(x: bigdecimal): bigdecimal; var v, d, sq, u, k, one, two :bigdecimal; i : integer; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; x := x.RoundToPrecision(dpcs); one := '1'; two := '2'; d := x; sq := x * x; sq := sq.RoundToPrecision(dpcs); u := d; k := '0'; i := 1; repeat v := u; k := k + two; d := d * (sq / k / (k + one)); d := d.RoundToPrecision(dpcs); u := v + d; u := u.RoundToPrecision(dpcs); inc(i) until (v = u) or (i > 10000); v := v.RoundToPrecision(pcsBack); // 有効桁数pcsBack result := v; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // cosh // x function cosh_big(x: bigdecimal): bigdecimal; var d, s, e, k, one, two : bigdecimal; i : integer; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; x := x.RoundToPrecision(dpcs); one := '1'; two := '2'; e := one; s := one; k := one; i := 0; repeat d := s; e := e * x * x; e := e.RoundToPrecision(dpcs); e := e / (k * (k + one)); s := s + e; s := s.RoundToPrecision(dpcs); k := k + two; inc(i); until (s = d) or (i > 10000); s := s.RoundToPrecision(pcsBack); // 有効桁数pcsBack result := s; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // tanh // x function tanh_big(x: bigdecimal): bigdecimal; var v, d, sq, u, k, one, two :bigdecimal; s, e: bigdecimal; i : integer; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; x := x.RoundToPrecision(dpcs); one := '1'; two := '2'; d := x; sq := x * x; sq := sq.RoundToPrecision(dpcs); u := d; k := '0'; i := 1; repeat v := u; k := k + two; d := d * (sq / k / (k + one)); d := d.RoundToPrecision(dpcs); u := v + d; u := u.RoundToPrecision(dpcs); inc(i) until (v = u) or (i > 10000); e := one; s := one; k := one; i := 0; repeat d := s; e := e * x * x; e := e.RoundToPrecision(dpcs); e := e / (k * (k + one)); s := s + e; s := s.RoundToPrecision(dpcs); k := k + two; inc(i); until (s = d) or (i > 10000); result := v / s; result := result.RoundToPrecision(pcsBack); // 有効桁数pcsBack BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // カールソンの楕円積分 RC // // ---------------- RC --------------------------------------------------------- function RC(x, y: bigdecimal): bigdecimal; label Ext; //const // r = 1E-24; // Qc = 872; // power(3 * r, -1 / 8) = 871.68554 // c1 = 3 / 10; // c2 = 1 / 7; // c3 = 3 / 8; // c4 = 9 / 22; // c5 = 159 / 208; // c6 = 9 / 8; var xt, yt, w, yb, zero : bigdecimal; lamda, A0, Q, dp : bigdecimal; s, Am, Qm : bigdecimal; c, c1, c2, c3, c4, c5, c6, Qc : bigdecimal; m, dpcs, dpcm : integer; d1, d2, d3, d4, d7, d8, d9, d10, d22, d159, d208: bigdecimal; begin dpcs := BigDecimal.DefaultPrecision; dpcm := dpcs + 5; BigDecimal.DefaultPrecision := dpcm; dp := '1E-' + intTostr(round(dpcm / 6)); zero := '0'; if y = zero then begin result := zero; goto Ext; end; d1 := '1'; d2 := '2'; d3 := '3'; d4 := '4'; d7 := '7'; d8 := '8'; d9 := '9'; d10 := '10'; d22 := '22'; d159 := '159'; d208 := '208'; c1 := d3 / d10; c2 := d1 / d7; c3 := d3 / d8; c4 := d9 / d22; c5 := d159 / d208; c6 := d9 / d8; Qc := d1 / dp; if y > zero then begin xt := x; yt := y; w := d1; end else begin xt := x - y; yt := - y; w := bigdecimal.sqrt(x / xt, dpcm * 2); // sqrt は有効桁数が1/2 end; yb := yt; A0 := (xt + yt + yt) / d3; Am := A0; Q := Qc * bigdecimal.abs(A0 - xt); // Q := power(3 * r, -1 / 8) * abs(A0 - xt); m := 0; Qm := d1; repeat lamda := d2 * bigdecimal.sqrt(xt, dpcm * 2) * bigdecimal.sqrt(yt, dpcm * 2) + yt; // sqrt は有効桁数が1/2 lamda := lamda.RoundToPrecision(dpcm); xt := (xt + lamda) / d4; yt := (yt + lamda) / d4; Am := (Am + lamda) / d4; m := m + 1; Qm := Qm / d4; until (Qm * Q < bigdecimal.abs(Am)) or (m > 10000); s := (yb - A0) * (Qm / Am); c5 := c5 + s * c6; c5 := c5.RoundToPrecision(dpcm); c4 := c4 + s * c5; c4 := c4.RoundToPrecision(dpcm); c3 := c3 + s * c4; c3 := c3.RoundToPrecision(dpcm); c2 := c2 + s * c3; c2 := c2.RoundToPrecision(dpcm); c1 := c1 + s * c2; c1 := c1.RoundToPrecision(dpcm); c := c1 * s * s + d1; c := c.RoundToPrecision(dpcm); Am := bigdecimal.sqrt(Am, dpcm * 2); // sqrt は有効桁数が1/2 result := w * (c / Am); result := result.RoundToPrecision(dpcs); // s := (yb - A0) / power(4, m) / Am; // result := w * (d1 + s * s * (c1 + s * (c2 + s * (c3 + s * (c4 + s * (c5 + s * c6)))))) / bigdecimal.sqrt(Am); Ext: BigDecimal.DefaultPrecision := dpcs; end; // arcsin // 0.7を超えたら反対側の角度を求めます(arccos計算になります) // arcsin(x) = π/2 - arccos(x) function arcsina_big(xi: bigdecimal): bigdecimal; label EXT; var i : integer; pib, ans, bans, fa, fb, d, ab, ta, tb, f, x, seven, zero, one, two, ax: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; xi := xi.RoundToPrecision(dpcs); zero := '0'; seven := '0.7'; one := '1'; two := '2'; i := 1; fa := one; x := BigDecimal.Abs(xi); ax := x; if ax = one then begin // 1だったらπ/2 90° ans := pib / two; goto EXT; end; if ax > seven then begin // 0.7を超えたらcosの計算 x := one - x * x; x := x.RoundToPrecision(dpcs); x := BigDecimal.Sqrt(x, dpcs * 2); // sqrt は有効桁数が1/2 end; fb := two; ta := fa; tb := fb; ab := '3'; d := x * x * x; d := d.RoundToPrecision(dpcs); ans := (x + d * (ta / tb / ab)); repeat bans := ans; d := d * x * x; d := d.RoundToPrecision(dpcs); fa := fa + two; fb := fb + two; ab := ab + two; ta := ta * fa; ta := ta.RoundToPrecision(dpcs); tb := tb * fb; tb := tb.RoundToPrecision(dpcs); f := (ta / tb / ab) * d; ans := ans + f; inc(i); until (bans = ans) or (i > 10000); if ax > seven then ans := pib / two - ans; EXT: ans := ans.RoundToPrecision(pcsBack); if xi < zero then ans := -ans; result := ans; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arcsin(z) // カールソンの楕円積分 RCによる function arcsin_big(z: bigdecimal): bigdecimal; var zero, one, two, xi, eight, pib, sq: Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; eight := '0.8'; z := z.RoundToPrecision(dpcs); zero := '0'; one := '1'; two := '2'; pib := Pi_big / two; if BigDecimal.Abs(z) < eight then begin xi := one - z * z; xi := xi.RoundToPrecision(dpcs); result := z * Rc(xi, one); end else begin xi := z * z; xi := xi.RoundToPrecision(dpcs); result := RC(xi, one); sq := one - xi; sq := Bigdecimal.Sqrt(sq, dpcs * 2); result := sq * result; result := pib - result; if z < zero then result := -result; end; result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arccos(x) // 0.7を超えたら反対側の角度を求めます(arcsin計算になります) // arccos(x) = π/2 - arcsin(x) function arccosa_big(xi: bigdecimal): bigdecimal; label EXT; var i : integer; pib, ans, bans, fa, fb, d, ab, ta, tb, f, x, seven, zero, one, two, ax: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; xi := xi.RoundToPrecision(dpcs); zero := '0'; seven := '0.7'; two := '2'; one := '1'; i := 1; fa := '1'; x := BigDecimal.Abs(xi); ax := x; if ax = zero then begin // ゼロだったらπ/2 90° ans := pib / two; goto EXT; end; if ax > seven then begin // 0.7を超えたらcosの計算 x := one - x * x; x := x.RoundToPrecision(dpcs); x := BigDecimal.Sqrt(x, dpcs * 2); // sqrt は有効桁数が1/2 end; fb := '2'; ta := fa; tb := fb; ab := '3'; d := x * x * x; d := d.RoundToPrecision(dpcs); ans := (x + d * (ta / tb / ab)); repeat bans := ans; d := d * x * x; d := d.RoundToPrecision(dpcs); fa := fa + 2; fb := fb + 2; ab := ab + 2; ta := ta * fa; ta := ta.RoundToPrecision(dpcs); tb := tb * fb; tb := tb.RoundToPrecision(dpcs); f := (ta / tb / ab) * d; ans := ans + f; inc(i); until (bans = ans) or (i > 10000); if ax <= seven then ans := pib / two - ans; EXT: if xi < zero then ans := pib - ans; ans := ans.RoundToPrecision(pcsBack); result := ans; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arccos(z) // カールソンの楕円積分 RCによる function arccos_big(z: bigdecimal): bigdecimal; var v, t, e, one, two, zero, z1, pib: Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; z := z.RoundToPrecision(dpcs); zero := '0'; z1 := '0.1'; one := '1'; two := '2'; if bigdecimal.Abs(z) > z1 then begin v := z * z; v := v.RoundToPrecision(dpcs); t := Rc(v, one); e := BigDecimal.Sqrt(one - v, dpcs * 2); // sqrt は有効桁数が1/2 result := e * t; end else begin v := one - z * z; v := v.RoundToPrecision(dpcs); result := bigdecimal.Abs(z) * Rc(v, one); result := pib / two - result; end; if z < zero then result := pib - result; result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arctanf // 1 を越えたら反対側の角度を計算 // レオンハルト・オイラーの計算 function arctana_bigf(xi, pib: bigdecimal; dpcs: integer): bigdecimal; var i : integer; x, ax, y, ans, bans, f, ib, zero, one, two, a2: bigdecimal; begin xi := xi.RoundToPrecision(dpcs); x := bigdecimal.Abs(xi); ax := x; zero := '0'; one := '1'; two := '2'; if x > one then ax := one / x; a2 := ax * ax; a2 := a2.RoundToPrecision(dpcs); y := a2 / (one + a2); // y := y.RoundToScale(dpcs, rmNearesteven); i := 0; ib := one; ans := one; f := one; repeat bans := ans; f := f * ib * (two / (ib * two + one)) * y; f := f.RoundToPrecision(dpcs); ans := ans + f; ans := ans.RoundToPrecision(dpcs); ib := ib + one; inc(i); until (bans = ans) or (i > 10000); if ax <> zero then ans := ans * (y / ax) else ans := zero; if x > one then ans := pib / two - ans; if xi < zero then ans := -ans; result := ans; end; // arctan function arctana_big(xi: bigdecimal): bigdecimal; var pib: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; result := arctana_bigf(xi, pib, dpcs); // レオンハルト・オイラーの計算 result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arctanf(z) // カールソンの楕円積分 RCによる // 1 を越えたら反対側の角度を計算 function arctan_bigf(z, pib: bigdecimal; dpcs: integer): bigdecimal; var v, t, one, two, za: Bigdecimal; fr : integer; begin z := z.RoundToPrecision(dpcs); za := BigDecimal.Abs(z); fr := 0; one := '1'; two := '2'; if za > one then begin za := one / za; fr := 1; end; v := (one + za * za); v := v.RoundToPrecision(dpcs); t := Rc(one, v); result := za * t; if fr = 1 then result := pib / two - result; if z < 0 then result := -result; end; // arctan(z) function arctan_big(z: bigdecimal): bigdecimal; var pib: bigdecimal; pcsBack, dpcs: integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; result := arctan_bigf(z, pib, dpcs); // カールソンの楕円積分 RCによる // result := result.RoundToScale(pcsback, rmNearesteven); result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arctan2 // x=0,y=0 の場合定義なし0を返します。 // レオンハルト・オイラーの計算 function arctan2a_big(yi, xi: bigdecimal): bigdecimal; label EXT; var x, y, ans, zero, two, pib: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; xi := xi.RoundToPrecision(dpcs); yi := yi.RoundToPrecision(dpcs); zero := '0'; two := '2'; if xi = zero then begin if yi = zero then ans := zero; // 値無し 取り合えず 0 if yi > zero then ans := pib / two; if yi < zero then ans := -pib / two; goto EXT; end; if yi = zero then begin if xi > zero then ans := zero; if xi < zero then ans := pib; goto EXT; end; x := bigdecimal.Abs(xi); y := bigdecimal.Abs(yi); x := y / x; ans := arctana_bigf(x, pib, dpcs); // tan レオンハルト・オイラーの計算 if (yi > zero) and (xi < zero) then ans := pib - ans; if (yi < zero) and (xi < zero) then ans := ans - pib; if (yi < zero) and (xi > zero) then ans := -ans; EXT: // ans := ans.RoundToScale(pcsback, rmNearesteven); ans := ans.RoundToPrecision(pcsBack); result := ans; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arctan2 // カールソンのRCによる計算 // x=0,y=0 の場合定義なし0を返します。 function arctan2_big(yi, xi: bigdecimal): bigdecimal; label EXT; var x, y, ans, zero, two, pib: bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pib := pi_big; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; xi := xi.RoundToPrecision(dpcs); yi := yi.RoundToPrecision(dpcs); zero := '0'; two := '2'; if xi = zero then begin if yi = zero then ans := zero; // 値無し 取り合えず 0 if yi > zero then ans := pib / two; if yi < zero then ans := -pib / two; goto EXT; end; if yi = zero then begin if xi > zero then ans := zero; if xi < zero then ans := pib; goto EXT; end; x := bigdecimal.Abs(xi); y := bigdecimal.Abs(yi); x := y / x; ans := arctan_bigf(x, pib, dpcs); // カールソンのRCによる計算 if (yi > zero) and (xi < zero) then ans := pib - ans; if (yi < zero) and (xi < zero) then ans := ans - pib; if (yi < zero) and (xi > zero) then ans := -ans; EXT: // ans := ans.RoundToScale(pcsback, rmNearesteven); ans := ans.RoundToPrecision(pcsBack); result := ans; BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arsinh(z) // カールソンの楕円積分 RCによる function arsinh_big(z: bigdecimal): bigdecimal; var one, xi: Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; z := z.RoundToPrecision(dpcs); one := '1'; xi := one + z * z; xi := xi.RoundToPrecision(dpcs); result := z * Rc(xi, one); result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // arcosh(z) // カールソンの楕円積分 RCによる function arcosh_big(z: bigdecimal): bigdecimal; label EXT; var v, t, e, one: Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; z := z.RoundToPrecision(dpcs); one := '1'; if z <= one then begin application.MessageBox('１より大きくして下さい。','注意',0); result := '0'; goto EXT; end; v := z * z; v := v.RoundToPrecision(dpcs); t := Rc(v, one); e := BigDecimal.Sqrt(v - one, dpcs * 2); // sqrt は有効桁数が1/2 result := e * t; result := result.RoundToPrecision(pcsBack); EXT: BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // artanh(z) // カールソンの楕円積分 RCによる function artanh_big(z: bigdecimal): bigdecimal; label EXT; var v, t, one, two: Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; z := z.RoundToPrecision(dpcs); one := '1'; if z >= one then begin application.MessageBox('１より小さくして下さい。','注意',0); result := '0'; goto EXT; end; two := '2'; v := (one - z * z); v := v.RoundToPrecision(dpcs); t := Rc(one, v); result := z * t; result := result.RoundToPrecision(pcsBack); EXT: BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // ln(x) // カールソンの楕円積分 RCによる // あまり大きな値は不可 1e62が限度 // 小さな値は1e-26 function log_biga(xi: bigdecimal): bigdecimal; var v, t, e, one, two: Bigdecimal; dpcs : integer; begin dpcs := BigDecimal.DefaultPrecision; // 指定精度 xi := xi.RoundToPrecision(dpcs); one := '1'; two := '2'; v := (xi + one) / two; v := v * v; v := v.RoundToPrecision(dpcs); t := Rc(v, xi); e := xi - one; result := t * e; result := result.RoundToPrecision(dpcs); end; // x = m * 2^n のnの値と2^nの値を戻します // result = n // ta = 2^n // m = x / 2^n function frexpa(x: bigdecimal; var ta: biginteger): integer; var tb: biginteger; xi, one: bigdecimal; n, j, s, pcs : integer; begin pcs := BigDecimal.DefaultPrecision; one := '1'; one := one.RoundToPrecision(pcs); xi := x; if xi < one then x := one / x; ta := 1; n := 0; if x < ta then begin result := n; exit; end; j := 0; ta := '1'; s := 4096; // s=4^i i = 6 repeat while x > ta do begin biginteger.ShiftLeft(ta, s, tb); ta := tb; inc(j); end; if (s > 1) and (j > 0) then begin dec(j); biginteger.Shiftright(ta, s, tb); ta := tb; end; n := n + s * j; j := 0; s := s div 4; until s = 0; if xi < one then begin n := -n + 1; biginteger.Shiftright(ta, 1, tb); ta := tb; end; result := n // ta = 2^n end; // loge(x) ln(x) // 級数展開 // 1e±500000 程度が限度 // FL True 有効桁数補正あり False 有効桁数補正無 function log_bigf(xi: bigdecimal; FL: boolean): bigdecimal; label EXT; const KK = 10000; var LOG2, SQRT2, two, last, s, one, one1, x2, i, k2: bigdecimal; k : integer; kb, ta : biginteger; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin result := '0'; if xi <= result then begin application.MessageBox('無効な値です。','注意',0); exit; end; s := '1e-5'; last := '1e10'; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack; one := '1'; if xi = one then begin // result := '0'; // 設定済み goto EXT; end; if (xi > s) and (xi < last) then begin // 指数の値が小さい場合はRC計算 if FL then begin dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; end; result := log_biga(xi); // RC計算 goto EXT; end; one1 := '1e500000'; if xi < one then begin last := one / xi; if last > one1 then begin application.MessageBox('値が小さすぎます。','注意',0); goto EXT; end; end; if xi > one1 then begin application.MessageBox('値が大きすぎます。','注意',0); goto EXT; end; if FL then begin dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; end; one1 := one.RoundToPrecision(dpcs); two := '2'; LOG2 := log_biga(two); // log(2) SQRT2 := bigdecimal.Sqrt(two, dpcs * 2); k := frexpa(xi / SQRT2, ta); // x / √2 = m * 2^k ta = 2^k kb := k; if k < 0 then begin // kが負数だったら // k := -k; // 符号反転 k2 := ta; // 2^k intege to decimal k2 := k2.RoundToPrecision(dpcs); k2 := one1 / k2; // 1/(2^k) = 2^-k end else k2 := ta; // 2^k intege to decimal k2 := k2.RoundToPrecision(dpcs); xi := xi / k2; // x / 2^k // Main.Form1.Memo1.Lines.Append('x=' + xi.ToPlainString); xi := (xi - one) / (xi + one); x2 := xi * xi; x2 := x2.RoundToPrecision(dpcs); i := one; s := xi; k := 0; // 長ループ脱出用 repeat xi := xi * x2; xi := xi.RoundToPrecision(dpcs); i := i + two; last := s; s := s + xi / i; s := s.RoundToPrecision(dpcs); inc(k); until (last = s) or (k > KK); result := LOG2 * kb + two * s; if k > KK then begin application.MessageBox('演算エラー。','注意',0); result := '0'; end; EXT: result := result.RoundToPrecision(pcsBack); if FL then begin BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; end; // log(x) function log_big(xi: bigdecimal): bigdecimal; begin result := log_bigf(xi, True); end; // sub e^x function expb(x: Bigdecimal; k: bigInteger): Bigdecimal; var w : bigdecimal; izero, ione, km : bigInteger; dpcs : integer; begin dpcs := BigDecimal.DefaultPrecision; izero := '0'; ione := '1'; if k >= izero then w := '2' else begin w := '0.5'; k := -k; end; while k <> izero do begin if k and ione = ione then begin x := x * w; x := x.RoundToPrecision(dpcs); end; w := w * w; w := w.RoundToPrecision(dpcs); bigInteger.ShiftRight(k, 1, km); k := km; end; result := x; end; // e^x function expa(x: bigdecimal): bigdecimal; var neg, dpcs : integer; km, a, e, prev, LOG2, zero, one, two : bigdecimal; k, i, ione: bigInteger; begin dpcs := BigDecimal.DefaultPrecision; zero := '0'; one := '1'; two := '2'; LOG2 := log_big(two); // log(2) km := (x / LOG2); km := km.Floor; k := km.AsBigInteger; x := x - k * LOG2; x := x.RoundToPrecision(dpcs); if x >= zero then neg := 0 // 負数フラグ 0 else begin neg := 1; // 負数フラグ 1 x := -x; end; e := one + x; a := x; i := '2'; ione := '1'; repeat prev := e; a := a * (x / i); a := a.RoundToPrecision(dpcs); e := e + a; e := e.RoundToPrecision(dpcs); i := i + ione; until e = prev; if neg <> 0 then e := one / e; // 負数なら逆数 result := expb(e, k); end; // e^x function exp_big(x: bigDecimal): bigdecimal; var d : bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin d := '100000000'; if x > d then begin result := '0'; application.MessageBox('値が大きすぎます。','注意',0); exit; end; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; x := x.RoundToPrecision(dpcs); result := expa(x); result := result.RoundToPrecision(pcsBack); BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // power(x, y) x^y e^(log(x) * y) function power_big(xi, yi: bigdecimal): bigdecimal; label EXT; var t, lx : Bigdecimal; d : Bigdecimal; pcsBack, dpcs : integer; rmback : BigDecimal.RoundingMode; begin result := '0'; if xi < result then begin application.MessageBox('無効な値です。','注意',0); exit; end; pcsBack := BigDecimal.DefaultPrecision; // 除算演算精度バックアップ rmBack := BigDecimal.DefaultRoundingMode; // 除算丸めモードバックアップ dpcs := pcsBack + 5; // 指定精度 + α BigDecimal.DefaultPrecision := dpcs; BigDecimal.DefaultRoundingMode := rmNearesteven; xi := xi.RoundToPrecision(dpcs); yi := yi.RoundToPrecision(dpcs); // log(x) t := log_bigf(xi, False); lx := t * yi; lx := lx.RoundToPrecision(dpcs); d := '100000000'; if BigDecimal.abs(lx) > d then begin result := '0'; application.MessageBox('値が大きすぎます。','注意',0); goto EXT; end; // e^x result := expa(lx); result := result.RoundToPrecision(pcsBack); EXT: BigDecimal.DefaultPrecision := pcsBack; // 除算演算精度復帰 BigDecimal.DefaultRoundingMode := rmBack; // 除算丸めモード復帰 end; // Bigdecimal -> Double function BigdecimalToDouble(b: bigdecimal): double; const maxdouble = 1.7976931348623158E308; // 1.7976931348623'2'E308 -> INF mindouble = 4.94065645841247E-324; var bigstr : string; absb : bigdecimal; begin absb := bigdecimal.Abs(b); if absb = bigdecimal.Zero then begin result := 0; exit; end; if absb >= maxdouble then begin result := maxdouble; if b < bigdecimal.Zero then result := -result; exit; end; if absb <= mindouble then begin result := mindouble; if b < bigdecimal.Zero then result := -result; exit; end; b := b.RoundToPrecision(30); bigstr := b.ToString; result := strtofloat(bigstr); end; end.
Big_Math.zip
三角関数、逆三角関数、その他関数 に戻る