[Algorithm] BOJ 13926 gcd(n, k) = 1 — φ(n) 계산 (C++/Python)

문제

문제: gcd(n, k) = 1 (BOJ 13926)
링크: https://www.acmicpc.net/problem/13926
설명: 자연수 n이 주어졌을 때, gcd(n, k) = 1을 만족하는 1 ≤ k ≤ n의 개수를 구합니다. 이는 오일러 파이 함수 φ(n)입니다.
제한: 1 ≤ n ≤ 10^18, 시간 2초, 메모리 512MB

아이디어

φ(n) = n × ∏(1 − 1/p) (p는 n의 서로 다른 소인수)
10^18까지 처리하려면 빠른 소인수분해가 필요합니다 → Miller–Rabin(64비트 결정론 세트) + Pollard Rho
분해된 서로 다른 소인수 집합을 이용하여 φ(n)을 계산합니다.

복잡도

소수 판정/분해: 보통 빠름 (경우에 따라 준선형), 구현상 평균적으로 10^18 범위에서도 충분
메모리: O(log n) 수준

정답 코드

C++

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
// 더 많은 정보는 https://42jerrykim.github.io 에서 확인하세요.
#include <bits/stdc++.h>
using namespace std;
using u128 = unsigned __int128;
using u64 = unsigned long long;
using u32 = unsigned int;

static inline u64 mul_mod(u64 a, u64 b, u64 m) {
    return (u128)a * b % m;
}

static inline u64 pow_mod(u64 a, u64 e, u64 m) {
    u64 r = 1;
    while (e) {
        if (e & 1) r = mul_mod(r, a, m);
        a = mul_mod(a, a, m);
        e >>= 1;
    }
    return r;
}

static bool is_prime(u64 n) {
    if (n < 2) return false;
    for (u64 p : {2ull,3ull,5ull,7ull,11ull,13ull,17ull,19ull,23ull,29ull,31ull,37ull}) {
        if (n % p == 0) return n == p;
    }
    u64 d = n - 1, s = 0;
    while ((d & 1) == 0) { d >>= 1; ++s; }
    // Deterministic bases for 64-bit
    for (u64 a : {2ull, 325ull, 9375ull, 28178ull, 450775ull, 9780504ull, 1795265022ull}) {
        if (a % n == 0) continue;
        u64 x = pow_mod(a % n, d, n);
        if (x == 1 || x == n - 1) continue;
        bool comp = true;
        for (u64 r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1) { comp = false; break; }
        }
        if (comp) return false;
    }
    return true;
}

static u64 rho(u64 n) {
    if ((n & 1ull) == 0) return 2;
    std::mt19937_64 rng((u64)chrono::high_resolution_clock::now().time_since_epoch().count());
    while (true) {
        u64 c = (rng() % (n - 2)) + 1;
        u64 x = (rng() % (n - 2)) + 2;
        u64 y = x;
        u64 d = 1;
        auto f = [&](u64 v){ return (mul_mod(v, v, n) + c) % n; };
        while (d == 1) {
            x = f(x);
            y = f(f(y));
            u64 diff = x > y ? x - y : y - x;
            d = std::gcd(diff, n);
        }
        if (d != n) return d;
    }
}

static void factor(u64 n, vector<u64>& fac) {
    if (n == 1) return;
    if (is_prime(n)) { fac.push_back(n); return; }
    u64 d = rho(n);
    factor(d, fac);
    factor(n / d, fac);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    unsigned long long n; if (!(cin >> n)) return 0;
    if (n == 0) { cout << 0 << '\n'; return 0; }
    if (n == 1) { cout << 1 << '\n'; return 0; }

    vector<u64> fac; fac.reserve(64);
    factor(n, fac);
    sort(fac.begin(), fac.end());
    fac.erase(unique(fac.begin(), fac.end()), fac.end());

    __int128 phi = n;
    for (u64 p : fac) {
        phi = phi / p * (p - 1);
    }
    unsigned long long ans = (unsigned long long)phi;
    cout << ans << '\n';
    return 0;
}

Python

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# 더 많은 정보는 https://42jerrykim.github.io 에서 확인하세요.
import sys
import random
from math import gcd

def mul_mod(a, b, mod):
    return (a * b) % mod

def pow_mod(a, e, mod):
    r = 1
    a %= mod
    while e:
        if e & 1:
            r = (r * a) % mod
        a = (a * a) % mod
        e >>= 1
    return r

def is_prime(n: int) -> bool:
    if n < 2:
        return False
    small_primes = [2,3,5,7,11,13,17,19,23,29,31,37]
    for p in small_primes:
        if n % p == 0:
            return n == p
    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        if a % n == 0:
            continue
        x = pow_mod(a, d, n)
        if x == 1 or x == n - 1:
            continue
        comp = True
        for _ in range(s - 1):
            x = (x * x) % n
            if x == n - 1:
                comp = False
                break
        if comp:
            return False
    return True

def rho(n: int) -> int:
    if n % 2 == 0:
        return 2
    while True:
        c = random.randrange(1, n - 1)
        x = random.randrange(2, n - 1)
        y = x
        d = 1
        def f(v):
            return (v * v + c) % n
        while d == 1:
            x = f(x)
            y = f(f(y))
            d = gcd(abs(x - y), n)
        if d != n:
            return d

def factor(n: int, fac: list):
    if n == 1:
        return
    if is_prime(n):
        fac.append(n)
        return
    d = rho(n)
    factor(d, fac)
    factor(n // d, fac)

def main():
    data = sys.stdin.read().strip().split()
    if not data:
        return
    n = int(data[0])
    if n == 0:
        print(0)
        return
    if n == 1:
        print(1)
        return
    fac = []
    factor(n, fac)
    fac = sorted(set(fac))
    phi = n
    for p in fac:
        phi = phi // p * (p - 1)
    print(phi)

if __name__ == "__main__":
    main()