Math Skill

定義與假設

如果沒有特別說明，討論的變數定義範圍是 0 或正整數 $\N_0$

定義與假設

如果有出現 pow ，指的是快速冪

ll qpow(ll a, ll e)
{
    ll mod = e;
    ll res = 1;
    
    while( e )
    {
        if( e&1 ) res = res * a % mod;
        e/=2;
        a=a*a%mod;
    }
    return res;
}

Prime

質數

因數

如果 $x=ab$

則 $a,b$ 是 $x$ 的因數

除了$1$，每一個整數都有兩個顯然因數

$1,x$

因數

整除符號

如果 $a$ 是 $x$ 的因數，則使用下方法紀錄

$a\mid x$

example

$15=5\times 3$

$\therefore 3|15$

$\therefore 5|15$

example

不整除符號

$15\nmid 3$

$2\nmid 7$

Hint : 數字大的放後面，所以前面是因數

質數

如果一個數字 $n$ 是質數，則 $n$ 只有兩個相異的因數

$7$ 是質數，因數有 $1,7$

$6$ 不是質數，因數有 $1,2,3,6$

$1$ 不是質數，因為沒有相異的因數

$0$ 不在討論範圍

質數判定

判斷 $n$ 是不是質數

質數判定

判斷 $n$ 是不是質數

暴力方法，用 $2\sim n-1$ 的數字試試看

bool isp(int n)
{
    for(int i=2;i<n;++i)
    {
        if ( n%i==0 )
            return false;
    }
    return true;
}

時間複雜度：$O(n)$

質數判定

判斷 $n$ 是不是質數

定理 1. 若 $n=ab$, 且 $a\leq b$ 則

$a\leq \sqrt{n}$

如果 $n$ 是合數，會有一個因數小於等於 $\sqrt{n}$

可用反證法輕易得證

質數判定

判斷 $n$ 是不是質數

bool isp(int n)
{
    int s = (int) ceil(sqrt(n));
  
    for( int i=2 ; i<=s ; ++i )
    {
        if ( n%i==0 )
            return false;
    }
    return true;
}

// Way 2  
    for( int i=2 ; i*i<=n ; ++i )

注意判斷條件 <=

注意溢位 (way2)

不要把 sqrt 放到 for 上

時間複雜度：$O(\sqrt{n})$

sqrt 可能的時間複雜度為 $O(\log {n})$

質數判定

判斷 $x$ 是不是質數

預處理法：透過預先計算的資訊加速答案查詢

當中最有名的方法為 κόσκινον Ἐρατοσθένους

又稱為埃氏篩法 (sieve of Eratosthenes)

可以預先計算 $1\sim n$ 之間的所有數字是否為質數

質數判定

判斷 $x$ 是不是質數

埃氏篩法

先假設所有數字都是質數，0, 1 除外

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

質數判定

判斷 $x$ 是不是質數

埃氏篩法

對於一個數字 $x$，它的倍數 $2x, 3x \dots$ 一定不是質數

