競技プログラミング練習会

2019 Normal

第６回

互除法, エラトステネスの篩, 繰り返し二乗法, 合同式

担当 : aotsuki

自己紹介

京都大学理学部 2回生
本名：宮嶋優大
AtCoder：水色
競プロでは主にC++を使用
ゲーム制作とweb開発をほんの少しだけ
ボードゲームとか麻雀やったりする
- 好きな人は今度一緒にやりましょう
楽しくやっていきましょう

KMC-ID : aotsuki

slack(内部チャット)

皆さんの自己紹介もお願いします

KMC-ID or 名前
所属(大学など)
プログラミング経験や競プロの経験
何か言いたいこと

ささっと、やりましょう！

今日の内容

ユークリッドの互除法
- 拡張ユークリッドの互除法
エラトステネスの篩
繰り返し二乗法
合同式

数学のお話

競プロでは数学の問題がよく出る
今日はその中でも整数に関連した内容を扱います
ちゃんと理解しなくても使うことは出来る
- しかし......
  - どういう時に使うか理解する
  - 何となく分かったレベルで良い
皆さん数学は得意ですか...？

数学のお話

ユークリッドの互除法

正の整数 , に対して最大公約数gcd(a,b) を求める。
- gcd (Greatest Common Divisor)

ユークリッドの互除法

ところでをで割ったあまりを r とすると(a>b)、 gcd(a,b)=gcd(b,r)
- 最大公約数をとして a=a'g , b=b'g とおく
- r は a÷b のあまりなので、整数で r=a−nb と表せる
- 代入すると r=a'g−nb'g=(a'−nb')g となり、で割り切れる
- とは互いに素より、とは互いに素
- よって g は最大公約数

ユークリッドの互除法

再帰関数を使うとできる
実は a>b でなくても大丈夫。試すと分かる
a % b == 0 となるときは b が最大公約数なので、それを返す
- 実際にはもう一回呼び出して返しているが、動作を追うと同じことが分かる

コード例

ll gcd(ll a, ll b) {
    if (b == 0) {
        return a;
    } else {
        return gcd(b, a % b);
    }
}

O(log max{a,b})
実際に展開すると、 2 回ごとに第一引数が半分以下になる
ラメの定理「最悪でも小さい方の十進法での桁数の 5 倍繰り返せば最大公約数に達する」

計算量

拡張
ユークリッドの
互除法

整数解 x',y' が分かっていたとする
r=a とし、ある n について r=a−nb とする
bx'+ry'=gcd(a,b)の整数解 x',y' が求まっているとする
- つまり、bx'+(a−nb)y'=gcd(a,b)
a と b について整理すると ay'+b(x'−ny')=gcd(a,b)
x=y', y=x'−ny' として構成できた。
どこかで b=0 となるので、そのときは (x,y)=(1,0) とすればよい。

ax+by=gcd(a,b)を解く

コード例

// 戻り値は gcd, xd と yd に解が入る。
ll extgcd(ll a, ll b, ll &xd, ll &yd) {
    if (b == 0) {
        xd = 1, yd = 0;
        return a;
    }
    ll gcd = extgcd(b, a % b, yd, xd);
    yd -= (a / b) * xd;
    return gcd;
}

やっていることはユークリッドの互除法と同じ

O(log max{a,b})

計算量

エラトステネスの篩

自然数以下の素数を高速に求める手法。
具体的な計算量はとなる。

エラトステネスの篩

までのループを回す。
各について、今までに登場した素数全部で割ってみる
一つでも割り切れたら素数でないので、諦める
どれでも割り切れなかったら素数として追加
計算量はたぶん \(\sum^N_{k=2}{\frac{k}{\log_e{k}}}\)くらい？

愚直にやると？

for (int i = 2; i <= N; i++) {
    bool is_prime = true;
    for (int prime : primes) {
        is_prime = is_prime && (i % prime != 0);
    }
    if (!is_prime) continue;
    primes.push_back(i);
}

