2020-08-05

FFT

何回見ても非再帰 FFTの話が覚えられないのでメモ

1のN乗根をgとする。
$f(x)=\sum_{i=0}^{N-1}a_i x^i$ を離散フーリエ変換したい。
その結果は
$\hat{f}(x)=\sum_{i=0}^{N-1}f(g^i) x^i$
になる。ということで、全てのiについての $f(g^i)$ が求められればよい。
$f(g^i)=\sum_{k=0}^{N-1}a_k g^{ik}$
これを頑張って計算する。
ところで、再帰 FFTの気持ちには2種類あることをご存じでしたか？ぼくは今知りました。

kの偶奇で分ける

$\begin{eqnarray*} f(g^i)&=&\sum_{k:even}a_k g^{ik}&+&\sum_{k:odd}a_k g^{ik}\\ &=& \sum_{k'=0}^{N/2-1} a_{2k'} g^{i2k'}&+&\sum_{k'=0}^{N/2-1} a_{2k'+1} g^{i(2k'+1)}\\ &=& \sum_{k'=0}^{N/2-1} a_{2k'} (g^{2i})^{k'}&+&g^i \sum_{k'=0}^{N/2-1} a_{2k'+1} (g^{2i})^{k'}\\ &=& f_e(g^{2i})&+&g^i f_o(g^{2i}) \end{eqnarray*}$
ということで、fの偶数次の項、奇数次の項のみからなるN/2次の多項式をFFTできればよいので再帰的にやればできる

iの偶奇で分ける

$\begin{eqnarray*} f(g^{2i})&=&\sum_{k=0}^{N/2-1}a_k g^{2ik}&+&\sum_{k=N/2}^{N-1} a_k g^{2ik}\\ &=&\sum_{k=0}^{N/2-1}a_k g^{2ik}&+&\sum_{k=0}^{N/2-1} a_{k+N/2} g^{2i(k+N/2)}\\ &=&\sum_{k=0}^{N/2-1}a_k g^{2ik}&+&\sum_{k=0}^{N/2-1} a_{k+N/2} g^{2ik}\\ &=&\sum_{k=0}^{N/2-1}(a_k+a_{k+N/2}) g^{2ik}& & \end{eqnarray*}$
$\begin{eqnarray*} f(g^{2i+1})&=&\sum_{k=0}^{N/2-1}a_k g^{(2i+1)k}&+&\sum_{k=N/2}^{N-1} a_k g^{(2i+1)k}\\ &=&\sum_{k=0}^{N/2-1}(a_k g^k)g^{2ik}&+&\sum_{k=N/2}^{N-1} (a_k g^k) g^{2ik}\\ &=&\sum_{k=0}^{N/2-1}(a_k g^k) g^{2ik}&+&\sum_{k=0}^{N/2-1} (a_{k+N/2} g^{k+N/2}) g^{2i(k+N/2)}\\ &=&\sum_{k=0}^{N/2-1}(a_k g^k) g^{2ik}&+&\sum_{k=0}^{N/2-1}(- a_{k+N/2} g^k) g^{2i(k+N/2)}\\ &=&\sum_{k=0}^{N/2-1}(a_k-a_{k+N/2}) g^k g^{2ik}& & \end{eqnarray*}$
ということで、iが偶数のときは $\sum_{k=0}^{N/2-1}(a_k+a_{k+N/2})x^k$ 、奇数のときは $\sum_{k=0}^{N/2-1}(a_k-a_{k+N/2})g^k x^k$ という多項式がそれぞれFFTできればよいので再帰的にやればできる

非再帰 FFT

iの偶奇で分けるFFTを非再帰にするのは簡単。 $a_k+a_{k+N/2}$ と $(a_k-a_{k+N/2}) g^k$ を作っていけばいいだけなので。
欲張りにinplaceでやりましょう。適当に、前半に $a_k+a_{k+N/2}$ を、後半に $(a_k-a_{k+N/2}) g^k$ をまとめていくことにする。
とりあえず1回

これで前半と後半が別の問題になったので、後半のことは忘れて前半だけに注目してもう1回

ところでこうやってやっていくと、前半には「下のbitが0のもの」が、後半には「下のbitが1のもの」が集まっていくことになる。
要するにbit reverseされているので、最後にその辺をうまく並び替えればOK

kの偶奇で分けるFFTも、実は頑張ると非再帰にできる。
解説は省略、原理はこう（再帰でやっていたものをボトムアップにやっていく気持ち）

参考にしました
https://tomorinao.blogspot.com/2018/10/various-fft.htmltomorinao.blogspot.com

2020-02-12

純粋培養競プロerが応用情報技術者試験に合格するまで

ポエム

Q.これはなに笑　
A.う笑　

・スペック
理系だが情報系ではない
競プロ以外のプログラミングをしたことはない

・2019年2月
基本情報技術者試験とかいうのを受けておくといいことがあるらしい
適当に参考書を1冊買う

