C++17でLegendre関数を使ってみる。

この記事の内容

C++17では特殊関数が標準の仲間入りした。自分の記事の中で少し触れたこともある散々いらない子扱いした。特殊関数の中でも比較的簡単でかつ有用性の高いLegendre関数(多項式)を使ってC++17の特殊関数の使い方等を解説する。

ガウスルジャンドル求積

特殊関数はある範囲での直交性を持っている。Legendre関数の場合$-1\sim1$で

\int_{-1}^{1}P_n(x)P_m(x)dx=\delta_{nm}

と表されます。この直交性を持つために補完関数として優秀なので数値積分する際の補完関数としてよく用いられる。特にLegendre関数を用いた補完は有名で数値積分を数ステップの計算で$\delta \sim10^{-8}$の誤差を達成できます。
台形則を用いるとその誤差は$\sim h^2$、シンプソン則を使うと$\sim h^4$なのでシンプソン則で$\delta \sim 10^{-8}$を得ようとするとだいたい100ぐらいの分割が必要になります。これを見るといかにガウスルジャンドル法が優秀か解ると思います。しかし一方で補間法なので高次項が積分範囲の端で影響を受けやすく、端の方に振動のような形を作り出すルンゲ現象というものもあるので過信はできません。
ここでは積分法の優劣は問題にせずにいかに実装していくかを見ていきます。
std::legendreはconstexprではないのでコンパイル時定数化はできません。しかしガウスルジャンドル法はlegendre関数の0点とそれに対応した重みがあれば計算できるのでlegendre関数自体がconstexprでなくてもconstexpre化できます。
倍精度doubleで計算するときに倍精度以下の計算誤差に抑えたいので計算自体には4倍精度のlong doubleを使います。
以下が計算部分です。

main.cc

#include <bits/stdc++.h>
#include "solver.hh"

