[BOJ 4798] 등차수열에 관한 디리클레의 정리

2021. 6. 2. 01:19PS/백준

굉장히 어려운 문제이다...

코드업에 비슷한 문제를 냈었는데 그거보다 어렵다 ㅜㅜ

하지만 무려 8일간 고민을 거듭해서 맞았다!

그 풀이에 대해 적어보려고 한다!

 

 

문제 

등차수열에 관한 디리클레의 정리는 서로소인 두 양의 정수 a와 b가 있을 때, 등차수열 t(n) = a*n + b (n ≥ 0)은 무한히 많은 소수를 포함한다는 내용이다. 소수는 2보다 큰 양의 정수로, 약수가 1과 자기자신 밖에 없는 수이다.

예를 들어, a=4, b=3인 경우 등차수열은 다음과 같다.

3, 7, 11, 15, 19, 23, 27, 31, 35, ...,

여기서 이 등차수열의 첫 부분에 많은 소수가 있음을 눈으로 볼 수 있다.

a > 0과 b ≥ 0, 그리고 U ≥ L ≥ 0이 주어졌을 때, t(n) = a*n+b에 소수가 몇 개 있는지 구하는 프로그램을 작성하시오.

(L ≤ n ≤ U)

 

 

1 단계 : $a, b$에 따른 $an+b$꼴 소수의 존재 가능성에 대해 생각하기

$gcd(a,b)=d>1$라고 하자

그럼 $an+b \equiv 0 \pmod{d}$

아쉽게도 소수가 아니다.

 

라고 하면 '틀렸습니다.'를 받는다.

$n=0$이면 $a$항이 의미가 없어진다!

따라서 이때는 $b$가 소수인지 꼭! 확인하자.

 

2 단계 : 소수 찾기

$n$ 자리에 들어올 숫자는 최대 '100만개'이다.

이건 배열에 저장해두고 소수를 에라토스테네스의 체로 거르자.

 

2.1 단계 : 에라토스테네스 체 구현

$an+b$ 꼴 수의 범위는 $10^{12}$까지이다.

그러니 $10^{6}$까지의 소수 배열을 만들도록 하자.

 

$an+b \equiv 0 \pmod{p}$를 만족하는 $n$을 찾는게 우리의 목표이다.

이것은 유클리드 호제법으로 찾자!

그럼 $n \equiv m \pmod{p}$을 얻는다.

위 식을 만족하는 $n$은 $an+b$가 합성수가 되도록 한다.

위의 $n$에 해당하는 수들을 모두 쳐낸다!

 

유클리드 호제법을 사용할 때 주의점이 있다.

만약 $p\ |\ a$이면 $n$항이 의미가 없어지므로 $p \not |\ a$인 경우에만 사용하자.

 

반면 $an+b=p$인 경우를 살펴야한다.

이를 만족하려면 $n$이 주어진 범위에서 가장 작은 수가 되어야한다.

즉, $n$이 *$L/p \cdot p+m$일 때 소수인지 살핀다.

 

*반면 $a=1$인 경우

$L/p \cdot p+p+m$도 살펴야한다.

이유를 알 수 없다.

 

3 단계 : 소수 판정

소수인지 의심되는 수를 확정적으로 소수인지 알아내자.

밀러-라빈 소수 판정법과 분할 정복을 이용한 거듭제곱 계산을 하면 된다!

 

소스 코드는 다음과 같다.

 

더보기

# 코드 보기

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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#include<stdio.h>
#include<math.h>
#include<string.h>
 
char s[1000003];
char rs[1000001];
static const int max = 1000000;
 
long long power(long long n, long long k, long long mod){
    long long s = 1;
    for( ; k; k >>= 1){
        if(k & 1) s = (__int128)s * n % mod;
        n = (__int128)n * n % mod;
    }
    return s;
}
 
int isPrime(long long n) {
    static const long long aPrimes[12= {23571113171923293137};
    static const long long bPrimes[12= {20471373653253260013215031751215230289874734747496603833415500717283210382512305654641305100, 1ll<<63};
    if(n <= 1return 0;
    for(int i = 0; i < 12; i++)
        if(n % aPrimes[i] == 0)
            return n == aPrimes[i];
    long long m = n - 1, s = __builtin_ctzll(m), d = m >> s, k;
    for(k = 0; n >= bPrimes[k]; k++);
    for(long long i = 0; i <= k; i++){
        long long x = power(aPrimes[i], d, n), r = 0
        if(x != 1) {
            for(r = 0; r < s && x != m; r++)
                x = ((__int128)x * x) % n;
            if(r == s)
                return 0;
        }
    }
    return 1;
}
 
long long gcd(long long a, long long b){return b == 0? a : gcd(b, a % b);}
 
long long EEA(long long a, long long p, long long b){
    long long r1 = a, r2 = p, r, q;
    long long s1 = 1, s2 = 0, s;
    long long t1 = 0, t2 = 1, t;
    
    while(r2 > 0){
        q = r1 / r2;
        r = r1 - q * r2;
        r1 = r2;
        r2 = r;
 
        s = s1 - q * s2;
        s1 = s2;
        s2 = s;
 
        t = t1 - q * t2;
        t1 = t2;
        t2 = t;
    }
    
    b =- b * s1 % p;
    if(b < 0) b += p;
    return b;
}
 
int main(){
    long long a, b, L, U, st;
    long long t = 1, i, j, q, c = 0, mod;
    
    for(i = 3; i <= max; i += 2){
        if(s[i]) continue;
        for(j = i * i; j <= max; j += i)
            s[j] = 1;
    }
     
    while(1){
        scanf("%lld"&a);
        if(!a) return 0;
        scanf("%lld %lld %lld"&b, &L, &U);
        
        if(gcd(a, b) > 1) {
            if(L == 0 && isPrime(b)) c = 1;
            printf("Case %lld: %lld\n", t++, c);
            c = 0;
            continue;
        }
        
        q = (int)sqrt(U * a + b) + 1;
        
        for(i = 2; i <= q; i += 2){
            if(a % i == 0 || s[i]) {
                if(i == 2) i--;
                continue;
            }
            
            
            mod = EEA(a, i, b);
            st = L / i * i + mod;
            if(st < L) st += i;
            
            for(j = st; j <= U; j += i) 
                rs[j - L] = 1;
            
            rs[st - L] = !isPrime(st * a + b);
            if(a == 1 && st + i <= U) rs[st + i - L] = !isPrime((st + i) * a + b);
 
            if(i == 2) i--;
        }
        
        rs[0= !isPrime(L * a + b);
        rs[1= !isPrime((L + 1* a + b);
        
        for(i = L; i <= U; i++
            c += !rs[i - L];
        
        printf("Case %lld: %lld\n", t++, c);
        memset(rs, 0, max + 1);
        c = 0;
    }
}
cs

# 닫기

 

이 코드는 아슬아슬하게 통과한다..

고수분들 코드를 보고 공부해야지

'PS > 백준' 카테고리의 다른 글

[BOJ 1770] 배수와 약수의 개수  (7) 2021.08.21
[BOJ 11778] 피보나치 수와 최대공약수  (1) 2021.08.21
[BOJ 10908] Phibonacci  (0) 2021.08.21
[BOJ 19577] 수학은 재밌어  (0) 2021.08.16
[BOJ 13430] 합 구하기  (2) 2021.08.15