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

2019 Normal

第６回

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

担当 :  aotsuki

自己紹介

自己紹介

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

KMC-ID : aotsuki

slack(内部チャット)

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

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

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

今日の内容

今日の内容

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

数学のお話

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

数学のお話

ユークリッドの互除法

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

ユークリッドの互除法

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

計算量

エラトステネスの篩

  • 自然数 N 以下の素数を高速に求める手法。
  • 具体的な計算量は O(n log ⁡log⁡n) となる。

エラトステネスの篩

  • 2≦i≦N までのループを回す。
  • 各 i について、今までに登場した素数全部で割ってみる
  • 一つでも割り切れたら素数でないので、諦める
  • どれでも割り切れなかったら素数として追加
  • 計算量はたぶん \(\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が素数だった時
}

繰り返し二乗法

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

普通の累乗

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 が、 二進数で d+1 桁で、\(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^?\) となる。
  • つまり前から i 番目の項は、
    \(x_i=0\) なら指数部が 0 になるから 1 、
    \(x_i=1\) なら指数部が \(2^i\) となるから \(a^{2^i}\) となる。

繰り返し二乗法

  • あとは何とかして高速に \(a^{2^i}\) が求められればよい
  • i−1 から 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 0 1 1 1 0 0
y 1 1 0 0 1 0 1 0
x&y 1 0 0 0 1 0 0 0
  • 両方のビットが1の時、結果が1になる

OR  (論理和)    x|y 

  • どちらかのビットが1の時、結果が1になる
x 1 0 0 1 1 1 0 0
y 1 1 0 0 1 0 1 0
x|y 1 1 0 1 1 1 1 0

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

  • 互いのビットが不一致の時、結果が1になる
x 1 0 0 1 1 1 0 0
y 1 1 0 0 1 0 1 0
x^y 0 1 0 1 0 1 1 0

右シフト    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 など) で割ったあまりを求めさせることが多い
  • 途中計算の過程でこのような数をどうやって扱っていくべきか？

剰余

  • 0,1,2,⋯,m−2,m−1 の m 個の数字しかない世界を考える
  • ここに和と積を定義する。基本的に通常の定義でよいが、その結果を m で割った余りとする
  • a+b は外の世界で言うところの (a+b)%m
  • a⋅b は外の世界で言うところの (a⋅b)%m
  • これらは、同様の計算を外の世界でしたときの最終結果を m で割った余りに等しくなる

和、積

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

引き算は？

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

割り算は？

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

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

特に、 a と p が互いに素なら

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

フェルマーの小定理

  • mが素数のとき、[0,m) 内の任意の整数 x に対しては
    　\(xx^{m-2}(=x^{m-1})\equiv1\pmod{m}\)
    が成り立つ
    • なぜなら[0,m) 内の整数は全て m と互いに素だから
  • 逆元あったじゃん！\(x\) の逆元は\(x^{m-2}\)
  • たとえばm=1000000007 のとき、これは素数なので逆元がある
  • 割り算するときは\(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から

終わり！

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

ごはん食べよう！！！

参考文献

Made with Slides.com