問題はこちら
No.456 Millions of Submits! - yukicoder
いろいろな知見が得られた
見出し
・愚直二分探索(TLE)
・W関数を使う二分探索
・ニュートン法(賢い)
・W関数をなんか凄い方法で求める
a,bの一方が0のときは自明。そうでないとする
とおく
b≠0の仮定からf(1)=0であり、
となることと、x>1でf(x)>0であることから、x>1でf(x)は狭義単調増加。
よって、t>0に対しはx>1にただ1つ解をもつので二分探索するだけ……?
・二分探索
区間の左端は1、右端は、x>1ではa,b≧1よりなのでサンプルケースから5.8などを取れば良いことがわかる
int main(){ int a,b,m; double t; scanf("%d",&m); while(m--){ scanf("%d%d%lf",&a,&b,&t); if(b==0)printf("%.10f\n",pow(t,1./a)); else if(a==0)printf("%.10f\n",exp(pow(t,1./b))); else{ int i; double nl=1,nr=5.8,nm,fm; for(i=0;i<33;i++){//5.8*2^-33≒6.7*10^-10 nm=(nl+nr)/2; fm=pow(nm,a)*pow(log(nm),b); fm<t?nl=nm:(nr=nm); } printf("%.10f\n",nm); } } return 0; }
これはTLE。
m=3*10^5で1600ms程度、m=10^6の最大ケースではTLE KO(制限時間が4.5秒なので5.5秒以上かかったということ)
は? 10^6*33だから間に合うんじゃないの?
……logとexpの計算が重いせいだが、この時点では思い至らず
しかたがないので、数学的に解が求まったりするのかなと思い式変形をしていると「これってランベルトのW関数っぽい」と気づく
ランベルトのW関数とは、の逆関数のこと。詳細はwikipediaなどを参照のこと
ランベルトのW関数 - Wikipedia
ということでゴリゴリ数学やるよ
====
問題:の(x>1における唯一の)解を求めよ
と置換した
両辺1/b乗した
と置換した
両辺にを乗じた
W関数の定義を適用した
置換を戻した
置換を戻した
====
ということで、W関数の値を素早く求める問題に帰着された
・W関数の二分探索
……帰着されたけどさあ、W関数の値も結局は二分探索するしかなくない?
はx>0で正の値を取る狭義単調増加関数だから、その逆関数Wもx>0で正の値をとる狭義単調増加関数
であり、
W(100)≒3.3856なので(参考)、区間は[0,3.4]などでよい
double W(double x){ double nl=0,nr=3.4,nm,fm; int i; for(i=0;i<50;i++){ nm=(nl+nr)/2; fm=nm*exp(nm); fm<x?nl=nm:(nr=nm); } return nm; } int main(){ double t; int a,b,m; scanf("%d",&m); while(m--){ scanf("%d%d%lf",&a,&b,&t); if(b==0)printf("%.10f\n",pow(t,1./a)); else if(a==0)printf("%.10f\n",exp(pow(t,1./b))); else printf("%.10f\n",exp(W(pow(t,1./b)*a/b)*b/a)); } return 0; }
二分探索の回数は適当に50回としたが通った。
なんで……さっきと同じく二分探索してるだけじゃん……と思ったのだが、powとlogの回数が違うせいだということを解説を読んでようやく理解した。
適当に設定したW関数の精度だが、このあとexpに食わせたりしているので、どのくらいの精度が必要かはよく考えないといけない
次のようにして評価できる
元の問題に対する答えは必ず1以上になるので、誤差は相対誤差のみを考えれば良い
任意のxとε>0に対してが成立するので
(∵平均値の定理とexpの単調増加性)
のとき
よって相対誤差は未満
さらには十分0の近いので、は1.1で真におさえられるとしてよく(実際ε<9*10^-3なら成立)
ということでε<1/(1.1*10^10)くらいならよいことになる
二分探索の回数を決め打ちでやるなら3.4/2^k<1/(1.1*10^10)を解いて36回以上やればよい。
・ニュートン法
tailsさんの提出をみて気づいた。
の対数をとり、と置換して
この方程式の解をニュートン法で求める
とおけばなので、次のように実装できる
double lt=log(t);//定数なので先に計算しておく double x1=0.1;//ニュートン法の初期値 for(i=0;i<9;++i)x1-=(a*x1+b*log(x1)-lt)/(a+b/x1); x=exp(x1);
(12/09 18:00初期値を修正)
2次収束なのでものの数回で十分な精度が得られる。
短い+速いの大正義である。
(対数を取らずに元のfに対してニュートン法をしてもほぼ同じものが得られる)
そもそも元の問題は「方程式の解を求めること」だったのに、なぜニュートン法にすぐ思い至らなかったのか。
研鑽が足りない。
この解法より短く出来る気がしないのでゴルフはしてない
・W関数のすごいやつ
解説で言及されてるHalley iterationというやつをググってみた
Halley's method - Wikipedia
次の漸化式によって方程式f(x)=0を高速に解くアルゴリズムであり3次収束するらしい。強い。
どこからこの漸化式が生えたのかと言えば、実はニュートン法の一般化らしい
ググって3番目に出てきたpdfがぱっと見で分かりやすかった
さて、これを使ってW関数を高速に求める。W(z)の値はの解だから、とおけばよい
このときとなるので、漸化式に投げ込むと
ここまで計算した上で、想定解のコードを眺めてみると、この漸化式を
と変形して使っていることが分かる。
初期値には、z<1では、冪級数展開を適当に打ち切ったもの(?)、z>1ではlog(z)を使っているようだ。
つまりライブラリはるだけゲーだった……?(過言)