2021-03-30

sum of multiplicative function

これはなに？

Min_25さんの以下の記事を自分が読みやすいように書き換えたものです。
min-25.hatenablog.com
Min_25さんによる実装はこちら
https://ideone.com/9F9amv

要求知識

・平方根の上下で分けるテク
・fenwick tree
・線形篩の"気持ち"を知っていると理解の助けになる箇所が一部ある

定義

定義域が正整数全体である複素数値関数 f が次を満たすとき、f は乗法的(multiplicative)であるという。
$\forall a,b\ ({\rm gcd} (a,b)=1 \Rightarrow f(ab)=f(a)f(b))$

定義から、f≠0 が multiplicativeなら f(1)=1 であることに注意。以下単に f とかいたとき、それは multiplicative であるとする。

以下では値域が複素数であることは使っていないので、演算がO(1)でできる環なら同様に計算できるはず。

記号の定義

・ $\displaystyle{ F_{{\rm prime}}(x):=\sum_{p\leq x,p:{\rm prime}}f(p)}$ 　素数のとこだけの和
・ $F(x):=\sum_{i=1}^{x} f(i)$ 　全部の和
・ ${\rm lpf}(x)$ ： x の最小素因子。 ${\rm lpf}(1)=\infty$
・ $p_k$ ： k番目の素数。1-based index
・ $\pi(x)$ ： x以下の素数の個数
・ $f_k(x):= \begin{cases} f(x) & ({\rm lpf}(x) \geq p_k) \\ 0 & (\text{otherwise}) \end{cases}$ 　最小素因子がでかいとこだけ生き残るやつ
・ $F_k(x):=\sum_{i=1}^{x} f_k(i)$ 　最小素因子がでかいとこだけの和。 $F_1=F$ に注意。
・ $H_N:=\{\lfloor \frac{N}{1}\rfloor, \lfloor \frac{N}{2}\rfloor,\ldots,\lfloor \frac{N}{\lfloor\sqrt{N}\rfloor}\rfloor\}$
・ $L_N:=\{1,2,\ldots,\sqrt{N}\}$
$H_N,L_N$ については N が明らかな時には単にH,Lと表す。

何ができるか

任意の素数べき $q=p^e$ に対し $f(q)$ は $q,p,e$ から $O(1)$ で計算可能であるとする。
全ての $x\in H_N\cup L_N$ に対して $F_{{\rm prime}}(x)$ が求まっているとき、全ての $x\in H_N\cup L_N$ に対する $F(x)$ を時間計算量 $O(N^{2/3})$ 、空間計算量 $O(N^{1/2})$ で求めることができる。

Q. 元記事に比べて仮定が強すぎない？（元記事では $F_{{\rm prime}}(x)$ の計算からやっているはず）
A. $F_{{\rm prime}}(x)$ の計算には本編と同じようなテクを使うので、説明の順番的に後回し

$F_k$ という概念を用意したことからお察しのとおり、kが大きいところから順に $F_k$ を計算していって、最終的に $F_1$ にたどり着きます。
前処理として $\sqrt{N}$ 以下の素数を列挙しておきます。
本編は3つのパートからなります。

お断り：この記事は気持ちを理解することを優先しているため、境界条件については適当です。正確な部分は元記事を見るか自分でよく考えて実装しましょう。

$x\in H\cup L$ に対して $F_{\pi(N^{1/3})+1}(x)$ を求める

要するに、 $N^{1/3}$ より大きな素因子しか持たない部分について計算しましょうという話です。 $k:=\pi(N^{1/3})+1$ 、 $p:=p_k$ とします。
${\rm lpf}(i) \geq p$ を満たすような i は「1」「素数」「素数の2乗」「相異なる2つの素数の積」のいずれかです。ということで、 $F_k$ は $F_{{\rm prime}}$ を使うことで簡単に計算できます。
具体的には x の値に応じて次の通りです。
・ $1\leq x < p$ のとき
$i\leq x$ の範囲には ${\rm lpf}(i)\geq p$ となる i は「1」しか存在しません。
よって $F_k(x)=1$ です。
・ $p\leq x < p^2$ のとき
$i\leq x$ の範囲には ${\rm lpf}(i)\geq p$ となる i は「1」「素数」しか存在しません。
よって $F_k(x)=1+(F_{{\rm prime}}(x)-F_{{\rm prime}}(p-1))$ です。
・ $p^2\leq x \leq N$ のとき
$i\leq x$ の範囲には ${\rm lpf}(i)\geq p$ となる i は「1」「素数」「素数の2乗」「相異なる2つの素数の積」しか存在しません。
よって $F_k(x)=1+(F_{{\rm prime}}(x)-F_{{\rm prime}}(p-1)) + \sum_{p\leq p'\leq \sqrt{x}}f(p'^2) + \sum_{p\leq p'\leq \sqrt{x}}f(p')(F_{{\rm prime}}(x/p')-F_{{\rm prime}}(p'))$ です。
第4項は「p'×(p'よりでかい素数p'')」の形で表される数について、p''で和を取ったものです。

