シュトルーベ関数

 Struve 関数は天文学種Karl Hermann Struve(1854-1920)が光の輝線スペクトルの強度を記述する関数として導入したものの様ですが、輝線との関係を示す部分は見けることが出来ませんでした。
計算は、ベッセル関数と同じで係数が少し違うだけなりで、プログラムは、ベッセル関数のプログラムを元に修正しています。

 Struve 関数は、虚数部の計算は不要なので、虚数部はゼロ(0)として、計算します。
虚数部が必用で有れば、虚数部無しのチェックボックスのチェックを外せば、複素数として虚数部の入力が可能となります。



// シュトルーベ関数の計算プログラムをベッセル関数の計算プログラムを基に作成したので
// 複素数の計算となっています。
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,
  VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, VCLTee.TeeProcs,
  VCLTee.Chart, System.Diagnostics;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    LabeledEdit1: TLabeledEdit;
    LabeledEdit2: TLabeledEdit;
    Chart1: TChart;
    Series1: TLineSeries;
    Series2: TLineSeries;
    Series3: TPointSeries;
    Series4: TPointSeries;
    Series5: TLineSeries;
    Series6: TLineSeries;
    CheckBox2: TCheckBox;
    LabeledEdit3: TLabeledEdit;
    CheckBox3: TCheckBox;
    NonComplex: TCheckBox;
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
    procedure NonComplexClick(Sender: TObject);
  private
    { Private 宣言 }
    procedure chart1coment;
  public
    { Public 宣言 }
  end;

var
  Form1: TForm1;

implementation

uses system.Math, Velthuis.BigIntegers,mp_real, mp_cmplx, mp_types, mp_base;

{$R *.dfm}

type
  DArray = array of double;

var
  BM : array of mp_float;       // ベルヌーイ数配列
  FA : array of mp_float;       // Γ(k+3/2)配列
  PVG : array of mp_float;      // +VΓ

  zero, one, two, four : mp_float;
  three, pai, log_2pis2 : mp_float;
  oneahf : mp_float;

const
  KMmax = 250;                  // KM max   250
  Vmax  = 100;                  // v 次数 max

//------------------------------------------------------------------------------
  NB = 120;                            // ベルヌーイ数  配列数 NB + 1

var
  NumeratorString : array[0..NB] of string;               // 分子
  DenominatorString  : array[0..NB] of string;            // 分母

// 最大公約数  ユークリッドの互除法 BigInteger
function gcd_BigInteger(x, y: BigInteger): BigInteger;
var
  t : BigInteger;
begin
  while y <> 0 do begin
    t := x mod y;
    x := y;
    y := t;
  end;
  result := x;
end;

// ベルヌーイ数
// Akiyama-Tanigawa algorithm
// BigInteger
procedure Bernoulli_number_BigInteger;
const
  n = (NB + 1) * 2;
var
  m, j, k : integer;
  a     : array of BigInteger;      // 分子
  b     : array of BigInteger;      // 分母
  tmpN  : BigInteger;               // 分子
  tmpD  : BigInteger;               // 分母
  gcd   : BigInteger;               // 最大公約数
  b0    : BigInteger;
begin
  setlength(a, n + 1);
  setlength(b, n + 1);
  k := 0;
  for m := 0 to n do begin
    a[m] := 1;                         // 分子 初期値
    b[m] := (m + 1);                                // 分母 初期値
    for j := m - 1 downto 0 do begin
      tmpN := a[j] * b[j + 1] - a[j + 1] * b[j];    // 分子
      tmpN := tmpN * (j + 1);                       // 〃
      tmpD := b[j] * b[j + 1];                      // 分母
      gcd := gcd_BigInteger(tmpN, tmpD);            // 最大公約数
      a[j] := tmpN div gcd;
      b[j] := tmpD div gcd;
    end;
    if (m > 0) and (m mod 2 = 0) then begin
      b0 := b[0];
      b0 := b0 * m * (m -1);                       // m ベルヌーイ数No
      NumeratorString[k] := a[0].tostring;
      DenominatorString[k] := b0.tostring;
      inc(k);
    end;
  end;
end;

//------------------------------------------------------------------------------

