読者です 読者をやめる 読者になる 読者になる

メモ

yukicoderで遊んでいる競プロゆるふわ勢

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)を使っているようだ。

つまりライブラリはるだけゲーだった……?(過言)