背景：蟻本をやります - 競プロはじめました
テキスト：
![プログラミングコンテストチャレンジブック [第2版] ～問題解決のアルゴリズム活用力とコーディングテクニックを鍛える～](https://m.media-amazon.com/images/I/41oruV+aJIL._SL500_.jpg "プログラミングコンテストチャレンジブック [第2版] ～問題解決のアルゴリズム活用力とコーディングテクニックを鍛える～")
ユークリッドの互除法
ユークリッドの互除法については，次で解説しています．代数の本を買えば，大体載っていると思います．Pythonでは，math.gcd(x, y)で最大公約数が得られます（Pythonで最大公約数と最小公倍数を算出・取得 | note.nkmk.me）．本文と同じように再帰関数で書くと
def gcd(a, b): if b==0: return a else: return gcd(b, a%b)
です（elseはなくてもよい（最初のreturnで抜ければ実行されないため））．
線分上の格子点の個数
【入力】1 11
5 3
【出力】
3
【コード】
import math x1, y1 = map(int, input().split()) x2, y2 = map(int, input().split()) if x1==x2 and y1==y2: exit(print('0')) else: print(math.gcd(abs(x1-x2), abs(y1-y2)) -1 )
【ポイント】
本文のように図を書いてみる（※）のがわかりやすいと思いますが，若干補足です．
※本文の方法：2点間の傾きを表す大きな三角形を，辺の長さが整数の相似な三角形に分割する方法を考える．
【式を使った考え方】
2点$(x_{1},y_{1})$，$(x_{2},y_{2})$を通る直線の傾きは
\frac{y_{2}-y_{1}}{x_{2}-x_{1}}
\end{aligned}
これを約分して既約分数にした結果を$q/p$とすれば，$x_{1}$から$\pm x$方向に$p$進んだ直線上の点，あるいは$y_{1}$から$\pm y$方向に$q$進んだ直線上の点が一番近い格子点になります．次の格子点も同じ間隔で存在します．
よって，$|\frac{x_{2}-x_{1}}{p}| -1$または$|\frac{y_{2}-y_{1}}{q}| -1$が2点間に存在する格子点の個数です．ここで，$|(y_{2}-y_{1})/q| = |(x_{2}-x_{1})/p|$は，（既約分数の定義から）$|y_{2}-y_{1}|$と$|x_{2}-x_{1}|$の最大公約数です．
【補足 (もう少し厳密に)】
直線の方程式は
y-y_{1}
=\frac{y_{2}-y_{1}}{x_{2}-x_{1}} (x-x_{1})
\end{aligned}
(y-y_{1})p = (x-x_{1})q
\end{aligned}
$x,y$が格子点上にあることと，$x-x_{1}, y-y_{1}$が整数であることは同値です．したがって，直線上の点$(x,y)$について，①「$x,y$が格子点上にある」とき，上式で$p,q$は互いに素であることから，②「$p$は$(x-x_{1})$の約数で，$q$は$(y-y_{1})$の約数」です．逆に，②が成り立つとき，$x-x_{1}, y-y_{1}$は整数なので，①が成り立ちます．
よって，$(x-x_{1})/p = (y-y_{1})/q$が整数となる$(x,y)$が格子点を意味します．
拡張ユークリッドの互除法 - 双六
これも，大概の代数の本でユークリッドの互除法と一緒に解説されていると思います．ちなみに，$ax+by=\gcd(a,b)$なる$x,y$の存在を示す際，与えられた$a,b$を予め最大公約数で割って考えれば，互いに素（i.e. $\gcd(a,b)=1$）な$a,b$に対し$ax+by=1$を満たす$x,y$の存在を示せば十分です．Wikipediaの行列演算による解説がわかりやすいですね：https://ja.wikipedia.org/wiki/ユークリッドの互除法
本と同じようにx, yを変数にすると上手くいかないなぁ，と思っていましたが，次のようにすればよいのですね：Pythonで理解する蟻本「2-6 双六」(p.108) - クルトンのプログラミング教室
「拡張ユークリッドの互除法」に関するコードは以下（双六の問題は，出力方法を少し書き換えるだけなので略します）．
【入力】
4 11
【出力】
4*3+11*-1=1
【コード】
a, b = map(int, input().split()) def extgcd(a, b): ''' ax+by=gcd(a, b)となるx,yに更新 返り値はgcd(a, b), x, y ''' d = a if b!=0: # (1) by + (a%b)x = gcd(b, a%b) の解x, yに更新 d, y, x = extgcd(b, a%b) # (1) <=> (2) ax + b(y - (a//b)x) = gcd(a, b) y -= (a//b) * x else: x, y = 1, 0 return d, x, y d, x, y = extgcd(a, b) print('{}*{}+{}*{}={}'.format(a, x, b, y, d))
素数
素数判定
iの初期値が2であることに注意．最後のreturn n != 1がスマートですね．$2, 3,... \leq\sqrt{n}$のどの自然数でも割り切れなかった$n$に対し，それが1ならFalse, 1以外ならTrueを返すことができます（本来は，2以上の値が入力されているかを，最初に判定すべきなのかもしれませんが）．
def is_prime(n): ''' 素数ならTrue, そうでないならFalse ''' i = 2 while i**2 <= n: if n % i == 0: return False i += 1 return n != 1 # --- テスト --- test = [1, 2, 4, 53, 295927] ans = [[x, is_prime(x)] for x in test] print(ans) # [[1, False], [2, True], [4, False], [53, True], [295927, False]]
約数の列挙
iの初期値が1であることに注意．def divisor(n): i = 1 res = [] while i**2 <= n: if n % i == 0: res.append(i) if i != n//i: # √nより大きいものも追加 res.append(n//i) i += 1 res.sort() return res # --- テスト --- test = [1, 2, 4, 53, 295927] ans = [[x, divisor(x)] for x in test] print(ans) # [[1, [1]], [2, [1, 2]], [4, [1, 2, 4]], [53, [1, 53]], [295927, [1, 541, 547, 295927]]]
素因数分解
collections.defaultdictを使いました．8.3. collections --- コンテナデータ型 — Python 3.6.13 ドキュメント
iの初期値が2であることに注意．
①辞書的に表現すると良い点，②n //= i で2つのwhile文のnを書き換えることで簡単に書いている点は，なかなか思いつけなぁ，という感想です．
from collections import defaultdict def prime_factor(n): i = 2 # res[i]=約数iの個数 res = defaultdict(int) while i**2 <= n: while n % i == 0: # i で割れる回数を記録 res[i] += 1 n //= i i += 1 if n != 1: res[n] = 1 return res # --- テスト --- test = [1, 2, 4, 53, 295927] ans = [[x, sorted(prime_factor(x).items())] for x in test] print(ans) # [[1, []], [2, [(2, 1)]], [4, [(2, 2)]], [53, [(53, 1)]], [295927, [(541, 1), (547, 1)]]]
エラトステネスの篩（ふるい）
素数の個数
ふるいって，こんな漢字なんですね．読めない．．．この方法は，小・中学生のころ「数の悪魔」という本で読んだ記憶があります．
- 作者:ハンス・マグヌス エンツェンスベルガー
- 発売日: 2000/04/01
- メディア: 単行本（ソフトカバー）
def sieve(n): p = 0 # 素数の個数 prime = ['']*n # prime[i]=i番目の素数 # 初期化 (0,1は素数でない) is_prime = [True]*(n+1) # is_prime[0]~is_prime[n]まで必要 is_prime[0] = is_prime[1] = False for i in range(2, n+1): # 調べていない最小のiは素数 if is_prime[i]: p += 1 prime[p] = i for j in range(2*i, n+1, i): # iの倍数は素数でない is_prime[j] = False return p, prime # --- テスト --- test = [11, 1000000] ans = [[x, sieve(x)] for x in test] print(ans) # [[11, (5, ['', 2, 3, 5, 7, 11, '', '', '', '', ''])], [1000000, (78498, ['', 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, ...
区間内の素数の個数
iの倍数でa以上となる最小のもの = (a+i-1)//i * iとなることは，次の記事にまとめました．床関数・天井関数と関係式 - Notes_JP
書籍中のコードでは，$[a,b)$のふるいで，max(2, )としている部分がありました．これがないと，$j=i$となるとき（$a \leq i < \sqrt{b}$の場合に起きる）に素数$i$が除かれてしまいます．例えば$[2,5)$のときに2が素数から除かれます（i=2が一番外側のwhileで弾かれないb > 2**2の場合に，i=a=2でj=2が素数から除かれるというバグが出る）．
import math def segment_sieve(a, b): """ 入力は 2<=a < b """ prime = ['']*(b-a) # prime[i-a]=i番目の素数 # 初期化 ## [a,b)の表．(b-1)-a+1=b-a要素 ## iが素数 <=> is_prime[i-a]=True is_prime = [True for _ in range(b-a)] ## [0,√b)の表 b2 = math.ceil(math.sqrt(b)-1) # √b未満の最大の整数 is_prime_small = [True for _ in range(b2+1)] i = 2 while i**2 < b: # [0,√b)で調べていない最小のiは素数 if is_prime_small[i]: j = 2*i while j**2 < b: # iの倍数は素数でない is_prime_small[j] = False j += i # iの倍数でa以上となる最小のもの = (a+i-1)//i * i # iが素数なので，素数でないjはiの2倍以上 j = max(2, (a+i-1)//i) * i while j < b: # [0,√b)の素数の倍数を[a,b)から除く is_prime[j-a] = False j += i i += 1 prime = [i+a for i in range(b-a) if is_prime[i]] return len(prime), prime # --- テスト --- test = [[2,5], [22, 37], [22801763489, 22801787297]] ans = [[x, segment_sieve(*x)] for x in test] print(ans) # [[[2, 5], (2, [2, 3])], [[22, 37], (3, [23, 29, 31])], [[22801763489, 22801787297], (1000, [22801763489, 22801763513, 22801763527, ...
べき乗の高速計算
Carmichael Numbers
def mod_pow(x, n, mod): res = 1 while n > 0: # n=0 なら1を返す if n & 1: # i(=0,1,...)番目のビットが1なら res * x^(2^i) res = res * x % mod x = x**2 % mod # x = x^(2^(i+1)) n >>= 1 return res # --- テスト --- def is_prime(n): ''' 素数ならTrue, そうでないならFalse ''' i = 2 while i**2 <= n: if n % i == 0: return False i += 1 return n != 1 def solve(n): if is_prime(n): return 'No ({}は素数)'.format(n) for x in range(2, n): res = mod_pow(x, n, n) if res != x: return 'No ({}^{}={} (mod {}))'.format(x, n, res, n) return 'Yes' test = [17, 561, 4] ans = [[n, solve(n)] for n in test] print(ans) # [[17, 'No (17は素数)'], [561, 'Yes'], [4, 'No (2^4=0 (mod 4))']]
mod_powとして，代わりに次を使うこともできます．
def mod_pow(x, n, mod): if n == 0: return 1 res = mod_pow(x**2 % mod, n//2, mod) if n & 1: # nが奇数ならxを掛ける res = res * x % mod return res
Pythonで**とpowのどちらが速いかは，一概には言えないようです：
