2022-03-24

project Euler 1-50 解説

間違ってたりより良い解法があったら教えて

計算量は四則計算の回数。桁が増えすぎると壊れる

1
包除原理
これくらい手で解け

2
愚直で解ける
項数を先に求めて行列累乗すればO(loglogN)
一般項を出して適当なmodを取ってCRTすると実質O(1)

3
素因数分解の計算量よくわからん
変な数を素因数分解したくなったら、そのへんのサイトに投げてみよう
https://www.alpertron.com.ar/ECM.HTM
http://factordb.com

4
乗数と被乗数を全探索計算量は謎
こういうのは、

int width=1e4;
for(int U=1e6;;U-=width){
  int L=U-width;
  for(int x=999;x*x>=L;x--){
    for(int y=min(U/x,x);x*y>=L;y--){
      //hoge
    }
  }
}

みたいにやるのが定石。意外と忘れやすい

5
lcm取るだけ O(NlogN)
素数の指数のmaxを考えると O(NloglogN)

6
1乗和・2乗和・3乗和のclosed-formを知っていますか？
ファウルハーバーの公式 - Wikipedia

7
篩をするとO(NlogN poly(loglogN))とかになるはず
O(N^(2/3))素数カウント+二分探索+区間篩でもっと早くなったりするのかな

8
愚直 O(NK)
0に注意して尺取するとO(N)

9
愚直 O(N^2)
ピタゴラス数の列挙式を使うと O(√N)
ピタゴラスの定理 - Wikipedia
正解者掲示板で1000ふぁぼついてる投稿が嘘解法で草。何が嘘かわからない人は1000じゃなくて36で解いてみて

10
愚直 O(NloglogN)
おまたせ　いつもの　親の顔より見た Lucy DP　O(N^(3/4)/logN)
実はO(N^(2/3))とかで出来るらしいです

11
愚直 O(NNK)
0に注意して尺取するとO(NN)
この前のABCで見た
C - Connect 6

12
愚直にやると約数カウントが大変そう
osa_k法で高速素因数分解すると早そう
計算量はわからん

13
やるだけ
ところで、例えば上n桁で概算すると、その下2桁は(100個の数を足すので)50くらいの誤差を持つことが期待される。ということは上13桁くらいで概算すればよさそう
この手の嘘計算量削減や、最後の1桁総当り作戦もオイラーでは重要なテクニックの一つ（？）

14
メモ化再帰

15
二項係数

16
多倍長整数

17
桁ごとに分けるテク
この前のABCで見た
C - Comma

18
愚直 O(2^N)
DP O(N^2)

19
ツェラーの公式
日付型ライブラリ
だいたい1/7くらいでしょ！ｗ→1200/7の前後を総当り

20
多倍長整数

21
エラトステネスみたいに約数和の配列を作ってO(NlogN)
線形篩と同様にすればO(N)

22
ファイル入出力

23
21と同様に過剰数を列挙して愚直でO(N^2)
FFTでO(NlogN)

24
貪欲でO(N^2)
セグ木で持てばO(NlogN)（実際には桁数が膨大になって四則演算がO(1)にならないが……）

25
ビネの公式

26
ラグランジュの定理
ラグランジュの定理 (群論) - Wikipedia

27
枝刈り全探索
最悪でもn=bで終わるので、2e6+αくらいまで篩えばOK

28
階差数列

29
多倍長整数 or logとってset O(N^2 logN)
実はN=10^16とかでも解けて466がそれ（え？）

30
左辺ではなく右辺を全探索

31
DP
制約強化版がここにある
No.137 貯金箱の焦り - yukicoder

32
全探索

33
全探索
実は手でも解けるらしい

34
30と同じようにやると見せかけて、桁数が多いので上下に分けるとより高速になる（こういうのもmeet in middleっていうのか？）

35
1e6まで篩っておいて4^kを全探索