・2019年3月
とりあえずノー勉で過去問を1回分解いてみると午前が5割弱、午後が1ミス(98点?)だった。
合格ラインは6割なのですぐ届く気がするが、4択問題なので、実は3割しか分かってなくても30%+70%×1/4=47.5%くらいの点が取れるんですよね～

・2019年4月
基本情報技術者試験ドットコムだかなんとかいうサイトの過去問をランダムに表示する機能を使ってとりあえず300問くらい解いてみる。
すると7割くらい解けるようになってきた気がするので過去問を2回分やってみるとどっちも7割前後とれている。
試験勉強終わり！　一度も開いてない参考書さん…………

・結果
午前が70ちょっと、午後満点。

・2019年8月
気付いたら申し込みの締め切りが翌日に控えていた。慌てて応用情報技術者試験に申し込みをする。

・2019年10月
気付いたら試験日が来週に迫っていた。
とりあえず午前の過去問を解いてみると7割くらい解ける。基本情報と難易度変わらね～　勝ったなガハハ

・試験本番
そういえば午後ノー勉だけど大丈夫？　→　配られた解答用紙を見たぼく「記述式マジ？初めて知った」

・結果
午前が70ちょっと、午後は80ちょっと

・所感
基本情報は英語と国語とプログラミングが出来ると合格できそう。
ぼくは英語ができないので英語の勉強をしました（はい？）
フォールトアボイダンスはフォールトをアボイドするという意味だし、フェイルセーフはフェイルしてもセーフという意味だし、フォールトトレランスはフォールトにトレランスだという意味。（トレランスという単語を初めて知りました）
応用情報は、基本情報との難易度の違いがわからなかった。出題の幅が広い気がするので、自分が取る予定の選択問題の難易度が爆発してたら死ぬってことがあるらしい？　ぼくは全部の問題をざっと解いてから、高点数が取れそうな問題を選ぶ立ち回りをしたのでよくわかりませんが……

・いかがでしたか？
巷で言われてるよりもずっと簡単ですよ。
みなさんも受けましょう。
ちなみに恩恵はいまのところありません。

・先行研究
degwer.hatenablog.com
ぼくは前述の通りノー勉だと5割にすら届かなかったんですが、賢い人がやると合格点の6割に届くらしいです

2019-05-31

位数Xのアーベル群

出典

さて、この問題はopが可換群であればBITを用いて解くことができます
数列aの全ての要素、i, j, kが全て1以上X以下の値であることが分かっている時、任意の数列aに対しBITで解くことのできる演算opが全部で何通りあるか求めてください
— 31536000 (@CuriousFairy315) 2019年5月28日

この問題は「互いに区別できるX要素の集合に、アーベル群の構造を入れる方法は何通りあるか？」と同じになる。
これはOEISで検索すると見つかる。
A034382 - OEIS
項の数が少なく(注:この記事を書いた当時は20項しかなかった)、求め方も載っていないので、最初のツイートの制約で解くにはまだ足りない。導出を考えてみよう。

以下、単に群といえばアーベル群を指すものとする。
各記号を次の通り定める。
pは素数とする
$C_n$ 位数nの巡回群　 $\mathbb{Z}/n\mathbb{Z}$ のこと
S(X) 位数Xのアーベル群の同型類の集合(各同型類とその代表元を同一視する)　これはここだけの記号であり一般的なものではない
Aut(G) Gの自己同型群
GL(k,R) 可換環Rのk×k行列で可逆なものの集合
単に「元gの位数がn」といった場合は「ちょうどn」とは限らないものとする。（例えばC4の元は全て位数4である）