// ログガンマ多倍長
procedure log_GammaMul(var x, ans : mp_float);
var
  v, w : mp_float;
  tmp, tmp0, s : mp_float;
  i : integer;
begin
  mpf_init2(v, w);
  mpf_init3(tmp, tmp0, s);

  mpf_set1(v);
  mpf_set_int(tmp, NB);
  while mpf_is_lt(x, tmp) do begin
    mpf_mul(v, x, v);
    mpf_add(x, one, x);
  end;
  mpf_mul(x, x, tmp);               // x^2
  mpf_div(one, tmp, w);             // w = 1 / x^2
  mpf_set0(s);
  for i := NB downto 1 do begin
    mpf_add(s, BM[i], tmp);         // tmp = s + B[i]
    mpf_mul(tmp, w, s);             // s = tmp * w
  end;
  mpf_add(s, BM[0], tmp);           // tmp = s + B[0]
  mpf_div(tmp, x, s);               // s = (s + B[0]) / x
  mpf_add(s, log_2pis2, s);         // s = s + ln(2π)/2
  mpf_ln(v, tmp);                   // ln(v)
  mpf_sub(s, tmp, s);               // s := s - ln(v)
  mpf_sub(s, x, s);                 // s := s - x
  mpf_div(one, two, tmp);           // tmp = 1/2
  mpf_sub(x, tmp, tmp0);            // tmp0 = x - 1/2
  mpf_ln(x, tmp);                   // ln(x)
  mpf_mul(tmp0, tmp, tmp0);         // tmp0 = (x - 1/2) * ln(x)
  mpf_add(s, tmp0, ans);            // ans = s + (x - 1/2) * ln(x)

  mpf_clear2(v, w);
  mpf_clear3(tmp, tmp0, s);
end;

// 多倍長ガンマ
// xの値が 0 と負整数の時Γは∞になるのですが処理をしていませんのでエラーになります。
// ケルビンの次数が整数の時は使用されません。
procedure gammaMul(var x, ans: mp_float);
var
  tmp, tmp0, logG : mp_float;
begin
  mpf_init3(tmp, tmp0, logG);

  if mpf_is_lt(x, zero) then begin
    mpf_mul(pai, x, tmp);            // x*π
    mpf_sin(tmp, tmp0);              // sin(πx);
    mpf_sub(one, x, tmp);            // 1-x
    log_GammaMul(tmp, logG);         // loggamma(1-x);
    mpf_exp(logG, logG);             // exp(logG)
    mpf_div(pai, tmp0, tmp);         // π/sin(πx)
    mpf_div(tmp, logG, ans);         // ans = π/(sin(πx) * logG(1-x))
  end
  else begin
    log_GammaMul(x, logG);           // logG
    mpf_exp(logG, ans);              // exp(logG)
  end;

  mpf_clear3(tmp, tmp0, logG);
end;


//------------------------------- complex power --------------------------------
// mpc_powは**.5時正しい答えを返しません
// 本来ゼロになる所にゼロで無い値が現れます
// x(-0.5, 0)   y(1.5, 0) 有効桁数70前後
// x(-0.5, -0.5) y(4, 0)  有効桁数60前後
// x(-0.5, 0.5)  y(2, 0) 有効桁数70前後
procedure mpc_powa(x, y : mp_complex; var ans : mp_complex);
var
  ay, ai, zero, harf, two : mp_float;