36
回文は上半分を決めると決まる
この前のABCで見た
E - (∀x∀)
2進数で回文かどうかの判定は、bit reverseしてtrailing zeroを削除するとビット演算だけでできる

37
素数を上に一桁伸ばすのと下に一桁伸ばすの、条件が厳しいのはどーっちだ？
条件が厳しいということは候補が少ないということです。そっちを全探索しましょう

38
これくらいなら紙と鉛筆で解きませんか？

39
9で解いた
ピタゴラス数の列挙式ではなく、生成行列を使うと、gcdのlogがつかない分早い
Pythagorean Triple -- from Wolfram MathWorld

40
再掲
C - Comma

41
一般に、2,3,5,7あたりの自明な条件は、最初に考察で弾くほうが良い
残りは全探索

42
set でクエリあたりO(log)＋前計算
平方根を計算すればクエリあたりO(1)

43
条件が厳しいということは候補が少ないということです（再）
手でも解けるそうです

44
実は適当にやっても通ってしまいますが、その解法が正しいことは誰も証明できていません（俺は思いつかないし、正解者掲示板にも書いていない）
計算量の保証がある解法を思いつくのはめちゃくちゃ難しい
Project Euler 44（1） - inamori’s diary

45
ペル方程式

46
答えが小さいことを祈る

47
答えが小さいことを祈る2

48
modの基礎
二分累乗するまでもない

49
並び替えを試す前に map< int, vector< int >> に分解しておくと早い

50
累積和＋尺取
尺取でバグらせたくなければ、区間DPみたいにlenでループするのもあり

for(int len=U;len>0;len--){
  for(int l=0,r=len;sum[r]-sum[l]<1e6;l++,r++){
    //hoge
  }
}

続き
sugarknri.hatenablog.com

2022-03-11

segment tree beats 覚書

ポエム

これを読め
rsm9.hatenablog.com
～終わり～

以下ポエム
ひとなのでさんによるこの記事を読むまでbeatsの理解が困難だったが、この記事を読むと全てを3秒で理解できた。
これ以前に存在していた記事で理解が困難だったのは、まずそもそもbeatsというデータ構造の概念の部分の話・実装に関する話・計算量解析に関する話がいっしょくたにされていることが多かったため。
beatsの本質は「mapping(f,x)が計算できないなら、op(left,right)で計算する」というだけ。冒頭の記事ではこれを一番最初に説明してくれており、それでようやくbeatsというものが何をするデータ構造なのか理解できた。この段階では計算量解析には触れずに「もしそういうことができるなら、これでうまくいくね」という話し方をしている。非常にわかりやすい。

作用 f(x) の計算に失敗してしまった場合，「とりあえず f を x の子に伝播させて子の計算を先に行い，子に関する計算結果を用いてボトムアップに x′ を再構築する」という方法で問題を回避できそうです．ここで計算の失敗回数が多すぎる（例えば列の長さ N やクエリ個数 Q に関する二乗のオーダーになってしまう）と計算量が悪化してしまいますが， S や F の構造などに由来する何らかの理由で「Segtree 上の全頂点における失敗回数の合計」に良い上界（例えば O((N+Q)log2⁡N)）が与えられるならば，計算量の本質的な悪化が回避でき，高速なクエリ処理が可能となります．
atcoder::lazy_segtree に1行書き足すだけの抽象化 Segment Tree Beats - ひとなので