從 $x=2$ 開始檢查數字

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i<=x;++i)
{
    for (int j=i*2;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i<=x;++i)
{
    for (int j=i*2;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

定理. 通過此方法後，每一個合數都會被標記為 false

prof. 因為每一個合數都會被它的因數給標記

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i<=x;++i)
{
    for (int j=i*2;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

時間複雜度

$n/2+n/3+n/4...n/n$

$< n(1/2+1/3+1/4+...1/n)$

$\in O(n\ln n)$

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i*i<=x;++i)
{
    for (int j=i*2;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

加速 1. 外迴圈只需要到 $\sqrt{x}$

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i*i<=x;++i)
{
    for (int j=i*i;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

加速 1. 外迴圈只需要到 $\sqrt{x}$

加速 2. 內迴圈只需要從 $i\times i$

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i*i<=x;++i)
{
    if (!isp[i])
        continue;
    for (int j=i*i;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

加速 1. 外迴圈只需要到 $\sqrt{x}$

加速 2. 內迴圈只需要從 $i\times i$

加速 3. 不是質數不用篩

質數判定

判斷 $x$ 是不是質數

埃氏篩法

vector<int> isp(n+1, true);

isp[0] = isp[1] = false;

for (int i=2;i*i<=x;++i)
{
    if (!isp[i])
        continue;
    for (int j=i*i;j<=x;j+=i)
    {
        isp[j] = false;
    }
}

加速 1. 外迴圈只需要到 $\sqrt{x}$

加速 2. 內迴圈只需要從 $i\times i$

加速 3. 不是質數不用篩

時間複雜度

$n/2+n/3+n/5...n/p$

$< n(1/2+1/3+1/5+...1/p)$

$\in O(n\ln\ln n)$

@質數倒數和

質數判定

判斷 $n$ 是不是質數

暴力判定 $O(\sqrt{n})$

質數篩法 $O(n\ln\ln n)$

歐拉篩法 $O(n)$

質數判定

判斷 $x$ 是不是質數

唬爛方法

利用一些質數才成立的公式來檢查 $x$ 是不是質數

example

是不是奇數

是不是 $6n\pm 1$

質數判定

判斷 $x$ 是不是質數

唬爛方法

警告：質數會成立，不代表非質數就不成立！

利用一些質數才成立的公式來檢查 $x$ 是不是質數

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

example

隨意取一個 $a$ ，則 $a^5\equiv a \mod{5}$

$2^5=32\equiv 2 \mod{5}$

$3^5=243\equiv 3 \mod{5}$

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

Prove.

考慮 $(b+1)^p \mod{p}$

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

Prove.

$$(b+1)^p =$$

$$\binom{p}{p}b^p+\binom{p}{p-1}b^{p-1}+\binom{p}{p-2}b^{p-2}+\binom{p}{p-3}b^{p-3}+\cdots+\binom{p}{2}b^{2}+\binom{p}{1}b^{1}+\binom{p}{0}b^{0}$$

因為 $\binom{p}{n}=\frac{p!}{n!(p-n)!}$，否則除非 $n=0,p$，分子的$ p$ 都無法被消除

所以 $\binom{p}{n}$ 會是 $p$ 的倍數

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

Prove.

$$(b+1)^p =$$

$$\binom{p}{p}b^p+\binom{p}{p-1}b^{p-1}+\binom{p}{p-2}b^{p-2}+\binom{p}{p-3}b^{p-3}+\cdots+\binom{p}{2}b^{2}+\binom{p}{1}b^{1}+\binom{p}{0}b^{0}$$

$$\equiv \binom{p}{p}b^p+\binom{p}{0}b^{0} = b^p+1\pmod{p}$$

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

Prove.

$$\equiv (b-1)^p+1+1\pmod{p}$$

$$\equiv (b-2)^p+1+1+1\pmod{p}$$

$$(b+1)^p \equiv b^p+1\pmod{p}$$

$$\equiv (b-3)^p+1+1+1+1\pmod{p}$$

$$\equiv\underbrace{1+1+1+1+1+1+1+1}_{b+1}\pmod{p}$$

$$\equiv b+1\pmod{p}$$

費馬小定理

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

Prove.

$$(b+1)^p \equiv b+1\pmod{p}$$

設 $a=b+1$ 得

$$a^p\equiv a \pmod{p} ~_\square$$

其實比較常用

$$a^{p-1}\equiv 1 \pmod{p}$$

費馬測試

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

隨機生成一些 $a$ 檢查有沒有滿足費馬小定理

Try ZJ a007

費馬測試

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

隨機生成一些 $a$ 檢查有沒有滿足費馬小定理

mt19937 mt(7122);
bool isp(ll x)
{

    ll a = mt() % (x-1) + 1;
    if( qpow( a, x ) != a ) return false;

    return true;
}

費馬測試

如果 $p$ 是質數，則對於正整數 $1\leq a\leq p-1$ 有

$$a^p\equiv a \pmod{p}$$

隨機生成一些 $a$ 檢查有沒有滿足費馬小定理

mt19937 mt(7122);
bool isp(ll x)
{

    ll a = mt() % (x-1) + 1;
    if( qpow( a, x ) != a ) return false;

    return true;
}

卡邁克爾數

一些不是質數，但是永遠滿足費馬測試的數字

$561,1105,1729,2465,2821,6601\cdots$ @

可以枚舉卡邁克爾數然後在判質數嗎？

UVA 10006

雖然這種數字很少，$1\sim n$ 之間約有 $n^{2/7}$ 個

但沒有簡單的列舉方法，long long 範圍內有 60 萬個

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1 \pmod p$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, 只有 $x=\pm 1$時成立 (方根測試)

對於第二點，帶入 $\pm 1$ 顯然成立

因此多挑一些 $x$ 同時檢驗兩件事情就好！

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1 \pmod p$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, 只有 $x=\pm 1$時成立 (方根測試)

prof.

$x^2\equiv 1$

$x^2-1 \equiv 0$

$(x-1)(x+1)\equiv 0$

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1 \pmod p$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, 只有 $x=\pm 1$時成立 (方根測試)

prof.

$x^2\equiv 1$

$x^2-1 \equiv 0$

$(x-1)(x+1)\equiv 0$

顯然的有平凡解 +1, -1

因為 $p$ 是質數

故不存在一個正整數 $1< t<p$

使得 $ p\equiv 0 \pmod t$

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1 \pmod p$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, 只有 $x=\pm 1$時成立 (方根測試)

prof.

$x^2\equiv 1$

$x^2-1 \equiv 0$

$(x-1)(x+1)\equiv 0$

顯然的有平凡解 +1, -1

所以，如果 (x-1), (x+1) 不等於 0

因為都沒有 $p$ 的因數

乘起來也不會是 $p$ 的倍數 $_\square$

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1 \pmod p$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, 只有 $x=\pm 1$時成立 (方根測試)

對於第二點，帶入 $\pm 1$ 顯然成立

因此多挑一些 $x$ 同時檢驗兩件事情就好！

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

優化：可以在算費馬測試時，順便測試第二點！

因為 $p-1$ 是偶數，可以將 $p-1$ 分解為

$d\times 2^u$

先算 $a^d$，再平方 $u$ 次，沿途產生的數字都來測試！

Note : $u\geq 1$

以下是加速技巧

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

因為 $p-1$ 是偶數，可以將 $p-1$ 分解為

$d\times 2^u$

計算 $a^d$

如果 $a^d \equiv \pm 1$ ，則一定通過費馬測試

則 $a^{d2^e}$一定通過方根測試

是質數

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

因為 $p-1$ 是偶數，可以將 $p-1$ 分解為

$d\times 2^u$

計算 $a^{d2^e}$

$a^{d}\neq\pm 1, e<u$

如果 $a^{d2^e}=1$ 則不通過方根測試

因為 $x=2^{d2^{e-1}}\neq\pm1$ 卻滿足方程式

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

因為 $p-1$ 是偶數，可以將 $p-1$ 分解為

$d\times 2^u$

計算 $a^{d2^e}$

如果 $a^{d2^e}=-1$ 則一定通過費馬測試

$a^{d2^e}$ 則一定通過方根測試

$a^{d}\neq\pm 1, e<u$

因為 $x=2^{d2^{e-1}} = -1$ ，是方程式預設的解

是質數

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試 $a^{p-1} \equiv 1$

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

因為 $p-1$ 是偶數，可以將 $p-1$ 分解為

$d\times 2^u$

計算 $a^{d2^u}$

如果以算到 $p-1$ 次方還沒判定是質數

就不是質數了

因為我們總是預測先預測下一步是否成立

米勒-拉賓測試

bool miller_rabin(int n, int a)
{
    if ( n==2 ) return true;
    if ( n < 2 || n % 2 == 0 ) return false;
    // 分解次方
    int u = n-1, t = 0;
    while (u % 2 == 0) u >>= 1, t++;
 
    // 先檢查 a^u
    int x = pow(a, u, n);   // x = a ^ u % n;
    if (x == 1 || x == n-1) return true;
 
    for (int i=0; i<t-1; i++) // 算密次剩下的部分，最後一次不用算
    {
        x = mul(x, x, n);   // x = x * x % n;
        if (x == 1) return false;  //方根測試失敗
        if (x == n-1) return true; //偷窺未來一定成功
    }
    return false; // for會偷窺歸下一步會不會成功，所以最後還沒成功就是失敗
}

米勒-拉賓測試

定理：如果 $p$ 不是質數，任選一個 $x$ 同時滿足米勒-拉賓測試的機率 $\leq 1/4$

Michael O Rabin 1980. Probabilistic algorithm for testing primality

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

定理：如果 $p$ 不是質數，任選一個 $x$ 同時滿足兩件事情的機率 $\leq 1/4$

因此，單次的正確率是 3/4

事實上，有人先找出了夠好了 $x$ 來檢驗

米勒-拉賓測試

如果 $p$ 是質數，且 $p>2$，則

1. 通過費馬測試

2. 對於方程式 $x^2\equiv 1 \pmod{p}$, $x=\pm 1$ 是唯一的解

int 範圍：選 $x=\{2,7,61\}$

long long 範圍：選 $x=$ 小於 $40$ 的所有質數

wiki = $\displaystyle 2,325,9375,28178,450775,9780504,1795265022$

Factor

因數分解

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

ex. 12 的因數有？

1,2,3,4,6,12

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

暴力法：$O(\sqrt{x})$

ll factor(ll x)
{
    for(ll i=2;i*i<=x;++i)
    {
        if( x%i==0 ) return i;
    }
    return x; //x is prime
}

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

建表法：$O( n\log\log n + \log{x})$

在做篩法的時候，幫每一個數字紀錄被誰除過

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

建表法：$O( n\log\log n + \log{x})$

iota(dp, dp+n, 0);// 初始化 dp[i] = i; 
for(ll i=2;i*i<=n;++i)
{
    if( isp[i] )
    {
        for( int j=i*i; j<=n; j+=i )
        {
            isp[j] = false;
            dp[j] = i; // j 的因數有 i 
        }
    }
}

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

混和大法：$O( \sqrt{n}\log\log \sqrt{n} + \pi(\sqrt{x}))$

如果數字大到建表建不下，可以至少建到$\sqrt{n}$

再用暴力法，但是只嘗試質數

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

唬爛法

random 一個數字來猜猜看他是不是因數

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

唬爛法

mt19937 mt(7122);
ll factor(ll x)
{
    int f;
    uniform_int_distribution<int> uni(2, x-1);
    do{
        f = uni(mt);
    }while( x%f!=0 );
    
    return f;
}

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

唬爛法

最慘的狀況下， $x$ 只有兩個因數，故命中率為

$$\frac{2}{x-2}$$

note. 1跟x 不用嘗試

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

唬爛法

如果 $x=pq$

$p,q$ 一定有一個數字小於等於 $\sqrt{x}$

所以可以縮小 random 的範圍提高機率

因數問題

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

唬爛法

mt19937 mt(7122);
ll factor(ll x)
{
    int f;
    uniform_int_distribution<int> uni(2, (int)sqrt(x+1));
    do{
        f = uni(mt);
    }while( x%f!=0 );
    
    return f;
}

命中率提高為

$$\frac{1}{\sqrt{x}}$$

生日攻擊

班上有 $23$ 位同學，當中有兩位同學同天生日的機率是多少？

$n$ 人選 $n^2$ 個東西，相同率為 $50\%$

$n$ 人選 $1\sim n^2$ 內的數字，相差 $\sqrt{n}$ 的機率約為 $50\%$

$n$ 人選 $1\sim n^2$ 內的數字，相同的機率約為 $50\%$

Pollard’s ρ Algorithm

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

透過亂數產生一個數列 $a_1,a_2,a_3,\cdots a_k$ $(1<a_i<x)$

判斷 $a_i-a_j$ 與 $x$ 的 $\gcd$ 是不是 $1$

只找 p 太浪費，連同 2p, 3p, 4p ... 一起找提高機率

至少 $\sqrt{n}$ 個選擇

如果不是就猜到了因數了

Pollard’s ρ Algorithm

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

有人說 $a_i = a_{i-1}^2+c \mod x$ 是一個很好的序列函數

玄學，不知道為什麼

因為這個公式會循環，要判斷循環節，如果已經循環還沒找到就失敗了 (換一個c繼續)

Pollard’s ρ Algorithm

給定一個數字 $x$ ，請回答一個 $x$ 的一個不等於 $1$ 或 $x$ 的因數

可以假設 $x > 2$ 且 $x$ 不是質數

Floyd 判環法

又稱龜兔賽跑法

$x=f(x), y=f(f(y))$

讓一個兔子，一次走兩步

讓一個烏龜，一次走一步

如果出發後兔子遇到烏龜就找到循環

Brent 判環法

比前一個快一些

$x=f(x), y=x \text{ if at $2^k$ state} $

讓一個烏龜，一次走一步

讓一個兔子，在烏龜走第$2^k$ 步時，直接飛到烏龜上

其他時間睡覺不動

Pollard’s ρ Algorithm

int c, mod;
ll f(ll x)
{
    return (x*x+c)%mod;
}

ll factor(ll x)
{
    c = mt()%x;
    mod = x;
    
    ll a=2, b=2;
    
    while(true)
    {
        a = f(a);
        b = f(f(b));
        
        if( a==b ) return -1; // fail
        if( __gcd(x, abs(a-b))!=1 )
            return abs(a-b);
    }
}

Pollard’s ρ Algorithm

期望複雜度：$O(\sqrt{\sqrt{x}})=O(x^{1/4})$

因為在 $1\sim n$ 的範圍內猜 $\sqrt{n}$ 種可能

相當於在 $1\sim \sqrt{n}$ 的範圍內猜 $1$ 種可能

根據生日攻擊，只需要猜$\sqrt{\sqrt{x}}$ 次就有 $50\%$ 命中率

很唬爛的證明，但完整證明太難，從略

Pollard’s ρ Algorithm

期望複雜度：$O(\sqrt{\sqrt{x}})=O(x^{1/4})$

如果失敗就多換幾個 $c$

如果 $x$ 是質數，會進入無窮迴圈

- 用米勒-拉賓預處理

Euler's totient function

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

定義 $a,b$ 互質 $\iff \gcd(a,b)=1$

$\varphi(1)=1$

$\varphi(5)=4$

$\varphi(6)=2$

$\varphi(30)=?$

$\varphi(7)=?$

$\varphi(2147483647)=?$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

性質1. 如果 $p$ 是質數，則

$\varphi(p)=p-1$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

性質2. 如果 $p$ 是質數，則

$$\varphi(p^k)=p^k-p^{k-1}=p^k(1-\frac{1}{p})$$

因為 $p$ 的倍數有 $p^k/p=p^{k-1}$ 這麼多個

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

性質2. 如果 $a,b$ 互質，則

$\varphi(ab)=\varphi(a)\varphi(b)$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

性質3. 如果 $a,b$ 互質，則

$\varphi(ab)=\varphi(a)\varphi(b)$

Prove. 若 $p,q$ 是質數，則

$\varphi(pq)=pq-1-(q-1)-(p-1)$

$=pq-p-q+1$

$=(p-1)(q-1)$

$=\varphi(p)\varphi(q)$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

性質3. 如果 $a,b$ 互質，則

$\varphi(ab)=\varphi(a)\varphi(b)$

Prove. 若 $p,q$ 是質數，則同理可得

$\varphi(p^xq^y)=\varphi(p^x)\varphi(q^y)$

將 $a,b$ 質因數分解後，將全部合併即可得證

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

如何計算 $\varphi(x)$

找自己的因數很麻煩

但是找自己的倍數很簡單

把自己的倍數都扣掉自己的 $\varphi(x)$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

如何計算 $\varphi(x)$

void phiTable(){
	for(int i=1;i<=N;i++)
            phi[i]=i;

	for(int i=1;i<=N;i++)    //對於每一個數由小到大
            for(x=i*2;x<=N;x+=i) // 每一個倍數
                phi[x]-=phi[i];  // 都扣掉自己的 phi
}

$O(n\log\log n) 打表$

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

如何計算 $\varphi(x)$

對 $x$ 做質因數分解後，用性質 2,3 暴力算

性質3. 如果 $a,b$ 互質，則

$\varphi(ab)=\varphi(a)\varphi(b)$

性質2. 如果 $p$ 是質數，則

$$\varphi(p^k)=p^k-p^{k-1}=p^k(1-\frac{1}{p})$$

快速冪！

$O(\sqrt{n}) $ 算一個

歐拉函數

歐拉函數 $\varphi(x)$ 的定義為

小於等於 $x$ 的正整數中，有多少數字與 $x$ 互質

小練習

https://toj.tfcis.org/oj/pro/427/

應該10~20行可ac !?

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

當 $n$ 是質數時就是費馬小定理

$$a^{p-1}\equiv 1 \pmod{p}$$

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

Prove：假設1~n之間與 n 互質的數字是

$\{x_1, x_2,x_3,\dots, x_{\varphi(n)}\}$

1. 這集合有所有與 n 互質的數字

2. 這集合任兩數字都不一樣

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

同乘 $a$

$\{ax_1, ax_2,ax_3,\dots, ax_{\varphi(n)}\}$

Prove：假設1~n之間與 n 互質的數字是

$\{x_1, x_2,x_3,\dots, x_{\varphi(n)}\}$

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

1. 這集合有所有與 n 互質的數字

2. 這集合任兩數字都不一樣

Prove：假設1~n之間與 n 互質的數字是

$\{x_1, x_2,x_3,\dots, x_{\varphi(n)}\}$

$\{ax_1, ax_2,ax_3,\dots, ax_{\varphi(n)}\}$

所以把新的序列 $\mod n$ 會是舊序列重新排列組合的一種

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

Prove：假設1~n之間與 n 互質的數字是

$\{x_1, x_2,x_3,\dots, x_{\varphi(n)}\}$

$\{ax_1, ax_2,ax_3,\dots, ax_{\varphi(n)}\}$

$$\prod_{\varphi(n)}{x_i}\equiv\prod_{\varphi(n)}{ax_i}\equiv a^{\varphi(n)}\prod_{\varphi(n)}{x_i} \mod n$$

歐拉定理

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

Prove：同除 $\prod_{\varphi(n)}{x_i}$

$$\prod_{\varphi(n)}{x_i}\equiv\prod_{\varphi(n)}{ax_i}\equiv a^{\varphi(n)}\prod_{\varphi(n)}{x_i} \mod n$$

$$1\equiv a^{\varphi(n)} \mod n ~$$

Hint : 因為與 $n$ 互質

歐拉降冪法

如果 $a,n$ 互質，則

$$a^{\varphi(n)}\equiv 1 \pmod{n}$$

如果 $a,n$ 互質，冪次結果會每 $\varphi(n)$ 循環一輪

如果 $a,n$ 不互質呢？

$$a^k=a^{k\%\varphi(n)} \pmod{n}$$

其實還是會每 $\varphi(n)$ 循環一輪，但第一輪可能不在循環節內

歐拉降冪法

對於任意的 $a,n$，降冪公式如下

$$a^k\equiv{\begin{cases}a^k & k<\varphi(n) \\a^{k\%\varphi(n)+\varphi(n)} & k\geq \varphi(n)\end{cases}}$$

證明太長，自行 Google

歐拉降冪法

對於任意的 $a,n$，降冪公式如下

有了降冪法，可以保證指數的部分最多只會是 $2n$

可以避免一些題目的大數 (雖然應該是要考降冪法)

$$a^k\equiv{\begin{cases}a^k & k<\varphi(n) \\a^{k\%\varphi(n)+\varphi(n)} & k\geq \varphi(n)\end{cases}}$$

歐拉降冪法

高冪次的求解

$$a^{b^{c^d}}\equiv a^{b^{c^d}\%\varphi(n)+\varphi(n)} \mod n$$

遞迴求

$$b^{c^d} \mod \varphi(n)$$

要小心 $b^{c^d} < \varphi(n)$ 時不能加 $\varphi(n)$ 有機會壞掉

歐拉降冪法

https://toj.tfcis.org/oj/pro/179/

UVA 10692

要小心 $b^{c^d} < \varphi(n)$ 時不能加 $\varphi(n)$ 有機會壞掉

很多測資沒測這點 (或正解根本寫錯)，自己小心

Multiplication

大數乘法

大數運算

如果一個數字 $n$ 的大小超過了硬體限制的範圍，就無法直接當作一個單位運算
此時只能透過軟體模擬，讓以往可以忽視複雜度的部份不再可以忽略

long long a, b;
cin >> a >> b;
cout << a * b << '\n';

簡易大數

using BigNum = vector<int>;
BigNum& trim(BigNum &b) {
    while(b.size() > 1 && b.back() == 0)
        b.pop_back();
    return b;
}

BigNum stob(const string &s) {
    BigNum b;
    for(char c:s | views::reverse) // C++20
    	b.emplace_back(c - '0');
    return trim(b);
}

void show(const BigNum &b) {
    for(auto i : b | views::reverse) cout << i;
    cout << '\n';
}

只考慮正整數

views::reverse

反向遍歷

#include<ranges>

vector<int> a {1,2,3,4,5};

//old C++ : 1 2 3 4 5
for(vector<int>::iterator it=a.begin(), e=a.end();it!=e;++it)
	cout << *it << ' ';
    
//old C++ : 5 4 3 2 1 
for(vector<int>::reverse_iterator it=a.rbegin(), e=a.rend();it!=e;++it)
	cout << *it << ' ';

// C++11 1 2 3 4 5
for(int i:a) cout << i << ' ';

// C++20 5 4 3 2 1
for(int i:a | views::reverse) cout << i << ' ';

樸素大數乘法

模擬直式乘法

   1 2 
x  3 4
-------
   4 8
 3 6
-------    (結果相加)
 3 10 8
-------    (處理進位)
 4 0 8

樸素大數乘法

BigNum& advance(BigNum& b) {
    for(int i=0;i+1<b.size();++i){
        b[i+1] += b[i]/10;
        b[i]%=10;
    }
    return b;
}

BigNum mul(const BigNum &a, const BigNum &b)
{
    BigNum c(a.size() + b.size());
    for(int i=0;i<a.size();++i)
        for(int j=0;j<b.size();++j)
            c[i+j] += a[i] * b[j];
    
    return trim(advance(c));
}

很明顯 $O(n^2)$

Karatsuba algorithm

Karatsuba 是第一個比樸素乘法還要快的演算法

把數字分解成數段，分別乘起來後再組合

$x = x_1B^m+x_2$

$y = y_1B^m+y_2$

Karatsuba algorithm

$x = x_1B^m+x_2$

$y = y_1B^m+y_2$

$xy=x_1y_1B^{2m}+(x_1y_2+x_2y_1)B^m+x_2y_2$

由於 $x_1,x_2,y_1,y_2$ 是原來 $x,y$ 長度的一半

經過變換後變成做 4 個比較小的乘法

Karatsuba algorithm

$xy=x_1y_1B^{2m}+(x_1y_2+x_2y_1)B^m+x_2y_2$

由於 $x_1,x_2,y_1,y_2$ 是原來 $x,y$ 長度的一半

經過變換後變成做 4 個比較小的乘法

複雜度 $T(n)=4T(\frac{n}{2})+O(n)$

根據主定理 $O(T(n))\in O(n^{\log_2 4})=O(n^2)$

一點也沒變好 !?

Karatsuba algorithm

$xy=x_1y_1B^{2m}+(x_1y_2+x_2y_1)B^m+x_2y_2$

$x_1y_2+x_2y_1 =$

$(x_1+x_2)(y_1+y_2)-x_1y_1-x_2y_2$

其中 $x_1y_1,x_2y_2$ 已經算過了 !

只需要再算 $(x_1+x_2)(y_1+y_2)$

Karatsuba algorithm

$xy=x_1y_1B^{2m}+(x_1y_2+x_2y_1)B^m+x_2y_2$

$x_1y_2+x_2y_1 =$

$(x_1+x_2)(y_1+y_2)-x_1y_1-x_2y_2$

把一個乘法變成很多加減法

複雜度 $T(n)=3T(\frac{n}{2})+O(n)$

根據主定理 $O(T(n))\in O(n^{\log_2 3})=O(n^{1.5826})$

賺!

Karatsuba algorithm

$xy=x_1y_1B^{2m}+(x_1y_2+x_2y_1)B^m+x_2y_2$

$x_1y_2+x_2y_1 =$

$(x_1+x_2)(y_1+y_2)-x_1y_1-x_2y_2$

把一個乘法變成很多加減法

複雜度 $T(n)=3T(\frac{n}{2})+O(n)$

根據主定理 $O(T(n))\in O(n^{\log_2 3})=O(n^{1.5826})$

賺!

Karatsuba algorithm

除了本次課程的拆解方法，後續的研究有提出更快更噁心的拆法，多用於科學計算上

Example :

Toom-Cook Algorithm$O(n^{1.46})$

這東西用 C/C++ 實做有點長度，參見卦長模板

不過我們還要更快!

Cooley–Tukey FFT algorithm

$O(n\log n\log\log n)$

回到回到樸素大數乘法

$C=AB$

處理進位問題之前， C 的第 x 位數的值為

$C_x=A_0B_x+A_1B_{x-1}+A_2B_{x-2}\cdots+A_xB_0$

c[i+j] += a[i] * b[j];

把數字轉換成一個多項式，那大數相乘的結果與多項式相乘是類似的!

如果 $A = 789$，設 $f_A(x)=7x^2+8x+9$

如果 $B = 345$，設 $f_B(x)=3x^2+4x+5$

那 $f_A(x)f_B(x)=f_C(x)=c_0+c_1x+c_2x^2+\cdots$

其中 $c_i=a_0b_i+a_1b_{i-1}+a_2b_{i-2}+\cdots+a_ib_0$

與前頁 $C_x=A_0B_x+A_1B_{x-1}+A_2B_{x-2}\cdots+A_xB_0$ 一樣!

捲積

定義大小為 $n$ 的離散卷積($*$)如下

$f(x), g(x)$ 為函數

則 $f*g=\sum\limits^{n}_{i=0} f(i)g(n-i)$

也就是做大數乘法的動作，等於對每一個系數求卷積

單純名詞介紹，但記住這個定義可以把這演算法推廣到不僅是乘法上，只要長得像卷積就可以

即便把問題變成求多項式相乘的係數好像也沒什麼用
但是如果把多項式變成更有趣的形式有沒有幫助呢 ?

插值多項式唯一性定理

給定 $n+1$ 個點 $(x_1,y_1), (x_2, y_2),\cdots (x_{n+1}, y_{n+1})$

當中 $x$ 皆不相等

則存在唯一的 $n$ 次多項式經過這些點

by 高中數學 : 牛頓法 / 拉格朗日插值

希望你還沒還給你的老師

插值多項式唯一性定理

設 $ f(x) = 3x^2+5x+7$

任取三個點

$f(1)=15$

$f(2)=24$

$f(3)=49$

拉格朗日插值

$f(x)=$

$15\times \frac{(x-2)(x-3)}{(1-2)(1-3)}+\\24\times \frac{(x-1)(x-3)}{(2-1)(2-3)}+\\49\times\frac{(x-1)(x-2)}{(3-1)(3-2)}=\\3x^2+5x+7$

插值多項式唯一性定理

$ f(x) = 3x^2+5x+7$

$f(1)=15$

$f(2)=24$

$f(3)=49$

插值法

取樣

再談多項式乘法

$f_C(x)=f_A(x)f_B(x)$

如果 $A,B$ 都是 $n$ 次多項式，那 $C$ 就是 $2n$ 次多項式

透過 $f_A$ , $f_B$ 採樣 $f_C$ 的 $2n+1 個點$

利用插值法拼回 $f_C$ 的樣子

就能避免樸素多項式乘法了!

那怎麼做取樣 / 差值 ?

離散傅立葉變換

$ f(x) = 3x^2+5x+7$

$f(\omega)=a$

$f(\omega^2)=b$

$f(\omega^3)=c$

逆離散傅立葉變換

離散傅立葉變換

選擇 "夠好的點" 採樣，讓採樣插值都很快!

夠好的點 - 單位根 $\omega$

$n$ 次單位根是 $x^n=1$ 的所有複數 (有 $i$ 的那個) 解
高中數學範圍內教過 3 次原根

$$1, \frac{1+\sqrt{3}i}{2},\frac{1-\sqrt{3}i}{2} $$

單位根 $\omega$

$n$ 次單位根剛好可以表示為

在複數平面上的一個正 $n$ 邊形，

頂點都在半徑為 1
原點在 $0+0i$ 上的圓

其中一個頂點在 $1+0i$

單位根 $\omega$

$n$ 次單位根剛好可以表示為

在複數平面上的一個正 $n$ 邊形，

頂點都在半徑為 1
原點在 $0+0i$ 上的圓

其中一個頂點在 $1+0i$

利用圖形的關係，來得到單位根的數值 !

$\omega = \exp(\frac{2\pi}{n}i)=\cos{\frac{2\pi}{n}}+i\sin\frac{2\pi}{n}$

離散傅立葉變換

對一個 $n-1$ 次多項式 $f(x)=a_0+a_1x^1+\cdots+a_{n-1}x^{n-1}$，

求 $f(1), f(\omega), f(\omega^2), f(\omega^3)\cdots ,f(\omega^{n-1})$

其中 $\omega$ 是 $n$ 次原根

{\begin{bmatrix} 1&1&1&1&\cdots &1\\ 1&\omega &\omega ^{2}&\omega ^{3}&\cdots &\omega ^{N-1}\\ 1&\omega ^{2}&\omega ^{4}&\omega ^{6}&\cdots &\omega ^{2(N-1)}\\ 1&\omega ^{3}&\omega ^{6}&\omega ^{9}&\cdots &\omega ^{3(N-1)}\\ \vdots &\vdots &\vdots &\vdots &\ddots &\vdots \\ 1&\omega ^{N-1}&\omega ^{2(N-1)}&\omega ^{3(N-1)}&\cdots &\omega ^{(N-1)(N-1)} \end{bmatrix}} {\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{bmatrix}}= {\begin{bmatrix} f(1) \\ f(\omega^1) \\ f(\omega^2) \\ \vdots \\ f(\omega^{n-1}) \end{bmatrix}}

對一個 $n-1$ 次多項式 $f(x)=a_0+a_1x^1+\cdots+a_{n-1}x^{n-1}$，

求 $f(1),f(\omega), f(\omega^2), f(\omega^3)\cdots,f(\omega^{n-1})$

利用矩陣表示前面的運算

逆離散傅立葉變換

由 $f(1), f(\omega^1), f(\omega^2)\cdots ,f(\omega^{n-1})$ 求出

一個 $n-1$ 次多項式 $f(x)=a_0+a_1x^1+\cdots+a_{n-1}x^{n-1}$
的所有係數

基本上什麼詭異差值法都很麻煩
可以利用前一頁的式子
把這個問題變成乘上那一大團東西的反矩陣

所以這一團東西的反矩陣是 ???

\frac{1}{N}

-1

證明：乘起來會變單位矩陣...

那東西的反矩陣是 ...

using CBigNum = vector<complex<double>>;
CBigNum btoc(const BigNum &b, int N) {
    CBigNum c;
    for(int i:b) c.emplace_back(i);
    while(c.size() < N) c.emplace_back(0);
    return c;
}

BigNum ctob(const CBigNum &c) {
    BigNum b;
    for(auto i:c) b.emplace_back((int)round(i.real()));
    return b;
}

使用該方法有小數精確度問題，要小心

BigNum mul2(const BigNum &a_, const BigNum &b_)
{
    int N = a_.size() + b_.size();

    CBigNum a = btoc(a_, N);
    CBigNum b = btoc(b_, N);

    auto fft = [=](const CBigNum &a, bool inv)-> CBigNum {
		/*next page*/
    };

    a = fft(a, false);
    b = fft(b, false);

    for(int i=0;i<N;++i)
        a[i] *= b[i];

    a = fft(a, true);
    auto res = ctob(a);
    return trim(advance(res));
}

auto fft = [=](const CBigNum &a, bool inv)-> CBigNum {
    //C++20, or PI = acos(-1)
    constexpr double PI = std::numbers::pi; 
    constexpr auto I = complex<double>(0, 1);

    int flag = inv ? -1 : 1;

    CBigNum res(a.size());
    for(int i=0;i<N;++i)
        for(int j=0;j<N;++j)
            res[i] += a[j] * exp( flag* i*j*(2*PI/N) * I);

    if (inv) for(auto &c:res) c/=N;

    return res;
};

就按照前面公式乘一乘，加一加

離散傅立葉變換

基本上來說，離散傅立葉/逆離散傅立葉變換幾乎一模一樣
矩陣乘法複雜度 $ O(n\times n \times 1) = O(n^2)$

完全沒有變好!!

但是這個矩陣夠怪，應該有可以加速的地方

快速傅立葉變換 FFT

在 $O(n\log n)$ 的時間乘上那一個奇怪的矩陣

快速傅立葉變換 FFT

假設 $f(x)$ 有 $2n$ 項
考慮 $f(w^x)=a_0+a_1\omega^x+a_2\omega^{2x}+\cdots+a_{2n-1}\omega^{(2n-1)x}$

把奇偶項分開

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$a_1\omega^{x}+a_3\omega^{3x}+a_5\omega^{5x}+\cdots+a_{2n-1}\omega^{(2n-1)x}$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

這個 $\omega$ 是 $2n$ 次原根

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

單位根的性質 1

2n 次單位根的偶數次方 $\omega^0, \omega^2, \omega^4...$

恰等於 n 次單位根的所有組合

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

單位根的性質 1

2n 次單位根的偶數次方 $\omega^0, \omega^2, \omega^4...$

恰等於 n 次單位根的所有組合

因此預先遞迴計算

$f_{even}(x)=a_0+a_2x+a_4x^2+a_{2n-2}{x^{n-1}}$ 的傅立葉變換

$f(1), f(\omega_{n}^1), f(\omega_{n}^2)\cdots f(\omega_{n}^{n-1})$

就能使用其結果組合答案 !

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

預先遞迴計算

$f_{even}(x)=a_0+a_2x+a_4x^2+a_{2n-2}{x^{n-1}}$ 的傅立葉變換 $f_{even}(\omega_n^x)$

$f_{odd}(x)=a_1+a_3x+a_5x^2+a_{2n-1}{x^{n-1}}$ 的傅立葉變換 $f_{odd}(\omega_n^x)$

那 $f(\omega_{2n}^x)=$

$f_{even}(w_{2n}^{2x})+\omega_{2n}^xf_{odd}(w_{2n}^{2x})$

$=f_{even}(w_{n}^x)+\omega_{2n}^xf_{odd}(w_{n}^x)$

對於所有 $x < n$ 成立

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

$f(\omega_{2n}^x)=f_{even}(w_{n}^x)+\omega_{2n}^xf_{odd}(w_{n}^x)$

對於所有 $x < n$ 成立

$f(\omega_{2n}^{x+n})=f_{even}(w_{2n}^{2(x+n)})+\omega_{2n}^{x+n}f_{odd}(w_{2n}^{2(x+n)})$

$=f_{even}(w_{n}^{x+n})+\omega_{2n}^{x+n}f_{odd}(w_{n}^{x+n})$

根據單位根的定義，$n$ 次單位根 $\omega_n$ 是 $x^n=1$ 的解 $(\omega_n^n=1)$

$=f_{even}(w_{n}^{x}\times 1)+\omega_{2n}^{x+n}f_{odd}(w_{n}^{x}\times 1)$

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

$f(\omega_{2n}^x)=f_{even}(w_{n}^x)+\omega_{2n}^xf_{odd}(w_{n}^x)$

對於所有 $x < n$ 成立

$f(\omega_{2n}^{x+n})=f_{even}(w_{n}^{x})+\omega_{2n}^{x+n}f_{odd}(w_{n}^{x})$

進一步的 $\omega^n_{2n}=\omega_{2}=-1$

$=f_{even}(w_{n}^{x})-\omega_{2n}^{x}f_{odd}(w_{n}^{x})$

快速傅立葉變換 FFT

$f(w^x)=a_0+a_2\omega^{2x}+a_4\omega^{4x}+\cdots+a_{2n-2}\omega^{(2n-2)x}+$

$\omega^{x}(a_1+a_3\omega^{2x}+a_5\omega^{4x}+\cdots+a_{2n-1}\omega^{(2n-2)x})$

$f(\omega_{2n}^x)=f_{even}(w_{n}^x)+\omega_{2n}^xf_{odd}(w_{n}^x)$

對於所有 $x < n$ 成立

$f(\omega_{2n}^{x+n})=f_{even}(w_{n}^{x})-\omega_{2n}^{x}f_{odd}(w_{n}^{x})$

於是所有的 $f(\omega^x)$ 就都能用 $f_{even}, f_{odd}$ 表示出來了!

如果令輸入長度是 $n=2^t$ ，就能完美的不斷遞迴此公式 !

BigNum mul3(const BigNum &a_, const BigNum &b_)
{
    int N = a_.size() + b_.size();
    int p2 = 1;
    while (p2 < N) p2*=2; 
    N = p2; // 要把 N 補成 2 的次方

    CBigNum a = btoc(a_, N);
    CBigNum b = btoc(b_, N);

	//... 下面跟前面一樣， fft_r 下頁解說
    a = fft_r(a, false, N);
    b = fft_r(b, false, N);
    //...
}

CBigNum fft_r(const CBigNum &a, bool inv, int N, bool first = true) {
    constexpr double PI = std::numbers::pi; //C++20, or acos(-1)
    constexpr auto I = complex<double>(0, 1);
    int flag = inv ? -1 : 1;

    CBigNum res(a);
    if (N==1) 
        return a;

    CBigNum even, odd;
    for(int i=0;i<N;++i) {
        if (i % 2 == 0) even.emplace_back(a[i]);
        else odd.emplace_back(a[i]);
    }

    even = fft_r(even, inv, N/2, false);
    odd  = fft_r(odd , inv, N/2, false);

    for(int i=0;i<N/2;++i) {
        res[i]     = even[i] + odd[i] * exp( flag * i * (2*PI/N) * I);
        res[i+N/2] = even[i] - odd[i] * exp( flag * i * (2*PI/N) * I);
    }

    if (inv&&first) for(auto &c:res) c/=N;
    return res;
}

快速傅立葉變換

快速傅立葉變換每次將輸入分解成大小為一半的兩個子問題
花費時間 $T(n)=2T(n/2)+O(n)$
根據主定理，$O(T(n)=O(n\log n)$

於是我們可以在 $O(n\log n)$的時間來處理大數乘法了!

不過常數有點大 !?

快速傅立葉變換 - 迴圈版

透過分析遞迴的過程，來把遞迴的程式改寫成迴圈的
當中關鍵點在 "奇偶分割" 的步驟

for(int i=0;i<N;++i) {
    if (i % 2 == 0) even.emplace_back(a[i]);
    else odd.emplace_back(a[i]);
}

快速傅立葉變換 - 迴圈版

我們可以重複使用陣列來擺放資料，讓要處理的資料連續就好

0	1	2	3	4	5	6	7

n = 8

n = 4

n = 2

0	4	2	6	1	5	3	7

n = 1

0	4	2	6	1	5	3	7

0	2	4	6	1	3	5	7

快速傅立葉變換 - 迴圈版

0	1	2	3	4	5	6	7

n = 8

0	4	2	6	1	5	3	7

n = 1

要用迴圈的話，要從 n=1 開始計算

問題是 n = 1 的那一層資料是怎麼放的 ?

方法 1. 模擬交換一遍 $O(n\log n)$

方法 2. 根據觀察
把第 $x$ 個資料放在第 $\operatorname{bit\_reverse}(x)$ 個格子上面

bit_reverse

把一個長度為 N 的數字二進位前後顛倒
像是 110010 變成 010011
C++ 沒有內建這東西，不難寫，看想寫得多花式...

unsigned int bit_rev_simple(unsigned int a, int N) {
    unsigned int ans = 0;
    for(int i=0;i<N;++i) {
        ans = ans*2 + a%2;
        a/=2;
    }
    return ans;
}
unsigned int bit_rev(unsigned int a,int N) {
    a=((a&0x55555555U)<<1)|((a&0xAAAAAAAAU)>>1);
    a=((a&0x33333333U)<<2)|((a&0xCCCCCCCCU)>>2);
    a=((a&0x0F0F0F0FU)<<4)|((a&0xF0F0F0F0U)>>4);
    a=((a&0x00FF00FFU)<<8)|((a&0xFF00FF00U)>>8);
    a=((a&0x0000FFFFU)<<16)|((a&0xFFFF0000U)>>16);
    return a>>(32-N);
}

CBigNum fft(const CBigNum &a, bool inv, int N) {
    constexpr double PI = std::numbers::pi;
    constexpr auto I = complex<double>(0, 1);
    int flag = inv ? -1 : 1;

    int bit_len = __lg(N); // log2
    CBigNum in(N);
    CBigNum out(N);

    for(int i=0;i<N;++i)
        in[bit_rev(i, bit_len)] = a[i];

    /* next page */

    if (inv) for(auto &c:in) c/=N;
    return in;
}

for(int n=2 ; n<=N ; n*=2)
{
    for(int base=0 ; base<N ; base+=n)
    {
        for(int i=0;i<n/2;++i)
        {
            auto wi = exp( flag * i * (2*PI/n) * I);
            out[base+i]     = in[base+i] + in[base+i+n/2] * wi;
            out[base+i+n/2] = in[base+i] - in[base+i+n/2] * wi;
        }
    }
    in.swap(out);
}

這個三層迴圈還有一些實作細節可以進一步優化，參見卦長模板

大數乘法投影片完整範例 : gist

練習題

https://asiahk16prel.kattis.com/problems/aplusb
hdu 1402
(C) https://codeforces.com/gym/100783/attachments
hdu 4609

進階主題 : 數論轉換，可以避免 FFT 的浮點誤差