エラトステネスの篩

まず「素数かどうか」を保持しておく要素数Nの配列 is_prime を持つ
- is_prime[i] = true ならiは素数
- is_prime[i] = false ならiは合成数
- 最初、中身は全て trueにしておく
変数 p を用意し、最初 p = 2 としておく
結果を入れる用の配列も用意しておく

エラトステネスの篩

is_prime[p] == true となるまで p を進める
- (最初から true のときはそのまま)
そのように進めて最初に出会った p は素数になっている
is_prime の、「p の倍数」番目の要素を全て false とする
p を結果用の配列に追加し、 p++ する

以下を繰り返す

エラトステネスの篩

コード例

vector<int> res;    //求めた素数を入れる配列
vector<bool> is_prime(N + 1, true);
for (int p = 2; p <= N; p++) {
    if (!is_prime[p]) continue;
    // p の倍数を全て消す。
    for (int k = p * 2; k <= N; k += p) is_prime[k] = false;
    res.push_back(p);
}

boolはtrueまたはfalseのみを持つ型
配列is_primeの型がboolなのはメモリを削減するため

素数判定

単に整数Zが素数かどうか判定するだけなら\(O(\sqrt Z)\) で求まる
- 2から\(\sqrt Z(+1)\)までの整数全てで割り切れなかったら素数

int z;    //zは何らかの整数
int S=sqrt(z);
int f=0;
for (int p = 2; p <= S; p++) {
    if(z%p==0)f=1;
}
if(f==0){
    //zが素数だった時
}

繰り返し二乗法

累乗の指数回かけ算するので、指数の最大値をとすると計算量。
多くの人に累乗するプログラムを書いてもらったら、こう書くのでは？

普通の累乗

ll pow_naive_for(ll a, int b) {
    ll res = 1;
    for (int i = 0; i < b; i++) { res *= a; }
    return res;
}

しかし...

繰り返し二乗法とは

\(N>10^{10}\) だと確実に死んでしまう...
どうしよう...
そんな時に繰り返し二乗法を使おう！

繰り返し二乗法とは

例えば、「\(3^{16}\)を求めよ」って言われたらどうしますか
3を16回かけますか？
多くの人は\(3\times 3=9, 9\times9=81, 81\times81=6561,\\6561\times6561=43046721\) とするのでは？
このように、繰り返し２乗することで高速に累乗計算する方法を繰り返し二乗法と言います

簡単な例で見ていきましょう

繰り返し二乗法とは

今ので、\(a^{2^x}\)についてなら高速に求められそう
じゃあ\(3^{15}\)ならどうするの？？
- \(3^{2^x}\)の積の形で表せる？
  - \(3^{15}=3^1\times3^2\times3^4\times3^8\)
  - 表せる！！
- 指数部分を2進数だと思ってやれば、任意の場合に大丈夫そう

指数が、 二進数で桁で、\(x_dx_{d-1}...x_2x_1x_0\) と表せるとする。
- \(x_0,x_1,..,x_d\) は二進数なので 全て 0 か 1
\(x=x_02^0+x_12^1+...+x_d2^d\)
指数法則より、\(a^x=a^{x_02^0}\times a^{x_12^1}\times a^{x_22^2} \times ...\times a^{x_d2^d}\)

繰り返し二乗法

\(a^x=a^{x_02^0}\times a^{x_12^1}\times a^{x_22^2}\times ... \times a^{x_d2^d} \)
\(a^x={(a^{2^0})}^{x_0}\times {(a^{2^1})}^{x_1} \times {(a^{2^2})}^{x_2} \times ... \times {(a^{2^d})}^{x_d} \)
こので区切られた部分を今だけ項と呼ぶ。
各項について、指数部は \(0\) か \(2^?\) となる。
つまり前から番目の項は、
なら指数部がになるから、
なら指数部が \(2^i\)となるから \(a^{2^i}\) となる。