実装パートも、beatsが何なのかを冒頭のように理解していれば1行書き足すだけなのは自明。非常にわかりやすい。
計算量解析パートの理解が困難だったのも、chmin/sumのような複雑な例を最初に出されることが多かったため。ひとなのでさんの記事はこの点でも完璧で、1番目の例は、gcdが本質的にchminであることから「failするたびその区間に属する値の"max"が真に小さくなるため(特に1回目でクエリの値まで落ちる)」と簡単に理解できる。次の例も「(Z/2Z)^64という束の"高さ"が64であり、failするたびにその区間に属する値の"min"と"max"の"高さ"の差が真に小さくなるため」と理解できる（ハッセダイアグラムを頭に思い浮かべながらchand/chorクエリの効果をイメージできますか？ぼくはできます）。わかりやすい。
データの持ち方パートについても、これらはともに「chminがうまく走るか見たいから区間の最大値もっとこ（最大値＜chminクエリの値なら無視できるので）」というだけなのでかなり自然。
それに比べてchmin/sumはめちゃくちゃ複雑。なぜかというと「failするたびにその区間に属する値のmaxが真に小さくなる」を言っても、それは1e9回起こりうるため。もっと小さくなる何かを持ちてえなあという気持ちになり、second maxを持って区間内の値の種類数を減らそうと思える（は？）。計算量解析は難しい（「failするたびにその区間に属する値の種類数が減るため」ではO(N^2)しか言えないので）（は？）　こんなもんを初心者に見せるな
以上です。

2022-02-14

a+b=a^b+2(a&b)の一般化

bit毎の排他的論理和を $\oplus$ 、論理積を $\odot$ と置くと、
任意の非負整数 $a, b$ に対して
$a+b=a\oplus b+2(a \odot b)$
が成り立ちます。
これは
$0\oplus 0+2(1\odot 1)=0+0=0\\ 1\oplus 1+2(1\odot 1)=0+2=2\\ 0\oplus 1+2(0\odot 1)=1+0=1$
を観察すれば成り立つことが示せます。
この等式は競技プログラミングで非常によく使われるテクニックの一つです。
Yukiちゃんはこの等式を一般化して、
$\displaystyle x_1+x_2+\cdots+x_N=\sum_{k=1}^{N}a_k\sigma_k$

の形で表せないか気になりました。ここで $\sigma_k$ は $\oplus ,\odot$ を演算とする $x_1,x_2,\ldots,x_N$ の $k$ 次対称式です。
例えば、 $N=3$ のときは、
$\displaystyle x_1+x_2+x_3=a_1(x_1\oplus x_2\oplus x_3)+a_2(x_1\odot x_2\oplus x_2\odot x_3\oplus x_3\odot x_1)+a_3(x_1\odot x_2\odot x_3)$
です。
実は帰納法により、
$\displaystyle x_1+x_2+\cdots+x_N=\sigma_1+2\sigma_2+4\sigma_4+8\sigma_8+16\sigma_{16}+\cdots$
であることが示せます。

37zigen.com

2021-07-05

線形漸化式と微分方程式

テクニック

微分方程式を微分作用素で解くやつ知ってる？

f''-3f'+2f=0　の解は？
微分作用素をDとして、与えられた方程式は(D-1)(D-2)f=0
ここで一般に(D-k)f=0の解はC*e^kxなので、f(x)=C1*e^x+C2*e^2x

全く同じようにして線形漸化式の一般項が出るんですよね

a[n+2]-3a[n+1]+2a[n]=0　の一般項は？
左シフト作用素をLとして、与えられた方程式は(L-1)(L-2)a=0 (※)
ここで一般に(L-k)a=0の解はC*k^nなので、a[n]=C1*1^n+C2*2^n

(※)微分方程式を解く時には、測度0の集合上での違いを無視する(ように同値類を定めて割る)のと同様に、高々有限個の項の違いは無視している

ね？簡単でしょ
ついでにもう1問

f'''-6f''+12f'-8f=0　の解は？
(D-2)^3 f = 0 なので、f(x)=(C2*x^2+C1*x+C0)e^2x

a[n+3]-6a[n+2]+12a[n+1]-8a[n]=0　の一般項は？
(L-2)^3 a =0 なので、a[n]=(C2*n^2+C1*n+C0)2^n

これでもう任意の線形漸化式の一般項を出せるね。めでたしめでたし

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

要求知識

定義

記号の定義

何ができるか

に対して を求める

に対して を求める

に対して を求める

結局 は？

$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}}$ は？