yukicoder No.456 Millions of Submits!

問題はこちら
No.456 Millions of Submits! - yukicoder

いろいろな知見が得られた

見出し
・愚直二分探索（TLE）
・W関数を使う二分探索
・ニュートン法（賢い）
・W関数をなんか凄い方法で求める

a,bの一方が0のときは自明。そうでないとする
$f(x)=x^a(\log x)^b$ とおく
b≠0の仮定からf(1)=0であり、
$f'(x)=(\frac{a}{x}+\frac{b}{x\log x})f(x)$ となることと、x>1でf(x)>0であることから、x>1でf(x)は狭義単調増加。
よって、t>0に対し $f(x)=t$ はx>1にただ1つ解をもつので二分探索するだけ……？

・二分探索
区間の左端は1、右端は、x>1ではa,b≧1より $f(x)\geqq x\log x$ なのでサンプルケースから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関数とは、 $xe^x$ の逆関数のこと。詳細はwikipediaなどを参照のこと
ランベルトのW関数 - Wikipedia
ということでゴリゴリ数学やるよ
＝＝＝＝
問題： $x^a(\log x)^b=t$ の（x>1における唯一の）解を求めよ
$x^a(\log x)^b=t$
$e^{ax_1}{x_1}^b=t$ 　　 $x=e^{x_1}$ と置換した
$e^{\frac{a}{b}x_1}x_1=t^{\frac{1}{b}}$ 　　両辺1/b乗した
$e^{x_2}\frac{b}{a}x_2=t^{\frac{1}{b}}$ 　　 $x_2=\frac{a}{b}x_1$ と置換した
$e^{x_2}x_2=\frac{a}{b}t^{\frac{1}{b}}$ 　　両辺に $\frac{a}{b}$ を乗じた
$x_2=W(\frac{a}{b}t^{\frac{1}{b}})$ 　　W関数の定義を適用した
$x_1=\frac{b}{a}W(\frac{a}{b}t^{\frac{1}{b}})$ 　　置換を戻した
$x=\exp(\frac{b}{a}W(\frac{a}{b}t^{\frac{1}{b}}))$ 　　置換を戻した
＝＝＝＝
ということで、W関数の値を素早く求める問題に帰着された

・W関数の二分探索
……帰着されたけどさあ、W関数の値も結局は二分探索するしかなくない？
$xe^x$ はx>0で正の値を取る狭義単調増加関数だから、その逆関数Wもx>0で正の値をとる狭義単調増加関数
$\frac{a}{b}t^{\frac{1}{b}}\leqq\frac{10}{1}10^{\frac{1}{1}}=100$ であり、
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に対して $|x-y|<\epsilon \Rightarrow |\exp(x)-\exp(y)|< \epsilon\exp(x+\epsilon)$ が成立するので
（∵平均値の定理とexpの単調増加性）
$|W(\cdot)-z|<\epsilon$ のとき
$|\exp(\frac{b}{a}W(\cdot))-\exp(\frac{b}{a}z)|<\frac{b}{a}\epsilon\exp(\frac{b}{a}(W(\cdot)+\epsilon))\\ =\frac{b}{a}\epsilon\exp(\frac{b}{a}\epsilon) \exp(\frac{b}{a}W(\cdot))$
よって相対誤差は $\frac{b}{a}\epsilon\exp(\frac{b}{a}\epsilon)$ 未満
さらに $\frac{b}{a}\epsilon$ は十分0の近いので、 $\exp(\frac{b}{a}\epsilon)$ は1.1で真におさえられるとしてよく（実際ε<9*10^-3なら成立）
$\frac{b}{a}\epsilon\exp(\frac{b}{a}\epsilon)<11\epsilon$
ということでε<1/(1.1*10^10)くらいならよいことになる
二分探索の回数を決め打ちでやるなら3.4/2^k<1/(1.1*10^10)を解いて36回以上やればよい。

・ニュートン法
tailsさんの提出をみて気づいた。
$x^a(\log x)^b=t$ の対数をとり、 $x_1=\log x$ と置換して $ax_1+b\log x_1=\log t$
この方程式の解をニュートン法で求める
$g(x)=ax_1+b\log x_1-\log t$ とおけば $g'(x)=a+\frac{b}{x_1}$ なので、次のように実装できる

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次収束するらしい。強い。
$x_{n+1}=x_n-\frac{2f(x_n)f'(x_n)}{2(f'(x_n))^2-f(x_n)f''(x_n)}$
どこからこの漸化式が生えたのかと言えば、実はニュートン法の一般化らしい
ググって3番目に出てきたpdfがぱっと見で分かりやすかった

さて、これを使ってW関数を高速に求める。W(z)の値は $xe^x=z$ の解だから、 $f(x)=xe^x-z$ とおけばよい
このとき $f'(x)=(x+1)e^x\ \ f''(x)=(x+2)e^x$ となるので、漸化式に投げ込むと
$x_{n+1}=x_n-\frac{2(xe^x-z)(x+1)}{2e^x(x+1)^2-(xe^x-z)(x+2)}$
ここまで計算した上で、想定解のコードを眺めてみると、この漸化式を
$x_{n+1}=x_n-\frac{xe^x-z}{e^x(x+1)-\frac{1}{2}(xe^x-z)\frac{x+2}{x+1}}$ と変形して使っていることが分かる。
初期値には、z<1では、冪級数展開を適当に打ち切ったもの(?)、z>1ではlog(z)を使っているようだ。

つまりライブラリはるだけゲーだった……？（過言）

メモ

yukicoderでゆるふわgolf

yukicoder No.456 Millions of Submits!