$X=\prod p_i^{e_i}$ とする。
S(X)から群を1つとりGとすると、それと同型になるようにするには、どの元にどれを割り当てるかX!通り。……ではなくX!/#Aut(G)通り。
例えば $C_3$ になるよう{0,1,2}に群構造を入れるとき、0を単位元にしたら、のこりの2つの元のどちらを1,2に割り当てても同じになるので。
ということで、求める答えは
$\displaystyle{ \sum_{G \in S(X)} \frac{X!}{\# Aut(G)} }$
ここで、a,bが互いに素ならS(ab)=S(a)×S(b)であること、
GとG'の位数が互いに素ならAut(G×G')=Aut(G)×Aut(G')であることから
$\displaystyle{ \sum_{G \in S(X)} \frac{X!}{\# Aut(G)}\\ =X!\prod_i \sum_{G \in S(p_i^{e_i})} \frac{1}{\# Aut(G)} }$
となる。
（これは $\sum_{i,j,k}a_i b_j c_k=(\sum_i a_i)(\sum_j b_j)(\sum_k c_k)$ と同じ式変形）
位数 $p^e$ の群は、必ず $\prod_i C_{p^{n_i}}^{k_i}$ の形で書くことができる。ここで $\{n_i\}$ は狭義単調増加するようとるものとする。
位数の条件から $e=\sum_i n_i k_i$ である。
巡回群は生成元の行き先さえ決めれば準同型が定まる。 $n_i$ の大きい方から考えることで、準同型が同型であるためには、 $C_{p^{n_i}}^{k_i}$ の生成元たちは送った先でもこの部分を生成していなければならないことがわかる。そのような割り当て方は
$\displaystyle{ \# GL(k_i,\mathbb{Z}/p^{n_i}\mathbb{Z})=p^{(n_i-1)k_i^2}\prod_{j=0}^{k_i-1} (p^{k_i}-p^{j}) }$ 通りである。
( $n_i=1$ のとき、 $\mathbb{Z}/p\mathbb{Z}$ が体であることから「k次元ベクトルk本を一次独立に選ぶ」ことと対応するので総積部分が得られる。 $n_i \geq 2$ のときは、"1の位"をそのように決め、"pの位"から上は自由に決められるので残りの部分が得られる)
また各iについて、送った先で $C_{p^{n_i}}^{k_i}$ にあたる部分以外は位数 $p^{n_i}$ であれば何を割り当ててもよいので、 $\displaystyle{ ( \prod_{j\neq i} p^{\min(n_i,n_j)k_j} )^{k_i} }$ 通りとなる。
（例えばC2×C4×C8のとき、C4の生成元はC2の部分は自由に取れ、C8の部分は位数4の元4つの中から自由にとれる）

よって以上より
$\displaystyle{ \# Aut(\prod_i C_{p^{n_i}}^{k_i}) =\prod_i \left( (p^{(n_i-1)k_i^2}\prod_{j=0}^{k_i-1} (p^{k_i}-p^{j})) ( \prod_{j\neq i} p^{\min(n_i,n_j)k_j} )^{k_i} \right) }$ 通りとなる。

ここで $\# S(p^e)$ はeの分割数に等しいことがわかるので、これは対して大きくならず、結局計算のボトルネックになるのはXの素因数分解部分であり、計算量は $O(\sqrt{X})$ である。

……と思ったのだが、プログラムを書いて計算すると、X=16のときがOEISの結果と一致しない。
上の計算式によれば
#Aut(C16)=8
#Aut(C2×C8)=(1×2)×(2×4)=16
#Aut(C4×C4)=(16×3×2)=96
#Aut(C2×C2×C4)=((3×2)×2^2)×(2^2×2)=192
#Aut(C2^4)=(15×14×12×8)=20160
であるから、求める答えは16!×(1/8+1/16+1/96+1/192+1/20160)=4250979532800となる。
一方、OEISでは4248755596800になっている。これは第4項が1/192ではなく1/196なら一致することがわかるが、プログラムを書いてみると、確かにAut(C2×C2×C4)=192であることが確かめられる。そもそも49の倍数になるわけないだろ、なんでだよ

ここまでの議論に誤りを見つけた方、または、OEISの記述が間違っていると確信が持てる方はご一報ください。
なお、上の議論に則って実装したコードはこれです。(mod 10^9+7で出力）
https://wandbox.org/permlink/9RLseRNaHU5CeQiU
自己同型群の位数が計算できてるので実装がバグってることはないはず

(2019/09/13追記)
mathoverflowできいた。
mathoverflow.net
OEISが間違っていた。（えー）
その結果、OEISにsugarknriの名前が載った（！！！？？？！？？！？！？）

2018-10-10

yukicoder No.712 赤旗

問題はこちら
No.712 赤旗 - yukicoder

'W'の数を数えるだけ。

s;
main(c){
	scanf("%*d%*d\n");
	while(c=getchar(),~c)if(c=='W')s++;
	printf("%d",s);
}

2行目以降は'\n'(10),'R'(82),'W'(87)のいずれかしかないので&1で識別可能。1行目の読み飛ばしとうまく噛み合わせる

w,s;main(c){for(;read(0,&c,1);w=c<11?:w)s+=c&w;printf("%d",s);}

63B

2018-10-10

yukicoder No.706 多眼生物の調査

問題はこちら
No.706 多眼生物の調査 - yukicoder

strlen(s)-2の最頻値を見るだけ

int a[1010];
int main(){
	int n;
	scanf("%d\n",&n);
	for(int i=0;i<n;i++){
		char s[1010];
		gets(s);
		a[strlen(s)-2]++;
	}
	
	int ans=2;
	for(int i=2;i<=1000;i++)if(a[i]>=a[ans])ans=i;
	printf("%d",ans);
}

ぐっと睨んでループをまとめる。文字列を読み込む変数をごまかしたり、端の処理をごまかしても通る

k,a[999],i;
main(m){
	for(;gets(&k)?++a[strlen(&k)]:i++<998?m=a[i]<a[m]?m:i:0;);
	printf("%d",m-2);
}

95B
これは996までしかチェックしていないので997や998が答えになるときWAになる

kの偶奇で分ける

iの偶奇で分ける

非再帰FFT

非再帰 FFT