begin
  mpf_init5(ay, ai, zero, harf, two);

  mpf_set_int(two, 2);                  // 2
  mpc_pow(x, y, ans);                     // exp(y*ln(x))
  // xの虚数部とyの虚数部0なら
  if mpf_is0(x.im) and mpf_is0(y.im) then begin
    mpf_set0(zero);                       // 0
    mpf_inv(two, harf);                   // 1/2  0.5
    mpf_abs(y.re, ay);                    // yの実数部の絶対値
    mpf_frac(ay, ai);                     // yの実数部の絶対値の小数部
    mpf_sub(ai, harf, ai);                // yの実数部の絶対値の小数部と0.5の差分
    // 差分がゼロ(yの実数部の小数部が0.5)
    if mpf_is0(ai) then
      // xの実数部が0と等しいか0より大きいなら
      if mpf_is_ge(x.re, zero) then mpf_set0(ans.im)   // 答えの虚数部0
                               // 0より小さいなら
                               else mpf_set0(ans.re);  // 答えの実数部0
    mpf_frac(y.re, ai);                                // 乗数の小数部
    if mpf_is0(ai) then mpf_set0(ans.im);              // 整数なら答えの虚数部0
  end;
  // xの整数部とyの虚数部が0なら
  if mpf_is0(x.re) and mpf_is0(y.im) then begin
    mpf_frac(y.re, ai);                                // 乗数の小数部
    if mpf_is0(ai) then begin                          // 整数なら
      mpf_div(y.re, two, ay);                          // 偶数奇数確認
      mpf_frac(ay, ai);                                // 小数部
      if mpf_is0(ai) then mpf_set0(ans.im)             // 偶数なら 虚数部0
                     else mpf_set0(ans.re)             // 奇数なら実数部0
    end;
  end;
  mpf_abs(x.re, ay);
  mpf_abs(x.im, ai);
  // xの実数部と虚数部の絶対値が等しくてyがゼロでなかったら
  if mpf_is_eq(ay, ai) and mpf_is0(y.im) and not mpc_is0(y) then begin
    mpf_frac(y.re, ai);                                // 乗数の小数部
    if mpf_is0(ai) then begin                          // 整数なら
      mpf_div(y.re, two, ay);                          // 偶数奇数確認
      mpf_frac(ay, ai);                                // 小数部
      if mpf_is0(ai) then begin
        mpf_add(two, two, two);                        // 4
        mpf_div(y.re, two, ay);                        // 4偶数奇数確認
        mpf_frac(ay, ai);                              // 小数部
        if mpf_is0(ai) then mpf_set0(ans.im)           // 偶数なら 虚数部0
                       else mpf_set0(ans.re);          // 奇数なら 実数部0
      end;
    end;
  end;

  mpf_clear5(ay, ai, zero, harf, two);
end;

//------------------------------------------------------------------------------
// Hv(x) Lhv 多倍長
// X 値 複素数
// v 次数
procedure Hvx(v : mp_float; var x, ri: mp_complex);
var
  k : integer;
  s, x24k, tmp, tmp0, tmp1 : mp_complex;
  xc, xs2: mp_complex;
  kd, nd, vk : mp_float;
  khg, kf : mp_float;
  tmf, tmf0 : mp_float;
  sb : mp_complex;
begin
  mpc_init5(s, x24k, tmp, tmp0, tmp1);
  mpc_init2(xc, xs2);
  mpf_init2(khg, kf);
  mpf_init3(kd, nd, vk);
  mpf_init2(tmf, tmf0);
  mpc_init(sb);

  mpc_set0(s);                              // Σ=0
  mpc_set0(sb);
  mpc_set1(x24k);                           // k=0 ((x^2)/4)^0 初期値

  mpc_div_mpf(x, two, xs2);                 // x / 2

  mpc_mul(xs2, xs2, xc);                    // (x^2)/ 4

  for k := 0 to KMmax do begin
    mpf_set1(nd);
    mpf_add_int(v, k, tmf0);                // k + v
    mpf_add(tmf0, oneahf, vk);              // vk = k + v + 3/2;
    if mpf_is_lt(vk, zero) then begin       // vk < 0 時 nxが整数か確認
      mpf_int(vk, tmf);                     // int(vk);
      mpf_sub(vk, tmf, nd);                 // vk - int(vk) vkが負の整数だったら vk = 0
    end;
    if mpf_is_ne(vk, zero) and mpf_is_ne(nd, zero) then begin  // vkが負の整数の時は計算しない
      mpf_mul(FA[k], PVG[k], khg);          // Γ(k+3/2)Γ(v+k+3/2) vkが0,負の整数の時±∞
      mpc_div_mpf(x24k, khg, tmp0);         // ((i(x^2)/4)^k)/(Γ(k+3/2)Γ(v+k+3/2))
      if Form1.CheckBox3.Checked = false then       // true  Lvx
        if k mod 2 <> 0 then mpc_chs(tmp0, tmp0);   // false Hvx (-1)^k
      mpc_add(s, tmp0, s);                  // Σ
      mpc_sub(s, sb, tmp0);
      if mpc_is0(tmp0) then break;
      mpc_copy(s, sb);
    end;
    mpc_mul(x24k, xc, x24k);                // (i(x^2)/4)^k
  end;
  mpf_add(v, one, tmf0);                    // V+1
  mpc_set_mpf(tmp, tmf0, zero);             // (v+1)+ 0i
  // x が0で次数vが負数の時power演算ゼロでの除算防止
  if mpc_is0(x) and mpf_is_lt(tmf0, zero) then // (V+1)<0 x=0 時は無限大になるので計算しない
  else begin
    mpc_powa(xs2, tmp, tmp0);               // (x/2)^(v+1)
    mpc_mul(s, tmp0, s);                    // (x/2)^(v+1) * Σ
  end;
  mpc_copy(s, ri);

  mpc_clear5(s, x24k, tmp, tmp0, tmp1);
  mpc_clear2(xc, xs2);
  mpf_clear2(khg, kf);
  mpf_clear3(kd, nd, vk);
  mpf_clear2(tmf, tmf0);
  mpc_clear(sb);
