解題區/解題區 - 其他題目

[Tmt514] Beverage Cup 2 Warmup C - Exploding CPU 2

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution

Problem

詢問一個區間有多少 鄒賽爾數 (Zeisel number)

鄒賽爾數是一種無平方數因數的數，而且至少有三個質因數可以用下式表示，質因數由小排到大滿足$p_{x} = A p_{x-1} + B$。

105, 1419, 1729, 1885, 4505, 5719, 15387, 24211, 25085, …

1
2
3
4
5
6
7
4505 = 1 * 5 * 17 * 53
A = 3, B = 2
p0 = 1
p1 = 3 * p0 + 2 = 5
p2 = 3 * p1 + 2 = 17
p3 = 3 * p2 + 2 = 53

找到 $10^{15}$ 以內的所有 Zeisel number。

Sample Input

1
2
3
2
4505 4505
0 5000

Sample Output

1
2
1
5

Solution

一開始設想，既然小於 $10^{15}$ 又要三個質數，那麼最小的質因數應該要小於 $10^5$，否則會超過指定範圍。

接著第一次嘗試時，窮舉變數 $A$ 針對第一個質因數$p_{1}$ 得到$B = p_{1} - A$，得到下一項在範圍內 … 這樣效率非常差。A 的範圍非常大，質因數也非常多，而且還要檢驗每一項是否質數，效能相當低落。

接著第二次嘗試時，窮舉$p_{1}, p_{2}$，這樣的計算量是 $10^5$ 內的質數個數 n，至少為 $O(n^2)$，嘗試一下發現速度還是太慢，藉由聯立方程式解出 A, B，卻有不少 A 是不存在的整數解。

1
2
3
A + B = p1
p1 A + B = p2
A = (p2 - p1)/(p1 - 1), B = p1 - A

接著第三次嘗試時，反過來窮舉，先讓 $A$ 一定有解，再看$p_{2}$ 是否為質數，因此$p_2 = p_1 + k \times (p_1 - 1)$，由於$p_{i} - 1$ 造成的增長幅度遠比質數個數快很多，效能上有顯著地提升。

現在還有一個疑問，保證$p_1 \le 10^5$，但是當$p_1$ 很小$p_2$ 可能會大於 $10^5$ 嗎，甚至必須窮舉到 $\sqrt{10^{15}}$？

這裡我不確定，放大範圍去嘗試時，發現找到的 Zeisel number 最多 9607 個，放大搜索範圍不再增加，那就不斷地縮減，讓找到的時間縮小不逾時。

1.980s TLE

不是說好限制 2.000s？1.980s 得到 TLE，第一次使用 NTUJ 真歡樂。

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
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <queue>
using namespace std;
#define MAXL (10000000>>5)+1
#define GET(x) (mark[x>>5]>>(x&31)&1)
#define SET(x) (mark[x>>5] |= 1<<(x&31))
int mark[MAXL];
int P[5500000], Pt = 0;
vector<long long> ret;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 10000000;
    for(i = 2; i <= n; i++) {
        if(!GET(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET(j);
            P[Pt++] = i;
        }
    }
}
long long MAXVAL = 1e+15;
int isprime(long long v) {
    if (v < 2)	
        return 0;
    if (v < 10000000)
        return !GET(v);
    for (int i = 0; i < Pt && (long long) P[i] * P[i] <= v; i++)
        if (v%P[i] == 0)
            return 0;
    return 1;
}
void make() {
    for (int i = 0; i < Pt && P[i] < 100000; i++) {
        for (int j = P[i] + P[i]-1; j < 1000000; j += P[i] - 1) {
            if (!isprime(j))
                continue;
            long long A = (j - P[i])/(P[i]-1);
            long long B = P[i] - A;
            long long m = P[i], cnt = 1, mm = P[i];
            while (mm <= MAXVAL) {
                if (A && MAXVAL / m / A <= 0)
                    break;
                if (m * A + B <= m)
                    break;
                m = m * A + B;
                if (MAXVAL / mm / m <= 0)
                    break;
                if (!isprime(m))
                    break;
//				printf("%lld ", m);
                mm = mm * m, cnt ++;
                if (cnt >= 3) {
                    ret.push_back(mm);
//					printf("%lld %lld %lld\n", mm, A, B);
                }
            }
        }
    }
    sort(ret.begin(), ret.end());
    for (int i = 0; i < 243; i++)
        printf("%d\n", ret[i]);
//	printf("%d\n", ret.size());
}
int main() {
    sieve();
    make();
    int testcase;
    long long L, R;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%lld %lld", &L, &R);
        int r = 0;
        for (int i = 0; i < ret.size(); i++)
            r += ret[i] >= L && ret[i] <= R;
        printf("%d\n", r);
    }
    return 0;
}
/*
2
4505 4505
0 5000
*/