三分探索の実装をライブラリ化した

二分探索と比べて、区間の内分点の取り方やら、下に凸か上に凸かやら、地味に考えることが多いので汎用性を目指してライブラリ化してみた。

前提

• 三分探索とは、ある区間[l, r]において、極値がただ一つだけ存在するとき、その極値の近似値を区間の長さをnとして、およそ O(log n) で求める手法
• 二分探索が単調グラフにおける最大値、最小値を求めるものに対して、三分探索は凸関数における最大値、最小値を求めるためのもの
• 三分探索の実装では、以下のことを行っている
  • 区間を等間隔に3つに内分する2点を取り、f(x) が極値とより離れている方を狭める、ということを繰り返す
  • なぜ上記によって区間内に極値が収まるのかは、以下の記事に図解がある

三分探索のわかりやすい解説

juppy.hatenablog.com

実装

• l, rに、対象の区間を指定する
• 第3引数fに、対象の関数を与える
• 上に凸な関数なら、第4引数is_downward_convexに false を指定する

/// 三分探索
double ternary_search(double l, double r, const function<double(double)> f, const bool is_downward_convex=true){
  int loop=500;
  if(is_downward_convex){ // 下に凸
    while(loop--){
      auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3;
      if(f(mid_l)>=f(mid_r)) l=mid_l;
      else r=mid_r;
    }
  }else{                  // 上に凸
    while(loop--){
      auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3;
      if(f(mid_l)<=f(mid_r)) l=mid_l;
      else r=mid_r;
    }
  }
  return l;
}

検証用問題

atcoder.jp

呼び出し元

上記問題を例に、以下に実際の使い方を載せる。

double p;

/* ternary_search */

signed main()
{
  cin>>p;
  auto f=[&](double x){
    return x+p/pow(2,x/1.5);
  };
  auto ans_x=ternary_search(0,1e18,f);
  printf("%.15lf\n", f(ans_x));
}