end;

// グラフバックグラウンド消去
procedure ChartBackclear;
begin
  with Form1.Chart1.BackImage.Bitmap.Canvas do begin
    Brush.Color := clBtnFace;
    Brush.Style := bsSolid;
    FillRect(rect(0, 0, Form1.Chart1.Width, Form1.Chart1.Height));
    Brush.Style := bsClear;
  end;
end;

// グラフバックグラウンド ラインの種類表示
procedure TForm1.chart1coment;
begin
  with Chart1.BackImage.Bitmap.Canvas do begin
    Font.Size := 8;
    Font.Style := [];
    Font.Color := clred;
    TextOut(chart1.Width - 30, 1,'実数');
    Font.Color := clblue;
    if NonComplex.Checked = false then
      TextOut(chart1.Width - 30, 13,'虚数');
    Pen.Width := 1;
    Pen.Color := clred;
    MoveTo(chart1.Width - 50, 6);
    LineTo(chart1.Width - 32, 6);
    Pen.Color := clblue;
    if NonComplex.Checked = false then begin
      MoveTo(chart1.Width - 50, 18);
      LineTo(chart1.Width - 32, 18);
    end;
  end;
end;

// 計算
// xの値が大きくなると誤差が大きくなります。
procedure TForm1.Button1Click(Sender: TObject);
label
  EXT;
const
  N = 200;
  x0m = 1e-50;        // ゼロ近傍値 infinty の符号設定計算値
  zeromg = '1E-90';   // 計算結果表示これより小さい値は0にします 演算桁数の1/2
var
  ch, i : integer;
  berxe, beixe: double;
  xin, vin, xv, rk, ik: double;
  xmin, xmax, dx : double;
  ymin, ymax, yre, yim : double;
  xm, vm, xvm, ixm: mp_float;
  nd, tmf : mp_float;
  ri : mp_complex;
  xc, xcb : mp_complex;
  k : integer;
  vk, avm, tmf0, zz : mp_float;

  redata : DArray;
  imdata : DArray;

  ix, x : double;
  istr : string;

  // double to mpc  グラフ計算用
  procedure xvtoXc(xv, ia : double; var xc : mp_complex);
  begin
    mpf_set_dbl(xvm, xv);
    mpf_set_dbl(tmf, ia);
    mpc_set_mpf(xc, xvm, tmf);
  end;