int main(){
  const size_t N=100;
  size_t n=1;
  std::vector<std::vector<long double> > vec_zero_points(1), vec_weights(1);
  for( ; n<100; n++ ){
    std::vector<long double> rough;
    const long double delta=1.0/(100*n);

    for( long double x=-1; x+delta<1.0; x+=delta ){
      if( std::legendrel(n, x)*std::legendrel(n, x+delta)<0 ) rough.push_back(x+0.5*delta);
    }
    if( rough.size()!=n ) return 0;

    auto func=[n](long double x)->long double { return std::legendrel(n, x); };
    std::vector<long double> zero_points;
    for( size_t i=0; i<rough.size(); i++ ){
      long double min=rough[i]-delta, max=rough[i]+delta;
      zero_points.push_back(math::solver::bynaryl(func, min, max, 1.0e-19));
    }
    std::vector<long double> weights;
    for( size_t i=0; i<zero_points.size(); i++ ){
      weights.push_back(2.0*(1.0-zero_points[i]*zero_points[i])/(n*n*std::legendrel(n-1, zero_points[i])*std::legendrel(n-1, zero_points[i])));
    }
    vec_zero_points.push_back(zero_points);
    vec_weights.push_back(weights);
  }

最初、ラフに0点を探させその後に2分法を使って高精度な0点を計算させています。std::legendre関数は$-1\sim1$の定義らしく少しでも$-1\sim1$から漏れるとstd::domain_errorを飛ばすので範囲外を踏まないようにラフサーチのfor文をfor( long double x=-1; x+delta<1.0; x+=delta ){としています。本当に余計なお世話です。
収束の速さだけで言えばニュートン(ラフソン)法のほうが優秀ですが微分が入るため誤差を小さくしたい場合は2分法のほうが優秀になります（もちろん数式微分などで微分自体の精度を上げれば同程度の精度は出るはずだが数式微分の実装の難しさを考えれば実行速度が遅くても2分法で書くほうが早い)。
2分法の部分です。C++14で代数方程式を解くで解説したものをlong double化しました。テンプレートにして引数で!とも考えましたが今回は標準ライブラリっぽく末尾にlをつけることでlong doubleであることを示しました。4倍精度でだいたい30桁ぐらいの精度はあるはずですが20桁ぐらいから精度が上がらなかったので19桁で打ち切ってます、それでも倍精度なら十分なので満足しています。

solver.hh

namespace math{
  namespace solver{
    long double bynaryl(const std::function<long double(long double)> &f, long double min, long double max, long double error=10e-10){
      if( min>max ) std::swap(min, max);
      if( f(min)*f(max)>0 ){
        std::cout<<"!!! bynarySearch  bad initial parameter ["<<min<<":"<<max<<" !!!"<<std::endl;
        return 0.f/0.f;
      }
      long double half=0.5;
      long double mid=half*(min+max);
      long double val_min=f(min);
      long double val_mid=f(mid);
      long double val_max=f(max);

      while( max-min>error ){
        if( val_mid*val_max>0 ) max=mid, val_max=val_mid;
        else min=mid, val_min=val_mid;
        mid=half*(max+min);
        val_mid=f(mid);                                                                                     
      }
      return mid;
    }
  }
}

#endif

ガウスルジャンドル法の重みですが$w_i=\frac{2}{((1-x_i^2)P'_n(x_i))^2}$で微分を含むので漸化式

(1-x^2)P_n'(x)=nxP_n(x)-nP_{n-1}(x)

を使い微分を消去する。$x_i$はルジャンドル関数の0点なので右辺第一項は消えて

w_i=\frac{2(1-x_i^2)}{(nP_{n-1}(x_i))^2}

隣ります。式で書くと簡単だがコードにするとめっちゃ長くて見づらい...orz、コードが長いのは悪だ！それで何回ポカをやったことか！！！
あとは計算した値をテキストダンプさせます。

main.cc

  std::ofstream ofs("gauss_legendre_param.hh");
  ofs<<"#ifndef GAUSS_LEGENDRE_PARAM_HH"<<std::endl;
  ofs<<"#define GAUSS_LEGENDRE_PARAM_HH"<<std::endl;

  ofs<<"namespace math{ namespace integral{ namespace gauss_legendre{"<<std::endl;
  ofs<<"constexpr double zero_points["<<n<<"]["<<n-1<<"] = {"<<std::endl;
  for( size_t i=0; i<n; i++ ){
    ofs<<"{ ";
    for( size_t j=0; j<n-1; j++ ){
      long double val=0;
      if( j<vec_zero_points[i].size() ) val=vec_zero_points[i].at(j);
      if( fabs(val)<1.0e-16 ) val=0;

      if( j==n-2 ) ofs<<std::setprecision(16)<<val;
      else ofs<<std::setprecision(16)<<val<<", ";
    }
    if( i==n-1 ) ofs<<" }"<<std::endl;
    else ofs<<" },"<<std::endl;
  }
  ofs<<"};"<<std::endl;

  ofs<<"constexpr double weights["<<n<<"]["<<n-1<<"] = {"<<std::endl;
  for( size_t i=0; i<n; i++ ){
    ofs<<"{ ";
    for( size_t j=0; j<n-1; j++ ){
      long double val=0;
      if( j<vec_weights[i].size() ) val=vec_weights[i].at(j);
      if( fabs(val)<1.0e-16 ) val=0;

      if( j==n-2 ) ofs<<std::setprecision(16)<<val;
      else ofs<<std::setprecision(16)<<val<<", ";
    }
    if( i==n-1 ) ofs<<" }"<<std::endl;
    else ofs<<" },"<<std::endl;
  }
  ofs<<"};"<<std::endl;
  ofs<<"}}}"<<std::endl;
  ofs<<"#endif"<<std::endl;

  return 0;
}

C++ヘッダーファイル形式で名前空間math::integral::gauss_legendreでくくっています。
使う側はこのダンプしたファイルをincludeして計算させます。2次元配列なので上半分はダミーデータが入っています。一応、$n=100$まで出していますがまず使うことはないので10程度で十分です(それで収束しない時は多分いくら精度を上げても無駄)。

使ってみた感想

今回は実際使うためのパラメーターの計算ということで半分書捨であまり綺麗には書いていませんがチャチャッと特殊関数を計算できるのはありがたいです。$n=100$まで計算しましたが規格的には$n=128$までが保証の範囲内でそれ以上は実装によるらしいです。例外を投げやすいのは気をつけるべき点です。
さて、$n=100$のルジャンドル関数ですが当然、100次の多項式です、99次はないが98、96はあるという複雑な形で$|x|=1$の周りだと変化が大きくなる結構やっかいな関数です。以前の記事で高次項のゼロ点が計算できなかったのですがその原因はどうやら$x=1$周辺でどこかの項が桁落ちしてるっぽいです。こういった関数は超高次多項式と言われ桁落ち問題が起こりやすいそうです。
ただが多項式、されど多項式というわけでそうったものに特化したアルゴリズムもいくつかありました。
今回はGCCを使いましたがC++の実装がなかったので初等関数同様にアセンブラで書いて高速化してるんでしょう。$-1\sim1$の範囲とはいえ-20桁の精度を安定して出せるのは結構大変なはずです。特に宗教上問題がなければ使えるんじゃないかと思いました。
ガウスホニャララの積分法はガウスラゲール、ガウス-チェビシェフ、ガウスエルミート等ありますがこの辺は全部実装されているので同様に使えるはずです。