bzoj4802 歐拉函數（附Millar-Rabin和Pollard-Rho講解）

傳送門：http://www.lydsy.com/JudgeOnline/problem.php?id=4802

【題解】

參考：http://www.matrix67.com/blog/archives/234

Millar-Rabin質數檢驗方法：

根據費馬小定理，如果p是素數，a<p，那么有a^(p-1) mod p = 1。

直觀想法我們直接取若干個a，如果都有一個不滿足，那么p就是合數。

遺憾的是，存在Carmichael數：你無論取多少個a，有一個不滿足，算我輸。

比如：561 = 11*51就是一個Carmichael數。

那么就很江了啊。。我們需要改進算法。

首先有：如果p是素數，x是小於p的正整數，且x^2 mod p = 1，那么要么x=1，要么x=p-1

（這個廢話，x=p-1模意義下等於x=-1）

然后我們可以展示下341滿足2^340 mod 341 = 1，卻不是素數（341=31*11）的原因：

2^340 mod 341 = 1

2^170 mod 341 = 1

2^85 mod 341 = 32

（32這個數很py啊怎么不等於340也不等於1啊。。這明顯有交易嘛）

那么就能說明這個數不是素數。

如果是素數，一定是從p-1變到1，或是把所有2的次冪去除完，本來就等於1（這樣平方完就一直是1了）

所以要么把所有2的次冪去除完，本來就等於1，要么存在某一個次冪=p-1（這樣就正常多了）

這就是Millar-Rabin素數驗證的二次探測。

應該來說Millar-Rabin算法也是挺好寫的

其中mul(a,b,c)表示a*b%c（因為a*b會爆longlong，所以用快速加）

namespace Millar_Rabin {
    const int Prime[14] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};    
    const int PN = 13;
    
    inline bool witness(int pr, ll res, int times, ll n) {
        ll p = pwr2((ll)pr, res, n);
        if(p == 1) return 0;   
        for (int i=0; i<times; ++i) {
            if(p == n-1) return false;
            if(p == 1) return false;
            p = mul(p, p, n);
        }
        return true;
    }
    
    inline bool main(ll n) {
        for (int i=1; i<=PN; ++i) {
            if(n == Prime[i]) return 1;
            if(n % Prime[i] == 0) return 0;
        }
        ll p = n-1;
        int times = 0;
        while(!(p&1)) {
            ++times;
            p >>= 1;
        }
        for (int i=1; i<=PN; ++i) 
            if(witness(Prime[i], p, times, n)) return false;
        return true;
    }
}

然后我們會檢驗素數了，現在要質因數分解。

好了下一個是Pollard-Rho算法：

如果現在拆分的是n：Pollard-Rho(n)

主要流程：Millar-Rabin判斷是否質數，是返回，否就試圖找出其中一個因子d，然后遞歸做Pollard-Rho(d)和Pollard-Rho(n/d)。

我猜你會說廢話這誰都會。問題在於：試圖找出其中一個因子d

參考：https://wenku.baidu.com/view/3db5c7a6ad51f01dc381f156.html?from=search

參考文章講的非常詳細了。。我就不細講了qwq

所以這題就是分解因數，按照歐拉函數定義式求解即可。

# include <stdio.h>
# include <string.h>
# include <iostream>
# include <algorithm>
// # include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
const int M = 5e5 + 10;
const int mod = 1e9+7;

# define RG register
# define ST static

inline ll mul(ll a, ll b, ll mod) {
    ll ret = 0;
    a %= mod, b %= mod;
    while(b) {
        if(b&1) {
            ret = ret + a;
            if(ret >= mod) ret -= mod;
        }
        a <<= 1;
        if(a >= mod) a -= mod;
        b >>= 1;
    }
    return ret;
}

inline ll pwr2(ll a, ll b, ll mod) {
    ll ret = 1;
    a %= mod;
    while(b) {
        if(b&1) ret = mul(ret, a, mod);
        a = mul(a, a, mod);
        b >>= 1;
    }
    return ret;
}

inline ll gcd(ll a, ll b) {
    return b==0 ? a : gcd(b, a%b);
}

namespace Millar_Rabin {
    const int Prime[14] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};    
    const int PN = 13;
    
    inline bool witness(int pr, ll res, int times, ll n) {
        ll p = pwr2((ll)pr, res, n);
        if(p == 1) return 0;   
        for (int i=0; i<times; ++i) {
            if(p == n-1) return false;
            if(p == 1) return false;
            p = mul(p, p, n);
        }
        return true;
    }
    
    inline bool main(ll n) {
        for (int i=1; i<=PN; ++i) {
            if(n == Prime[i]) return 1;
            if(n % Prime[i] == 0) return 0;
        }
        ll p = n-1;
        int times = 0;
        while(!(p&1)) {
            ++times;
            p >>= 1;
        }
        for (int i=1; i<=PN; ++i) 
            if(witness(Prime[i], p, times, n)) return false;
        return true;
    }
}

namespace PollardRho {
    const int N = 110;
    ll q[N]; int qn;    
    
    inline void PR(ll n) {
        if(Millar_Rabin::main(n)) {
            q[++qn] = n;
            return ;
        }
        ll a, b, c, del;
        while(1) {
            c = rand() % n;
            a = b = rand() % n;
            b = (mul(b, b, n) + c) % n;
            while(a != b) {
                del = a-b;
                del = gcd(abs(del), n);
                if(del > 1 && del < n) {
                    PR(del); PR(n/del);
                    return ;
                }
                a = (mul(a, a, n) + c) % n;
                b = (mul(b, b, n) + c) % n;
                b = (mul(b, b, n) + c) % n;
            }
        }
    }
    
    inline ll getphi(ll n) {
        if(n == 1) return 1ll;
        sort(q+1, q+qn+1);
        ll res = q[1] - 1;
        for (int i=2; i<=qn; ++i) {
            if(q[i] != q[i-1]) res = res * (q[i] - 1);
            else res = res * q[i];
        }
        return res;
    }
    
    inline void main(ll n) {
        qn = 0;
        PR(n);
        cout << getphi(n) << endl;
    }
}

int main() {
    srand(19260817);
    ll n; cin >> n;
    if(n == 1) {
        puts("1");
        return 0;
    }
    PollardRho::main(n);
    return 0;
}