繰り返し二乗法

あとは何とかして高速に \(a^{2^i}\) が求められればよい
からへの遷移でどうすればよいか？と考えると、\(a^{2(i-1)}=a^{2^i2^{-1}}=(a^{2^i})^{\frac{1}{2}}=\sqrt{a^{2^i}}\)から \(a^{2^i}\)なので、二乗すればよい
二乗くらいなら \(r\times r\) として一瞬
ちなみに初項は、\(a^{2^0}=a^1=a\) でよい

繰り返し二乗法

コード例

ll pow_ans(ll x, int n) {
    ll res = 1;
    while (n > 0) {
        if (n & 1) { res *= x; }
        x *= x;
        n >>= 1;
    }
    return res;
}

累乗はすぐオーバーフローしてしまうので、 long long など大きい型を使おう
大きい数 (1000000007 など) で割ったあまりを求めさせる問題も多い
- 剰余は 最後だけではだめ！途中計算の段階で既にオーバーフローしちゃってるので意味のない数字になる
- 途中計算の過程でもこまめに剰余をとって、どの計算でもオーバーフローしないように

繰り返し二乗法

ビット演算

（復習）

ビット演算とは

各のビット全てに対する論理演算をいっぺんに行う演算操作

いろんなところで使える

Binary Indexed Tree (BIT)
bit 全探索
bit DP
...

ビット演算の利点

データ量が少なくて済む
多少速い
一つのデータに複数の情報を詰め込める(フラグ管理など)

ビット演算子の種類

NOT (否定) ~x
AND (論理積) x & y
OR (論理和) x | y
XOR (排他的論理和) x ^ y
右シフト x >> n
左シフト x << n

NOT (否定) ~x

x	0	1	1	0	0	1	0	1
~x	1	0	0	1	1	0	1	0

1の補数
ビットを反転させる（0→1/1 →0）

AND (論理積) x&y

x	1	0	1	1	1	0
y	1	1	0	1	0	1
x&y	1	0	0	1	0	0

両方のビットが1の時、結果が1になる

OR (論理和) x|y

どちらかのビットが1の時、結果が1になる

x	1	0	1	1	1	0
y	1	1	0	1	0	1
x\|y	1	1	1	1	1	1

XOR (排他的論理和) x ^ y

互いのビットが不一致の時、結果が1になる

x	1	0	1	1	1	0
y	1	1	0	1	0	1
x^y	0	1	1	0	1	1

右シフト x >> n

各ビットを右にn個ずらす
　　　と同じ

x	1	0	0	1	1	1	0	0

x>>n	0	0	1	0	0	1	1	1