計算量解析
上2つのケースは各xに対しO(1)です。そもそもxの候補が全体で $O(\sqrt{N})$ 個しかないので、この部分は合わせて $O(\sqrt{N})$ です。
3つ目のケースは各xに対し $O(\pi(\sqrt{x}))$ です。xとして $\lfloor \frac{N}{1}\rfloor, \lfloor \frac{N}{2}\rfloor,\ldots,\lfloor \frac{N}{\lfloor\sqrt[3]{N}\rfloor}\rfloor$ を計算するので、
$O(\sum_{i=1}^{N^{1/3}}\pi(\sqrt{N/i})) = O(\int_1^{N^{1/3}} (N/i)^{1/2} / \log (N^{2/3}) di) = O(N^{2/3}/\log N)$ です。
空間計算量は、 $F_k$ の値の保持以外に領域を使っていないため、明らかに $O(N^{1/2})$ です。

$x\in H\cup L$ に対して $F_{\pi(N^{1/6})+1}(x)$ を求める

ここが一番難しい。
さっきは $N^{1/3}$ まででしたが今度は $N^{1/6}$ までやります。 $F_k$ から $F_{k-1}$ を作るというステップを順番に刻んでいて $N^{1/6}$ まで行きます。更新はin-placeにやることにします。
説明のためだけに新たな記号を用意します（えぇ……）。
$x\in H\cup L$ に対し、 $x':=\max\{y \mid y\in H\cup L, y< x\}$ とします。 $H\cup L$ に入ってるやつで x の次に小さい値ですね。ただし、 $x=1$ のとき $x'=0$ とします。
$x$ が $N^{2/3}$ 以上か以下かで更新の仕方を変えるんですが、そもそも $N^{2/3}$ 以上か以下かで値の持ち方も変えます。
つまり、今までは、各 $x\in H\cup L$ に対して $F_k(x)$ を持っていましたがそれをやめて、 $x< N^{2/3}$ のときには、一旦 $F_k(x)-F_k(x')$ を持つことにします。 $x\geq N^{2/3}$ のときは今まで通り $F_k(x)$ を持ちます。
そうして、 $x< N^{2/3}$ の部分をfenwick treeで管理します。これで今まで通り $F_k(x)$ を取ってこれます。
要するに、「点更新クエリが飛んでくるなら累積和じゃなくてfenwick treeかなーやっぱりｗ」です。

更新の仕方を考えましょう
・ $N^{2/3}\leq x \leq N$ のとき
「もらうDP」の気持ちで更新します。
$F_k(x)=\sum_{e=0}^{\infty} f(p_k^e) F_{k+1}(x/p_k^e)$ を使って愚直に計算するだけです。
$x/p_k^e< N^{2/3}$ のときの $F_{k+1}(x/p_k^e)$ の値はfenwick treeにrange sumクエリを投げて取ってくることになるのでlogがつくことに気をつけましょう。
・ $1\leq x < N^{2/3}$ のとき
「配るDP」の気持ちで更新します。
まずは x=1から初めて $xp_k,xp_k^2,\ldots$ に配ります。そうしたら、配った先の地点 $y$ で、今度は更に $yp_{k+1},yp_{k+1}^2,\ldots,yp_{k+2},yp_{k+2}^2,\ldots,$ に配ります。線形篩と同じですね。
配る時にはfenwick treeにpoint addクエリを投げることになるのでlogがつくことに気をつけましょう。
$N^{2/3}$ 以上の範囲は上述の通り別に計算しているので、配る範囲は $N^{2/3}$ 未満までであることにも気をつけましょう。
配り元としては $N^{2/3}/ p_k < N^{1/2}$ までを考えればいいので、 $f_{k+1}(x)=F_{k+1}(x)-F_{k+1}(x')$ の値が問題なくとれることにも注意しましょう。
初手で配り先は $p_k$ 以上になるので、2個目以降の素数としては、 $N^{2/3}/p_k\leq N^{1/2}$ 以下のみを考慮すれば良く、事前に列挙した範囲で収まることにも気をつけましょう。
気をつけることがいっぱい。

計算量解析
ここも難しい
・ $N^{2/3}\leq x \leq N$ の範囲について
面倒なのでfenwick treeにクエリを投げなくていいところも投げてしまうことにすると、クエリを投げる回数の総和は
$O(\sum_{N^{1/6}< p\leq N^{1/3}, p:{\rm prime}}\sum_{i=1}^{N^{1/3}}\log_p(N/i))$ となります。
$\sum_{i=1}^{N^{1/3}}\log_p(N/i) =O(N^{1/3}\log N/\log p)=O(N^{1/3})$ に気をつけて変形すると $O(N^{1/3}/\log N \cdot N^{1/3})$ となるので、クエリ1回に $O(\log N)$ かかることから、全体の計算量は $O(N^{2/3})$ です。
・ $1\leq x < N^{2/3}$ の範囲について
配る先がいくつあるかがクエリを投げる回数になります。
配る先は、lpf が $N^{1/6}$ 以上であるような $N^{2/3}$ 未満の数全てです。そのような数は $O(N^{2/3}/\log N)$ 個らしいので*1、クエリ1回に $O(\log N)$ かかることから、全体の計算量は $O(N^{2/3})$ です。
$F_k$ の計算時には $F_k,F_{k+1}$ の値の一部をfenwick treeで持っていますが、定数倍の違いしかないので空間計算量は $O(N^{1/2})$ です。

$x\in H\cup L$ に対して $F(x)$ を求める

$F_k(x)=\sum_{e=0}^{\infty} f(p_k^e) F_{k+1}(x/p_k^e)$ を使って愚直に計算するだけです。やるだけ。

計算量解析
k を1つ進めるために、各 xに対し、 $\log_{p_k}(x)$ 回の計算が必要なので、全体で
$\sum_{p\leq N^{1/6}, p:{\rm prime}} \sum_{i=1}^{N^{1/2}}(\log_p(N/i)+\log_p(i))$ くらいの更新が必要で
$O(\sum_{2\leq p\leq N^{1/6}} \sum_{i=1}^{N^{1/2}}\log_p(N))=O(N^{1/2}\log N \sum_{2 \leq p\leq N^{1/6}} 1/\log p) =O(N^{2/3})$ になる。
$F_k$ の計算時には $F_k,F_{k+1}$ 以外の値を保持していないので空間計算量は $O(N^{1/2})$ です。

結局 $F_{{\rm prime}}$ は？

TODO：書く

2<=p<=N^{1/6}は (Lucy) DP で O(N^{2/3}/log(N))。
N^{1/6}<p<=N^{1/3} は、 N^{2/3} 未満は Fenwick で愚直更新で O(N^{2/3}/log(N)*log(N))、N^{2/3} 以上は DP + Fenwick で O(N^{1/3}/log(N)*N^{1/3}*log(N))。
N^{1/3}<p<=N^{1/2} は、Lucy DP で O(N^{2/3}) (N^{2/3}/log(N)) (?) 想定でした。
— Min_25 (@min_25_) 2020年4月7日

Lucy DPというのは $O(N^{3/4})$ prime counting でおなじみのやつです（多分）

*1:Buchstab function - Wikipedia

2021-02-11

A Sequence of Permutations

a[i]=a[i-1]-a[i-2]の特性方程式が1の6乗根だからmod6で見れる、なるほどなぁという気持ちに
— お気持ち表明を控える (@tk0_math) 2021年2月10日

ほんまか？

$a_{n}=a_{n-1}-a_{n-2}$ は $a_{n}=a_{n-6}$ を満たすので、mod 6で見たいという気持ちで式変形をすると
$a_n\\ = a_{n-1} a_{n-2}^{-1}\\ = a_{n-2} a_{n-3}^{-1} a_{n-2}^{-1} \\ = (a_{n-4} a_{n-5}^{-1} a_{n-4}^{-1}) (a_{n-5} a_{n-6}^{-1} a_{n-5}^{-1})^{-1} (a_{n-4} a_{n-5}^{-1} a_{n-4}^{-1})^{-1}\\ =(a_{n-4} a_{n-5}^{-1} a_{n-4}^{-1} a_{n-5}) a_{n-6} (\ldots)^{-1}$
となって
$a_{n-4} a_{n-5}^{-1} a_{n-4}^{-1} a_{n-5}\\ =(a_{n-5} a_{n-6}^{-1}) a_{n-5}^{-1} (a_{n-5} a_{n-6}^{-1})^{-1} a_{n-5}\\ = a_{n-5} a_{n-6}^{-1} a_{n-5}^{-1} a_{n-6}\\ = a_1 a_0^{-1} a_1^{-1} a_0 =:X$
として
$a_n = X a_{n-6} X^{-1}$
を得る。よって $O(N\log K)$ で解けた。

ところでなんでこれうまくいくんすかね。

例えば $a_{n}=a_{n-1}a_{n-2}^{-1}a_{n-3}a_{n-4}^{-1}$ で同じ問題を考えると、同様にある定数 X により $a_{n}=X a_{n-10} X^{-1}$ と書けるが、途中の計算は違っていて、例えば $a_{n}=X(n) a_{n-5}^{-1} X(n)^{-1}$ とは書けない。

証明or反例を募集しています

2020-12-21

yukicoder 1322

お断り：提出してないので間違ってるかも

O(N^0.75) prime countingを完全に理解していてこれが解けないのは端的に言ってカスなんだよね

・そもそもφはmultiplicativeなので、当然素因子ごとにわけて考えたくなる
・dp[i][j]=#{n | nの全ての素因子はPi以下、φ(n)=j}としてDPできる
遷移はdp[i-1][j]からdp[i][j*(Pi-1)*Pi^k]に配れば良い
・第2添字には乗算しか登場しないので、「あと何倍したらNを超えるか」しか情報がいらない。平方根で分けるテクで状態数がO(√N)になる
・in-placeに更新することにすると、O(N^0.75) prime counting とほとんど同じ計算量解析により、p*(p-1)≦Nを満たす素数p=O(√N)までO(N^0.75)で計算できることがわかる
・それ以降の素数は高々1個しか含まれないので、各pに対してφ(x)≦N/(p-1)となるxの個数がわかればよい（そしてこれは既に計算済み）
・N/(p-1)は小さいので、pごとではなく、N/(p-1)ごとにまとめて計算すれば高々O(√N)。N/(p-1)=kとなるようなpの個数はO(N^0.75) prime countingを実行すると副産物として手に入っている。（境界が1個ずれてる分はいちいちチェックすれば良い）
・よって解けた

2020-10-28

純粋培養競技プログラマが就職して１年

ポエム

■スペック
国立大理系院卒
就職まで競プロ以外でのプログラミング経験なし
今はAtCoder橙
社会性破滅・就活破滅
英語壊滅
IT中小企業にどうにか潜り込む

就活編はこちら
sugarknri.hatenablog.com

■1年経過時点での環境
pythonで機械学習やってる
git使ってない(バージョン管理はSVN)
Linuxわからない
Docker？ AWS？なんすかそれ

上の方の役職にいる人や意欲的な人は、機械学習の最新の論文を追ったり、Kaggleで遊んでいたりするらしい。俺よりすごいというのはわかる。
そうでない大多数の人たちはどうなんでしょう

周りの人たちと比べてコーディング速度が3倍くらい違うので異世界転生感があってよい。
そういう環境なので向上心は持っておかないとそのまま怠惰で破滅しそう（持ってますか？）

■仕事
最初の数カ月は研修。残りはPoC案件2つを主にやっていたら1年終わった。ソースコード書く人は俺と先輩1,2人だけ、みたいな状況なので、チーム開発とかいうのはよくわからず。PEP8に従ったコーディングをしろと言われたのでそれは守ってる。
コードレビューも最初は熱心な先輩に丁寧に見てもらっていたけど、それ以外の人からは……。ところで、前任者から引き継いだ納品済みのコードが無限にバグっていました。大丈夫ですかこれ？（自明なバグから非自明なバグまで大量に指摘して、私の評価が大変に上がったのでありがたいといえばありがたい（いいえ））

■今後
どうするんでしょうね

2020-10-18

指数型母関数

テクニック

この記事は 37zigenさんによる以下の記事をもとに、お気持ちパートを自分にとってわかりやすいようにメモしたものです。
真面目な解説は書いていないので、そういうのを読みたい人は元記事の方をみてください。
37zigen.com

・物を数える

気持ち：「互いに区別できるN要素のものがあります。いい感じの条件を満たすような構造の個数を考えます」というタイプの問題を考えます
構造ってなんやねん。なんらかのなんかです。例をみて

自明１：互いに区別できるN要素のものからなる集合はいくつありますか？
1個。それはね、そう

自明２：互いに区別できるN要素のものからなる順列はいくつありますか？
N!個。はい

自明３：互いに区別できるN個の箱に互いに区別できるK個のボールを入れる方法はいくつありますか？
N^K　はい

自明４：互いに区別できるN要素のものからなる有向サイクルはいくつありますか？
(N-1)!　円順列って知ってる？

非自明：互いに区別できるN要素のものからなる木はいくつありますか？
実はN^(N-2)なんですね

こういう各問題について、Nに対する答えをa_nとしたとき、 $\sum(a_n \frac{x^n}{n!})$ というものを考えます。
これ{a_n}の指数型母関数(EGF)っていうらしいですよ。

以下では構造Aの指数型母関数をそのままA(x)と書きます。

・母関数の積

互いに区別できるN要素のものがあります。そのうちいくつかを使って構造Aのものを、残りの全部を使って構造Bのものを作ります。そのような方法は何通りありますか。

適当に計算をすると、こいつの指数型母関数はA(x)B(x)になることが示せます。
2種類以上の構造が登場するときに役に立ちます。

問題：互いに区別できるN個の要素を、aの倍数個の要素からなる集合とbの倍数個の要素からなる集合の2つに分ける方法はいくつありますか？　ただしa≠b
「互いに区別できるN要素のものからなる、aの倍数個の要素の集合はいくつありますか？」は、Nがaの倍数のとき1、そうでないとき0なのでEGFもそんな感じで求まる。
それぞれA(x)、B(x)としてA(x)B(x)が答え

他にも構造Bを「集合」にすると包除もできます。

問題：N要素の撹拌順列はいくつありますか？
撹拌順列のEGFをA(x)、集合のEGFをB(x)、順列のEGFをC(x)とすると、A(x)B(x)=C(x)です。
B(x)、C(x)の正体はわかっているので、A(x)は頑張ると計算できます。

問題：互いに区別できるK個のボールを互いに区別できるN個の箱に入れる方法はいくつありますか？　ただし、どの箱も空になってはいけない
Kを固定して、「i個の箱に入れる。空はダメ」というもののEGFをA(x)、集合のEGFをB(x)、「i個の箱に入れる。空でもOK」というもののEGFをC(x)とすると、A(x)B(x)=C(x)です。
B(x)、C(x)の正体はわかっているので、A(x)は頑張ると計算できます。

・母関数の合成

構造Aをなす各要素を、構造Bをなしているもので置き換えます。全体でN要素になるような方法は何通りありますか。

構造を「代入」しているように、実はEGFの方も代入してA(B(x))になるんですね。なにそれすごい

問題：互いに区別できるN頂点の連結functional graphの個数を求めよ
連結functional graphは、有向サイクルの各頂点を根付き木に置き換えたものです。
有向サイクルのEGFをA(x)、根付き木のEGFをB(x)とすると、求めるもののEGFはA(B(x))です。

問題：互いに区別できるN頂点の森の個数を求めよ
森は、集合の各要素を木に置き換えたものです。
集合のEGFをA(x)、木のEGFをB(x)とすると、求めるもののEGFはA(B(x))です。

問題：互いに区別できるN頂点の単純連結無向グラフの個数を求めよ
連結性の条件を落とすと、明らかに2^(N(N-1)/2)個です。EGFをC(x)とします。
連結とは限らない単純無向グラフは、集合の各要素を単純連結無向グラフに置き換えたものです。
集合のEGFをA(x)、求めるもののEGFをB(x)とすると、A(B(x))=C(x)です。
A(x)、C(x)の正体はわかっているので、B(x)は頑張ると計算できます。

・余談

重みつきで和をとりたかったら、当然係数に重みをつければいいわけですね。

問題：互いに区別できるN頂点の森について、「2^(連結成分の個数)」の和を求めよ
集合のEGFの代わりに a_n=2^n に関するEGF をA(x)、木のEGFをB(x)として、A(B(x))です。

問題：互いに区別できるN頂点の森について、「各連結成分の大きさの積」の和を求めよ
木の大きさをどうやってEGFに載せようかという話になるわけです、a_nの代わりにn*a_nを使いたいという話なので、木のEGF B(x)のかわりにx(d/dx)B(x) を使えばよくて、終わりです。

要求知識

定義

記号の定義

何ができるか

に対して を求める

に対して を求める

に対して を求める

結局 は？

$x\in H\cup L$ に対して $F_{\pi(N^{1/3})+1}(x)$ を求める

$x\in H\cup L$ に対して $F_{\pi(N^{1/6})+1}(x)$ を求める

$x\in H\cup L$ に対して $F(x)$ を求める

結局 $F_{{\rm prime}}$ は？