メモ

yukicoderでゆるふわgolf

yukicoder No.477 MVP

問題はこちら
No.477 MVP - yukicoder

(「N人が選挙権を持っている時、上位K位以内に入るには何票獲ればよいか」という問題と同値)
「上位K位以内に入る」⇔「自分より上にK人以上いる、ということはない」なので、N/(K+1)より真に大きなダメージを与えればよい
N/(K+1)より真に大きな整数値はfloor(N/(K+1))+1で与えられるので

int main(){
	long N,K;
	scanf("%ld%ld",&N,&K);
	printf("%ld",N/(K+1)+1);
	return 0;
}


K+1は-~KなのでN/(K+1)は-N/~Kとなる

long N,K;
main(){
	scanf("%ld%ld",&N,&K);
	N=!printf("%ld",1-N/~K);
}

63B

yukicoder No.476 正しくない平均

問題はこちら
No.476 正しくない平均 - yukicoder

やるだけ
入力が全て整数なので、平均×個数=合計 となっているかどうかで調べると、小数を使わずに済む
10^9*10=10^10>2^32なので32bitでは足りない事に注意(撃墜ケースも用意されている)

int main(){
	int a,n,s,i,x;
	scanf("%d%d",&n,&a);
	s=0;
	for(i=0;i<n;i++){
		scanf("%d",&x);
		s+=x;
	}
	puts((long)n*a==s?"YES":"NO");
	return 0;
}

さて、できるだけ変数を減らすことを考えると、n,aを特別扱いせずにまとめてしまいたい。
やるべきことは n*a-(残りの和) が0かどうかを判断することなので
1回目   s=x
2回目   s*=x
3回目以降 s-=x
とすればよさそう

long s,x;
main(i){
	for(;~scanf("%d",&x);--i?~i?s-=x:(s*=x):(s=x));
	x=!puts(s?"NO":"YES");
}

最初の代入をちょっと工夫するとすこし縮む

long s,x;
main(i){
	for(;~scanf("%d",--i?&x:&s);~i?s-=x:(s*=x));
	x=!puts(s?"NO":"YES");
}

やはり3分岐が長く見える
そこでこれをどうにかまとめる方法がないか考えると
最初    s=1
1,2回目  s*=x
3回目以降 s-=x
とすれば良いことに気付く

long s=1,i;
main(x){
	for(;~scanf("%d",&x);++i<3?s*=x:(s-=x));
	i=!puts(s?"NO":"YES");
}

82B
ちなみにカッコを外してみても同じ長さになる

++i<3?s*=x:(s-=x)
s+=++i<3?~-x*s:-x

sがintでよければ、mainの第1引数にすることで初期化が省略できるので更に縮むのだけど。

ゴルフに際しては、過去にどこかでにたような問題を見たことがあったので、わりと簡単に縮めることができた

yukicoder No.478 一般門松列列

問題はこちら
No.478 一般門松列列 - yukicoder

要するに「長さnの数列であって、その長さ3の連続部分列のうち、門松列でないものがちょうどk個であるようなもの」を答えれば良い
サンプル3を眺めると、どうやら「1,3,2,4」を繰り返したものは必ず門松列列になるようだ。これでkが0の場合は解決できた。
これの先頭に1をつけると、最初の「1,1,3」は門松列で無くなるが、それ以降は影響なく門松列なので(先頭に1つ増やした分、末尾を1つ削れば)k=1の場合も解決できる。
同様に先頭に1をk個つけることで一般の場合も解決できる

int main(){
	int n,k,i;
	int a[]={1,3,2,4};
	scanf("%d%d",&n,&k);
	for(i=0;i<k;i++)printf("1 ");
	n-=k;
	for(i=0;i<n;i++)printf("%d ",a[i%4]);
	return 0;
}

「1,3,2,4」を出力する部分だが、これは大小関係さえこうなればよいので、例えばi%4*3%5により達成することが出来る

i%4 i%4*3%5
0 0
1 3
2 1
3 4

余計な変数を使わないよう、nをデクリメントしながらループするという形にする
そうすると「先頭k個を~」という処理はめんどくさいので、「末尾k個を~」と書き換え、以下のようなものができる

n;main(k){for(scanf("%d%d",&n,&k);n--;)printf("%d ",(n>k?n:k)%4*3%5);}

70B

yukicoder No.457 (^^*)

問題はこちら
No.457 (^^*) - yukicoder


・O(|S|^3)
各カッコ対に対して、その間に'^^*'があるかどうか判定する
当然TLE

