マクローリン展開によるsin(x)の計算です。角度の範囲は、±90°になっているので、360°の角度で計算する場合は、別途角度変換が必要です。
計算は、Double、多倍長、BigDecimal、多倍長sin->BigDecimal変換でプログラムを作成しました。
計算の精度指定桁数はDoubleを除き10000桁までです。
通常の計算では、そこまで必要はないかもしれません。
BigDecimalによるSin(x)の計算
BigDecimalによるπの計算プログラムを作成したので、Sin(x)の計算についても検討してみました。
sin(x)の中で使用されるπの値は、BigDecimalによるπの計算にあるBigDecimalのπ計算ルーチンを利用します。
Double程度の有効桁数であれば、固定値で良いのですが、有効桁数が非常に大きい場合を考慮して、計算により求めます。
問題無く計算出来れば、BigDecimal用の各種関数サブルーチン作成に使用します。
BigDecimalの組み込みについては、ベルヌーイ数 その４を参照してください。
マクローリン展開によるsin(x)の計算です。
角度の範囲は、±90°になっているので、360°の角度で計算する場合は、別途角度変換が必要です。
計算は、Double、多倍長、BigDecimal、多倍長sin->BigDecimal変換でプログラムを作成しました。
計算の精度指定桁数はDoubleを除き10000桁までです。
通常の計算では、そこまで必要はないかもしれません。
プログラム
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)
BitBtn1: TBitBtn;
LabeledEdit1: TLabeledEdit;
Memo1: TMemo;
BitBtn2: TBitBtn;
BitBtn3: TBitBtn;
BitBtn4: TBitBtn;
BitBtn5: TBitBtn;
LabeledEdit2: TLabeledEdit;
CheckBox1: TCheckBox;
procedure BitBtn1Click(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
procedure BitBtn3Click(Sender: TObject);
procedure BitBtn4Click(Sender: TObject);
procedure BitBtn5Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ Private 宣言 }
function inputcheck(DF: Boolean): boolean;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses mp_real, mp_types, mp_base, Velthuis.BigDecimals;
var
DRD : integer; // 有効桁数
dsr : string; // πの値string
const
DEG = '180'; // 180°
// πの計算 ガウス＝ルジャンドルのアルゴリズム
// pcs πの精度
// np true ループ数 表示
function pi_big(pcs: integer; np: boolean): 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; // 指定精度 + α
BigDecimal.DefaultPrecision := dpcs;
BigDecimal.DefaultRoundingMode := rmNearesteven;
one := BigDecimal.One;
two := BigDecimal.Two;
EPSILON := one / BigDecimal.IntPower(two, pcs shl 1); // 収束判定値
one := one.RoundToPrecision(dpcs); // oneの有効桁数をdpcsに設定
two := two.RoundToPrecision(dpcs); // twoの有効桁数をdpcsに設定
four := two + two;
a := one;
b := one / BigDecimal.Sqrt(two);
t := one / four; // 除算精度
x := '1';
n := 0;
repeat
an := (a + b) / two; // 平均値 除算有効桁数での丸め
b := a * b;
b := b.RoundToPrecision(dpcs); // pcs + α 丸め
b := BigDecimal.Sqrt(b);
t := t - x * (an - a) * (an - a);
t := t.RoundToPrecision(dpcs); // pcs + α 丸め
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; // 除算丸めモード復帰
if np then
Form1.Memo1.Lines.Append('πの計算 ガウス＝ルジャンドル Big n = ' + intTostr(n));
end;
// 入力チェック
// 有効桁数の読み込みセット DRD
// πの値をmp_float又はBigDecimalからString変換 dsr。
// DF double フラグ
function TForm1.inputcheck(DF: Boolean): boolean;
var
x : double;
c : integer;
paif : mp_float;
paib : BigDecimal;
begin
result := false;
val(labelededit1.Text, x, c);
if c <> 0 then begin
memo1.Text := 'xの値に間違いがあります。';
exit;
end;
if abs(x) > 90 then begin
memo1.Text := 'xの値は90以下にして下さい。';
exit;
end;
val(labelededit2.Text, DRD, c);
if c <> 0 then begin
memo1.Text := 'xの値に間違いがあります。';
exit;
end;
if (DRD < 30) or (DRD > 1000) then begin
memo1.Text := '精度桁数の値は30以上 又は 1000迄にして下さい。';
exit;
end;
if not DF then
if checkbox1.Checked = False then begin // πの値計算選択
paib := pi_big(DRD, true); // 有効桁数 (DRD + 10)
dsr := paib.ToPlainString;
end
else begin
mpf_set_default_decprec(DRD + 10); // 有効桁数 DRD + 10に設定
mpf_init(paif);
mpf_set_pi(paif); // πの値
dsr := string(mpf_adecimal_alt(paif, DRD + 10)); // 有効桁数 DRD + 10 String変換
mpf_clear(paif);
Form1.Memo1.Lines.Append('πの値 多倍長');
end;
result := true;
end;
// sin(x)計算 double
// マクローリン展開
// delphiには標準でsin(x)が組み込まれています参考用です。
function sindouble(xi: double): double;
var
x, v, d, sq, u : double;
k : integer;
begin
x := abs(xi);
d := x;
sq := x * x;
u := d;
k := 0;
repeat
k := k + 2;
v := u;
d := -d * sq / (k * (k + 1));
u := v + d;
until (v = u) or (k > 100);
if xi < 0 then v := -v;
result := v;
Form1.memo1.Lines.Append('Double sin 変換ループ n=' + intTostr(k)); // Loop表示
end;
procedure TForm1.BitBtn1Click(Sender: TObject);
var
x : double;
begin
if not inputcheck(True) then exit;;
x := strTofloat(labelededit1.Text);
x := x / 180 * pi;
x := sindouble(x);
memo1.Lines.Append(floatTostr(x));
memo1.Lines.Append('');
end;
// sin(x)計算 多倍長
// マクローリン展開
// mp_arith には、sin関数がありここのプログラムは参考用です。
procedure sinmpf(var xi, v : mp_float);
var
tmp0, tmp1 : mp_float;
x, K: mp_float;
d, sq, u, one, two, vd : mp_float;
c : integer;
begin
mpf_init2(tmp0, tmp1);
mpf_init3(x, K, d);
mpf_init5(sq, u, one, two, vd);
mpf_abs(xi, x); // 絶対値
mpf_set1(one); // 1
mpf_add(one, one, two); // 2
mpf_copy(x, d); // x = d
mpf_sqr(x, sq); // sq = x * x
mpf_copy(x, u); // x = u
mpf_set0(k); // k = 0
c := 0;
repeat // repeat
mpf_add(k, two, k); // k = k + 2
mpf_copy(u, v); // u = v
mpf_add(k, one, tmp0); // tmp0 = k + 1
mpf_mul(tmp0, k, tmp1); // tmp1 = k(k + 1)
mpf_div(sq, tmp1, tmp0); // tmp0 = sq / k(k + 1)
mpf_chs(d, d); // d = -d
mpf_mul(d, tmp0, d); // d = d * tmp0 d = -d * sq / k(k + 1)
mpf_add(v, d, u); // u := v + d
inc(c); // c = c + 1
until mpf_is_eq(v, u) or (c > 500); // until (v = u) or (c > 200)
if not mp_real.s_mpf_is_ge0(xi) then mpf_chs(v, v); // 符号設定
Form1.memo1.Lines.Append('Mp_float sin 変換ループ n=' + intTostr(c)); // Loop表示
mpf_clear2(tmp0, tmp1);
mpf_clear3(x, K, d);
mpf_clear5(sq, u, one, two, vd);
end;
procedure TForm1.BitBtn2Click(Sender: TObject);
var
x, v, pai : mp_float;
begin
if not inputcheck(False) then exit;;
mpf_set_default_decprec(DRD + 10); // 精度DRD + 10桁に設定
mpf_init3(x, v, pai);
if checkbox1.Checked = False then
mpf_read_decimal(pai, pansichar(ansistring(dsr))) // π BigDecimal計算値
else
mpf_set_pi(pai); // π
mpf_read_decimal(x, pansichar(ansistring(Form1.LabeledEdit1.Text))); // 入力
mpf_div_int(x, 180, x); // Rad変換 x = x / 180
mpf_mul(x, pai, x); // x = x * pi
sinmpf(x, v); // sin(x)
memo1.Lines.Append(string(mpf_adecimal_alt(v, DRD))); // DRD桁表示
memo1.Lines.Append('');
mpf_clear3(x, v, pai);
end;
// sin(x)計算 BigDecimal
// マクローリン展開
// BigDecimalは、
// 除算時、有効桁数がDefaultより大きくて、分子の有効桁数が分母の有効桁数より大きい場合
// 除算に先だって分子に小数部がある時はその小数点以下桁数を、RoundToで小さくするか、
// Defaultの有効桁数を大きくする必要があります、多少は余裕があるようですが、注意が必要です。
// 整数部の桁数が有効桁数より大きい場合、有効桁数以外には"0"がはいります。
// 演算での有効桁数丸めは除算でしか発生しません。
function sinBig(xi : BigDecimal): BigDecimal;
var
x, v, d, sq, u, k : BigDecimal;
c: integer;
begin
x := BigDecimal.Abs(xi); // 絶対値
d := x;
sq := x * x; // sqは小数点以下 (DRD + 10) * 2 桁
u := d;
k := 0;
c := 0;
repeat
k := k + 2;
v := u;
d := -d * sq; // 小数点以下(DRD + 10)×2 桁
d := d.RoundToPrecision(DRD + 10); // 有効桁数DRD + 10桁に丸めないと次の除算不可分母は小数点以下無し整数
d := d / (k * (k + 1)); // 除算で有効桁数DRD + 10桁 除算Default値
u := v + d;
u := u.RoundToPrecision(DRD + 10); // 有効桁数DRD + 10桁に丸めないと次の除算不可分母は小数点以下無し整数
inc(c);
until (v = u) or (c > 500);
if xi < 0 then v := -v; // 符号設定
result := v;
Form1.memo1.Lines.Append('BigDecimal sin 変換ループ n=' + intTostr(c)); // Loop表示
end;
procedure TForm1.BitBtn3Click(Sender: TObject);
var
x , v : BigDecimal;
pai, d180 : BigDecimal;
begin
if not inputcheck(False) then exit;;
BigDecimal.DefaultPrecision := DRD + 10; // 除算でのデフォルトの有効桁数 整数部が有効桁数を超える場合有効桁数以外は"0"
BigDecimal.DefaultRoundingMode := rmNearestDown; // デフォルトの丸め指定
if checkbox1.Checked = False then // πの値設定
pai := pi_big(DRD + 10, False) // πの計算 ガウス＝ルジャンドルのアルゴリズム
else
pai := dsr; // πの値設定 pai := dsr 多倍長->big
x := labelededit1.Text; // 角度入力 誤差を防ぐ為string直接変換
x := x * pai; // radに変換 xは小数点以下100 + x 小数点以下桁数入力値次第
x := x.RoundTo(-DRD - 2); // xを小数点以下DRD + 2桁に丸めないと次の除算出来ない時があります
d180 := DEG;
x := x / d180; // xは有効桁数DRD + 10桁 除算Default値 x < pi / 2
v := sinBig(x); // sin(x)
v := v.RoundTo(-DRD + 1); // 小数点以下DRD-1桁表示DRD桁で丸め |V| <=1
v := BigDecimal(v).RemoveTrailingZeros(0); // 末尾のゼロ消去
memo1.Lines.Append(v.ToPlainString); // 値表示
memo1.Lines.Append('');
end;
// 多倍長浮動小数点の利用によるsin(x)の計算
// mp_floatとBigDecimalの値はStringで変換
function sinmpg(xi : BigDecimal): BigDecimal;
var
x, v : mp_float;
begin
mpf_set_default_decprec(DRD + 10); // mp_float 精度DRD + 10桁に設定
mpf_init2(x, v);
mpf_read_decimal(x, pansichar(ansistring(xi.ToPlainString))); // BigDecimal -> mp_float
mpf_sin(x, v); // 多倍長sin(x)
result := string(mpf_adecimal_alt(v, DRD)); // mp_float -> BigDecimal
Form1.memo1.Lines.Append('多倍長->BIG sin 変換ループ無し'); // Loop表示
mpf_clear2(x, v);
end;
procedure TForm1.BitBtn4Click(Sender: TObject);
var
x , v : BigDecimal;
pai, d180 : BigDecimal;
begin
if not inputcheck(False) then exit;;
BigDecimal.DefaultPrecision := DRD + 10; // 除算でのデフォルトの有効桁数 有効桁数を超える場合は先頭からの桁数
BigDecimal.DefaultRoundingMode := rmNearesteven; // デフォルトの丸め指定
pai := dsr; // πの値設定
x := labelededit1.Text; // 角度入力 誤差を防ぐ為string直接変換
x := x * pai; // radに変換 小数点以下桁数入力値次第
x := x.RoundToPrecision(DRD + 10); // xを有効桁数DRD + 10桁に丸めないと次の除算出来ない時があります
d180 := DEG;
d180 := d180.RoundToPrecision(DRD + 10);
x := x / d180; // xは有効桁数DRD + 10桁 除算Default値 x < pi / 2
v := sinmpg(x); // sin(x)
v := v.RoundTo(-DRD + 1); // 小数点以下DRD桁で丸め
v := BigDecimal(v).RemoveTrailingZeros(0); // 末尾のゼロ消去
memo1.Lines.Append(v.ToPlainString); // 値表示
memo1.Lines.Append('');
end;
procedure TForm1.BitBtn5Click(Sender: TObject);
begin
memo1.Clear;
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
memo1.Width := 820; // 一行100桁表示に設定
memo1.Font.Size := 11; // 一行100桁表示にフォントサイズ設定
end;
end.
sin_x.zip
三角関数、逆三角関数、その他関数 に戻る