begin
  memo1.Clear;
  val(labelededit1.Text, vin, ch);
  if ch <> 0 then begin
    application.MessageBox('次数vの値に間違いがあります。','注意', 0);
    exit;
  end;
  if abs(vin) > Vmax then begin
    application.MessageBox('次数vの値が計算範囲外です。','注意', 0);
    exit;
  end;
  val(labelededit2.Text, xin, ch);
  if ch <> 0 then begin
    application.MessageBox('xの値に間違いがあります。','注意', 0);
    exit;
  end;
  if abs(xin) > 100 then begin
    application.MessageBox('xが100を越えると条件によって誤差が大きくなります。','注意', 0);
  end;
  if NonComplex.Checked = false then begin
    val(labelededit3.Text, ix, ch);
    if ch <> 0 then begin
      application.MessageBox('i の値に間違いがあります。','注意', 0);
      exit;
    end;
    if abs(ix) > 100 then begin
      application.MessageBox('iが100を越えると条件によって誤差が大きくなります。','注意', 0);
    end;
  end
  else ix := 0;

  mpf_init4(xm, vm, xvm, ixm);
  mpf_init3(nd, tmf, zz);
  mpf_init3(vk, avm, tmf0);
  mpc_init3(ri, xc, xcb);

  mpf_set0(zero);
  mpf_set0(ixm);

  mpf_read_decimal(vm, PAnsiChar(ansistring(labelededit1.Text + #00)));
  mpf_read_decimal(xm, PAnsiChar(ansistring(labelededit2.Text + #00)));
  if NonComplex.Checked = false then
    mpf_read_decimal(ixm, PAnsiChar(ansistring(labelededit3.Text + #00)));

  mpf_read_decimal(zz, PAnsiChar(ansistring(zeromg + #00)));

  setlength(redata, N + 1);
  setlength(imdata, N + 1);

  mpc_set_mpf(xcb, xm, ixm);                   // xcb 計算用 xの複素数

  series1.Clear;
  series2.Clear;
  series3.Clear;
  series4.Clear;
  series5.Clear;
  series6.Clear;

  ChartBackclear;

  berxe := 0;
  beixe := 0;
  if NonComplex.Checked = false then begin
    if ix >= 0 then istr := ' +' + floatTostr(ix) + 'i'
               else istr := ' ' + floatTostr(ix) + 'i';
    memo1.Lines.Append('v = ' + floatTostr(vin) + '  x = ' + floatTostr(xin) + istr);
  end
  else
    memo1.Lines.Append('v = ' + floatTostr(vin) + '  x = ' + floatTostr(xin));

  application.ProcessMessages;

  rk := 0;
  ik := 0;

  // ガンマ値配列計算
  if mpf_is_ge(vm, zero) then begin
    for k := 0 to KMmax do begin
      mpf_set_int(tmf, k);                    // k
      mpf_add(tmf, vm, tmf0);                 // v + k
      mpf_add(tmf0, oneahf, vk);              // vk= v + k + 3/2
      gammaMul(vk, PVG[k]);                   // Γ(n+k+3/2)
    end;
  end
  else begin
    for k := 0 to KMmax do begin
      mpf_set1(nd);
      mpf_set_int(tmf, k);                    // k
      mpf_add(tmf, vm, tmf0);                 // -v + k
      mpf_add(tmf0, oneahf, vk);              // vk= -v + k + 3/2
      if mpf_is_lt(vk, zero) then begin       // vk < 0 時 nxが整数か確認
        mpf_int(vk, tmf);                     // int(vk);
        mpf_sub(vk, tmf, nd);                 // vk - int(vk) vkが負の整数だったら nd = 0
      end;
      if mpf_is_ne(vk, zero) and mpf_is_ne(nd, zero) then // vkが負の整数の時は計算しない
        gammaMul(vk, PVG[k])                  // Γ(n+k+3/2)
      else
        mpf_set0(PVG[k]);                     // Γ(n+k+3/2) = 0
    end;
  end;

  // ゼロ近傍の計算 無限大符号設定
  // x=0  v < -1   ±∞
  mpf_chs(one, tmf0);                         // -1
  if mpc_is0(xcb) and mpf_is_lt(vm, tmf0) and mpf_is_ne(nd, zero) then begin
    xvtoxc(x0m, ix, xc);                      // x double to mpc
    Hvx(vm, xc, ri);                          // ゼロ近辺のプラス側x計算
    if mpf_is_gt(ri.re, zero) then rk := infinity;      // x=0時の無限大±符号設定
    if mpf_is_lt(ri.re, zero) then rk := -infinity;
    if mpf_is_gt(ri.im, zero) then ik := infinity;
    if mpf_is_lt(ri.im, zero) then ik := -infinity;
  end;

  // 表示値の計算
  Hvx(vm, xcb, ri);                         // xcb 複素数
  mpf_div(one, two, tmf);                   // 0.5
  mpf_abs(vm, avm);                         // |v|
  mpf_frac(avm, avm);                       // avm - int(avm)
  // v = *.5 x=0 の時 Hvx Lvx = 0
  if mpc_is0(xcb) and mpf_is_eq(avm, tmf) then begin
    if checkbox3.Checked = false then
      memo1.Lines.Append( string('Hv(x) = ' + mpf_decimal(ri.re, 50)))
    else
      memo1.Lines.Append( string('Lv(x) = ' + mpf_decimal(ri.re, 50)));
    if NonComplex.Checked = false then  begin
      if mpf_is_ge(ri.im, zero) then
        memo1.Lines.Append( string('          +' + mpf_decimal(ri.im, 50) + 'i'))
      else
        memo1.Lines.Append( string('           ' + mpf_decimal(ri.im, 50) + 'i'));
    end;
    berxe := mpf_todouble(ri.re);           // グラフ用double値
    beixe := mpf_todouble(ri.im);
  end
  else
  // x=0  v < -1  vk = 整数  の場合  ±∞
    if mpc_is0(xcb) and mpf_is_lt(vm, tmf0) and mpf_is_ne(nd, zero) then begin
      berxe := rk;                          // ±無限大時
      beixe := ik;
      if checkbox3.Checked = false then
        memo1.Lines.Append('Hv(x) =  ' + floatTostr(berxe))       // 無限大表示
      else
        memo1.Lines.Append('Lv(x) =  ' + floatTostr(berxe));      // 無限大表示
      if NonComplex.Checked = false then
        if beixe >= 0 then
          memo1.Lines.Append('            +' + floatTostr(beixe) + 'i');
    end
    else begin
      if checkbox3.Checked = false then
        memo1.Lines.Append( string('Hv(x) = ' + mpf_decimal(ri.re, 50)))
      else
        memo1.Lines.Append( string('Lv(x) = ' + mpf_decimal(ri.re, 50)));
      if NonComplex.Checked = false then begin
        if mpf_is_ge(ri.im, zero) then
          memo1.Lines.Append( string('          +' + mpf_decimal(ri.im, 50) + 'i'))
        else
          memo1.Lines.Append( string('           ' + mpf_decimal(ri.im, 50) + 'i'));
      end;
      berxe := mpf_todouble(ri.re);         // グラフ用double値
      beixe := mpf_todouble(ri.im);
    end;

  if checkbox2.Checked = true then begin          // グラフ無だったら終了
    Chart1.Canvas.Font.Style := [fsBold];
    Chart1.Canvas.Font.size := 8;
    Chart1.Canvas.TextOut(170, 115,'グラフ無し');
    goto EXT;
  end;

  // グラフ表示
  ymin := 0;
  ymax := 0;
  if checkbox3.Checked = false then begin
    xmin := round(xin) - 10;
    xmax :=  xmin + 20;
    dx := 20 / N;
  end
  else begin
    xmin := round(xin) - 5;
    xmax :=  xmin + 10;
    dx := 10 / N;
  end;

  for i := 0 to N do begin
    x := dx * i + xmin;
    xvtoXc(x, ix, xc);
    if not mpc_is0(xc) then begin
      Hvx(vm, xc, ri);                              // xcb 複素数
      yre := mpf_todouble(ri.re);
      yim := mpf_todouble(ri.im);
      if ymin > yre then ymin := yre;
      if ymax < yre then ymax := yre;
      if ymin > yim then ymin := yim;
      if ymax < yim then ymax := yim;
      redata[i] := yre;
      imdata[i] := yim;

      if x < 0 then begin
        series1.AddXY(x, yre);
        if NonComplex.Checked = false then
          series2.AddXY(x, yim);
      end
      else begin
        series5.AddXY(x, yre);
        if NonComplex.Checked = false then
          series6.AddXY(x, yim);
      end;
    end;
  end;

  for i := 0 to N do begin
    x := dx * i + xmin;
    if x < 0 then begin
      series1.AddXY(x, redata[i]);
      if NonComplex.Checked = false then
        series2.AddXY(x, imdata[i]);
    end;
    if x = 0 then begin
      xv := abs(berxe);
      if xv = infinity then begin
        if berxe < 0 then series5.AddXY(0, ymin - 3)
                     else series5.AddXY(0, ymax + 3);
      end;
    end;
    if x > 0 then begin
      series5.AddXY(x, redata[i]);
      if NonComplex.Checked = false then
        series6.AddXY(x, imdata[i]);
    end;
  end;

  xv := abs(berxe);
  if xv <> infinity then begin
    if NonComplex.Checked = false then
      series4.AddXY(xin, beixe);
    series3.AddXY(xin, berxe);
  end
  else begin
    if berxe < 0 then series3.AddXY(0, ymin - 3)
                 else series3.AddXY(0, ymax + 3);
  end;

  application.ProcessMessages;
  chart1coment;
EXT:
  mpf_clear4(xm, vm, xvm, ixm);
  mpf_clear3(nd, tmf, zz);
  mpc_clear3(ri, xc, xcb);
  mpf_clear3(vk, avm, tmf0);
end;


procedure TForm1.FormCreate(Sender: TObject);
var
  i : integer;
  N, D, tmp : mp_float;
  bt : Tbitmap;
begin
  mpf_set_default_decprec(180);   // 有効桁数180桁 50桁の精度に必要です。
  setlength(BM, NB + 1);          // ベルヌーイ数配列
  setlength(FA, KMmax + 1);       // Γ(i+3/2)配列
  setlength(PVG, KMmax + 1);      // +vΓ


  for i := 0 to NB do mpf_init(BM[i]);
  for i := 0 to KMmax do mpf_init(FA[i]);
  for i := 0 to KMmax do mpf_init(PVG[i]);


  mpf_init3(N, D, tmp);
  mpf_init4(zero, one, two, four);
  mpf_init3(three, pai, log_2pis2);
  mpf_init(oneahf);

  mpf_set0(zero);
  mpf_set1(one);
  mpf_set_int(two, 2);
  mpf_set_int(three, 3);
  mpf_set_int(four, 4);
  mpf_set_pi(pai);
  mpf_div(three, two, oneahf);    // 3/2

  mpf_mul(pai, two, tmp);         // 2π
  mpf_ln(tmp, tmp);               // ln(2π)
  mpf_div(tmp, two, log_2pis2);   // ln(2π)/2


  Bernoulli_number_BigInteger;    //  ベルヌーイ数作成

  for i := 0 to NB do begin
    mpf_read_decimal(N, PAnsiChar(ansistring(NumeratorString[i] + #00)));
    mpf_read_decimal(D, PAnsiChar(ansistring(DenominatorString[i] + #00)));
    mpf_div(N, D, BM[i]);
  end;

  for i := 0 to KMmax do begin    // Γ(i+3/2)配列作成
    mpf_add_int(oneahf, i, tmp);
    gammaMul(tmp, FA[i])
  end;

  memo1.Clear;

  bt := Tbitmap.Create;
  bt.Width := Chart1.Width;
  bt.Height := Chart1.Height;
  Chart1.BackImage.Bitmap := bt;
  bt.Free;

  ChartBackclear;
  if NonComplex.Checked = true then begin
    Labelededit3.Visible := False;
    Labelededit2.Left := 460;
  end;

  mpf_clear3(N, D, tmp);
end;

procedure TForm1.NonComplexClick(Sender: TObject);
begin
  if NonComplex.Checked = true then begin
    Labelededit3.Visible := false;
    Labelededit2.Left := 460;
  end
  else begin
    Labelededit3.Visible := True;
    Labelededit2.Left := 432;
  end;
end;

procedure TForm1.FormClose(Sender: TObject; var Action: TCloseAction);
var
  i : integer;
begin
  for i := 0 to NB do mpf_clear(BM[i]);
  for i := 0 to KMmax do mpf_clear(FA[i]);
  for i := 0 to KMmax do mpf_clear(PVG[i]);
  mpf_clear4(zero, one, two, four);
  mpf_clear3(three, pai, log_2pis2);
  mpf_clear(oneahf);
end;

end.


download Struvesche_funktion.zip

  三角関数、逆三角関数、その他関数 に戻る