・O(|S|^2)
各'('について、そこから右に見て'^^*'があった以降にある')'の数を数える
(*^^)も同様
17/01/31修正

char s[10010];i,j,L,R;
int main(){
	gets(s);
	for(i=0;s[i];i++)if(s[i]=='('){
		
		//(^^*)
		j=i;
		while(s[j]!='\0'&&s[++j]!='^');
		while(s[j]!='\0'&&s[++j]!='^');
		while(s[j]!='\0'&&s[++j]!='*');
		while(s[j]!='\0'){
			if(s[j]==')')L++;
			j++;
		}
		
		//(*^^)
		j=i;
		while(s[j]!='\0'&&s[++j]!='*');
		while(s[j]!='\0'&&s[++j]!='^');
		while(s[j]!='\0'&&s[++j]!='^');
		while(s[j]!='\0'){
			if(s[j]==')')R++;
			j++;
		}
		
	}
	printf("%d %d",L,R);
	return 0;
}


・O(|S|)
よく考えれば上でやったことをDPにすればO(|S|)だった。
左から順に見ていって ( , (^ , (^^ , (^^* ,(^^*) が各々いくつあるかを覚える
右も同様

int main(){
	int i,dpL[5]={},dpR[5]={};
	char s[10010];
	gets(s);
	for(i=0;s[i];i++){
		switch(s[i]){
			case '(':
				dpL[0]++;
				dpR[0]++;
				break;
			case '^':
				dpL[2]+=dpL[1];dpL[1]=dpL[0];dpL[0]=0;
				dpR[3]+=dpR[2];dpR[2]=dpR[1];dpR[1]=0;
				break;
			case '*':
				dpL[3]+=dpL[2];dpL[2]=0;
				dpR[1]+=dpR[0];dpR[0]=0;
				break;
			case ')':
				dpL[4]+=dpL[3];
				dpR[4]+=dpR[3];
		}
	}
	printf("%d %d",dpL[4],dpR[4]);
	return 0;
}


とりあえずゴルフした限りではO(|S|^2)の方が短くなった

・O(|S|)
超素直に

i,d[10];
main(){
	for(;read(0,&i,1);)
	//')'のとき
	i&1?d[4]+=d[3],d[9]+=d[8]:
		//'^'のとき
		i&4?d[2]+=d[1],d[1]=*d,*d=0,d[8]+=d[7],d[7]=d[6],d[6]=0:
			//'('のとき
			~i&2?++*d,d[5]++:
				//'*'のとき
				i>11?d[3]+=d[2],d[2]=0,d[6]+=d[5],d[5]=0:
					0;
	i=!printf("%d %d",d[4],d[9]);
}

206B


・O(|S|^2)
まずざっくりと。
tL,tRで状態を保存して遷移をチェックする

char s['~~'];i,j,L,R,tL,tR;
main(){
	for(gets(s);s[i];tL=tR=0,i++)for(j=i;s[++j];)if(s[i]=='('){
		tL<2&s[j]=='^'|tL==2&s[j]=='*'?++tL:s[j]==')'?tL+=tL>2,L+=tL>3:0;
		tR<3&tR>0&s[j]=='^'|!tR&s[j]=='*'?++tR:s[j]==')'?tR+=tR>2,R+=tR>3:0;
	}
	i=!printf("%d %d",L,R);
}

これをぎゅっとする
'^'が94、'*'が42、'('が40、')'が41

tL<2&s[j]=='^'|tL==2&s[j]=='*'?++tL:s[j]==')'?tL+=tL>2,L+=tL>3:0;
tL>1|s[j]-94&&tL-2|s[j]%7:s[j]==')'&&tL>2?L+=++tL>3:0:++tL;
tL>1|s[j]-94&&tL-2|s[j]%7:s[j]&tL>2?L+=++tL>3:0:++tL;

if(s[i]=='(')tL>1|s[j]-94&&tL-2|s[j]%7:s[j]&tL>2?L+=++tL>3:0:++tL;
tL>1|s[j]-94&&tL-2|s[j]%7:s[j]&tL>2?L+=++tL>3:0:s[i]%8||++tL;

ということでループ圧縮も一緒にして

#define f(X,Y,Z,W)X|s[j]-94&&Y|s[j]%7?s[j]&Z>2?W+=++Z>3:0:s[i]%8||++Z,
char s['~~'];
i,j,p,q,L,R;
main(){
	for(gets(s);s[++j]?f(p>1,p-2,p,L)f(q-1&q-2,q,q,R)1:(p=q=0,s[j=++i]););
	i=!printf("%d %d",L,R);
}

197B

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

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