{

n=2

x : unsigned char

空いた桁は0で埋まる

このbitは捨てられる

x/{2^n}

(ほぼ)

左シフト x << n

x	1	0	0	1	1	1	0	0

x<<n	0	1	1	1	0	0	0	0

{

n=2

x : unsigned char

空いた桁は0で埋まる

このbitは捨てられる

各ビットを左にn個ずらす
桁があふれなければ　　　と同じ

{x}*{2^n}

合同式

競プロで剰余を目にすることはかなり多い！
フィボナッチ数や組み合わせの数などを求めさせる系の問題では答えが爆発的に大きくなる
大きい素数 (1000000007 など) で割ったあまりを求めさせることが多い
途中計算の過程でこのような数をどうやって扱っていくべきか？

剰余

の個の数字しかない世界を考える
ここに和と積を定義する。基本的に通常の定義でよいが、その結果をで割った余りとする
は外の世界で言うところの
は外の世界で言うところの
これらは、同様の計算を外の世界でしたときの最終結果をで割った余りに等しくなる

和、積

実は引き算も同様にできる
- ただし負の数の余りを求めてしまわないように注意
例 : の世界で -> ...？
繰り下がりのある引き算と一緒で、先にだけ借りてくる
つまり (7+(3−5))%7=5

引き算は？

「で割る」というのは、「乗法の逆元 \(x^{-1}\) をかける」と一般化できる。
の逆元は、\(xx^{-1}=1\) となるような \(x^{-1}\) のこと。
例えば実数の世界では、の逆元は \(\frac{1}{x}\)
今考えている世界ではの逆元は \(xx^{-1}\equiv0\pmod{m}\\\)となるような \(x^{-1}\) のこと。
しかし、この世界に逆元なんてあるのか？どうやって求めるのか？

割り算は？

を任意の自然数、を素数とする。このとき、

\(a^p\equiv a\pmod{p}\)

特に、とが互いに素なら

\(a^{p-1}\equiv 1\pmod{p}\)

フェルマーの小定理

mが素数のとき、内の任意の整数に対しては
　\(xx^{m-2}(=x^{m-1})\equiv1\pmod{m}\)
が成り立つ
- なぜなら内の整数は全てと互いに素だから
逆元あったじゃん！\(x\) の逆元は\(x^{m-2}\)
たとえばのとき、これは素数なので逆元がある
割り算するときは\(x^{1000000005}\) を計算して掛ければよい
これには先の繰り返し二乗法が使える

つまり

二項係数

二項係数とは
- 組み合わせの数
- \({}_n\mathrm{C}_k\)
- \(=\frac{n!}{k!(n-k)!}\)
- \(=n!(k!)^{-1}{(n-k)!}^{-1}\)
累積積とフェルマーの小定理を使えばいける！

二項係数

コンテスト

コンテストはAizu Online Judgeで行います
まだしてない人は会員登録をしましょう
下記のURLにアクセスしてください

コンテスト

http://judge.u-aizu.ac.jp/onlinejudge/index.jsp?lang=ja

もしくは

「 AOJ 」で検索

下記のURLにアクセスしてください
パスワードは　"y"　です
コンテスト開始までしばしお待ちください

コンテスト

https://onlinejudge.u-aizu.ac.jp/beta/room.html#kmc2019_n_6

もしくは

部内wikiから

終わり！

解説して欲しいものありますか？

ごはん食べよう！！！

競技プログラミング練習会2018 Normal 第７回スライド

参考文献

第６回：互除法, エラトステネスの篩, 繰り返し二乗法, 合同式

By procon2019

第６回：互除法, エラトステネスの篩, 繰り返し二乗法, 合同式

発表日時　2019年5月25日(金) 18:30-21:00 https://onlinejudge.u-aizu.ac.jp/beta/room.html#kmc2019_n_6

1,000

procon2019

このアカウントは京大マイコンクラブ(KMC)の2019年度の競技プログラミング練習会Normal(初心者向け)での説明用に作成したスライドです。閲覧・参照はご自由にどうぞ。

競技プログラミング練習会

2019 Normal

自己紹介

自己紹介

皆さんの自己紹介もお願いします

今日の内容

今日の内容

数学のお話

数学のお話

ユークリッドの互除法

ユークリッドの互除法

ユークリッドの互除法

コード例

計算量

拡張 ユークリッドの 互除法

ax+by=gcd(a,b)を解く

コード例

計算量

エラトステネスの篩

エラトステネスの篩

愚直にやると？

エラトステネスの篩

エラトステネスの篩

エラトステネスの篩

コード例

素数判定

繰り返し二乗法

普通の累乗

繰り返し二乗法とは

繰り返し二乗法とは

繰り返し二乗法とは

繰り返し二乗法

繰り返し二乗法

繰り返し二乗法

コード例

繰り返し二乗法

ビット演算

合同式

剰余

和、積

引き算は？

割り算は？

フェルマーの小定理

つまり

二項係数

二項係数

コンテスト

コンテスト

コンテスト

終わり！

解説して欲しいものありますか？

ごはん食べよう！！！

参考文献

第６回：互除法, エラトステネスの篩, 繰り返し二乗法, 合同式

More from procon2019

拡張
ユークリッドの
互除法