luogu 1463 - HAOI2007 - 反素數 - 打表，數論

題目描述：

對於任何正整數x，其約數的個數記做g(x)。例如g(1)=一、g(6)=4。

若是某個正整數x知足：g(x)>g(i) 0<i<x，則稱x爲反質數。例如，整數1，2，4，6等都是反質數。

如今給定一個數N，你能求出不超過N的最大的反質數麼？

N<=2e9。

（暴力出奇跡(´థ౪థ)σ）

不打表的AC作法見洛谷題解區。

設將X質因數分解後爲X=p1^n1 * p2^n2 * ... * pk^nk，那麼X的約數個數爲(n1 + 1) * (n2 + 1) * ... * (nk + 1)。

這道題若是要打表的話，最暴力的算法是：枚舉for X = 1 to 2e9，將X質因數分解後統計約數個數。但這樣作的話5個小時恐怕都打不出表來。

不妨設X=p^n * ...，其中p是某個質數，省略號部分是其餘質因子。不可貴出g(X) = (n + 1) * g(X / p^n)，利用這個式子能夠大幅優化枚舉過程。

考慮到空間有限，咱們能夠將前N個數的約數個數g(1)~g(N)保存起來（個人代碼中取了N=2.4e8）。在質因數分解過程當中，若是X / p^n <= N，那麼就不用繼續試除，直接利用保存的值便可。

不過這樣優化以後耗時可能依然很長（N=3e8時在個人筆記本上耗時約35分鐘），因此咱們還須要一點「黑科技」——多線程。

多線程的思想就是把一個大任務分割成若干個小任務，而後讓它們同時進行。

對於本題，咱們能夠在預處理g(1)~g(N)以後，將剩餘的部分（N＋1 ~ 2e9）分割成若干塊，而後讓多個線程分別處理。

然而有個麻煩，就是用多線程處理的話，計算g(X)時咱們並不能知道g(1)~g(X-1)的最大值是多少。

對策：預處理時咱們能夠得出g(1)~g(N)的最大值M，而後在處理剩餘部分時，將全部知足g(X)>M的值X保存到數組vec，其他的值丟棄（絕大多數均可以被丟掉）。而後在vec中篩選反質數。

打表程序以下（clang++ -std=c++14 -O3編譯，在個人筆記本上預處理部分用了60s，多線程部分用了500s，加起來不到10分鐘）：

  1 #include <iostream>
  2 #include <cmath>
  3 #include <thread>
  4 #include <mutex>
  5 #include <atomic>
  6 #include <vector>
  7 #include <cstring>
  8 #include <algorithm>
  9 #include <boost/timer.hpp>
 10 
 11 std::vector<int> antiprimes;
 12 
 13 namespace prime
 14 {
 15     const int length = 100'000;
 16     bool primeFlag[length + 1];
 17     std::vector<int> primes;
 18 
 19     void init()
 20     {
 21         std::memset(primeFlag, 1, sizeof(primeFlag));
 22         for (int i = 2; i < 330; i++)
 23         {
 24             if (!primeFlag[i])
 25                 continue;
 26             for (int j = i << 1; j <= length; j += i)
 27                 primeFlag[j] = false;
 28         }
 29         for (int i = 2; i <= length; i++)
 30             if (primeFlag[i])
 31                 primes.push_back(i);
 32     }
 33 }
 34 
 35 namespace pre
 36 {
 37     const int length = 240'000'000;
 38     int factorN[length + 1];
 39     int maxFactorN = 0;
 40 }
 41 
 42 int getFactorN(int x)
 43 {
 44     int ans = 1;
 45     auto sqrtX = (int)std::sqrt(x);
 46     for (int p: prime::primes)
 47     {
 48         if (x == 1 || p > sqrtX)
 49             break;
 50         int cnt = 0;
 51         while (x % p == 0)
 52         {
 53             cnt += 1;
 54             x /= p;
 55         }
 56         if (cnt > 0)
 57         {
 58             ans *= (cnt + 1);
 59             if (x <= pre::length)
 60                 return ans * pre::factorN[x];
 61         }
 62     }
 63     if (x > 1)
 64         ans <<= 1;
 65     return ans;
 66 }
 67 
 68 namespace pre
 69 {
 70     void init()
 71     {
 72         maxFactorN = 0;
 73         for (int i = 1; i <= length; i++)
 74         {
 75             factorN[i] = getFactorN(i);
 76             if (factorN[i] > maxFactorN)
 77             {
 78                 maxFactorN = factorN[i];
 79                 antiprimes.push_back(i);
 80             }
 81         }
 82     }
 83 }
 84 
 85 struct VecNode
 86 {
 87     int val;
 88     int factorN;
 89     bool operator < (const VecNode& rhs) const {
 90         return val < rhs.val;
 91     }
 92 };
 93 std::vector<VecNode> selections;
 94 
 95 std::atomic_int taskN(0);
 96 
 97 //[first, last)
 98 void processEachThread(int first, int last, std::mutex& mtx)
 99 {
100     boost::timer timer;
101 
102     taskN += 1;
103     mtx.lock();
104     std::cout << "Thread #" << std::this_thread::get_id() << " is running.\n";
105     mtx.unlock();
106 
107     for (int cur = first; cur != last; ++cur)
108     {
109         int fn = getFactorN(cur);
110         if (fn > pre::maxFactorN)
111         {
112             mtx.lock();
113             selections.push_back({cur, fn});
114             mtx.unlock();
115         }
116     }
117     mtx.lock();
118     std::cout << "Thread #" << std::this_thread::get_id()
119               << " finished. first = " << first << ", last = " << last
120               << ", time used = " << timer.elapsed() << " sec.\n";
121     mtx.unlock();
122     taskN -= 1;
123 }
124 
125 void processThreads()
126 {
127     constexpr int maxTaskN = 6;
128     constexpr int taskRangeLength = 20'000'000;
129 
130     static_assert((2'000'000'000 - pre::length) % taskRangeLength == 0, "");
131 
132     std::vector<std::thread> threads;
133     std::mutex mtx;
134     int first(2'000'000'001 - taskRangeLength);
135     while (first >= pre::length)
136     {
137         while (taskN == maxTaskN) {}
138         threads.emplace_back(processEachThread, first, first + taskRangeLength, std::ref(mtx));
139         threads.back().detach();
140         first -= taskRangeLength;
141     }
142     while (taskN > 0) {}
143 }
144 
145 void pickSelections()
146 {
147     std::sort(selections.begin(), selections.end());
148     int maxFactorN = pre::maxFactorN;
149     for (auto node: selections)
150     {
151         if (node.factorN <= maxFactorN)
152             continue;
153         antiprimes.push_back(node.val);
154         maxFactorN = node.factorN;
155     }
156 }
157 
158 int main()
159 {
160     boost::timer timer;
161 
162     prime::init();
163     pre::init();
164     std::cout << "Initialization finished. Time used: " << timer.elapsed() << " sec.\n";
165 
166     timer.restart();
167     processThreads();
168     pickSelections();
169     std::cout << "Process finished. Time used: " << timer.elapsed() << " sec.\n";
170 
171     std::cout << "const int antiprimes[] = {";
172     for (int x: antiprimes)
173     {
174         if (x > 1)
175             std::cout << ", ";
176         std::cout << x;
177     }
178     std::cout << "};";
179 
180     return 0;
181 }