luogu 1463

题目描述:

对于任何正整数x,其约数的个数记作g(x)。例如g(1)=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.
";
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.
";
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.
";
165 
166     timer.restart();
167     processThreads();
168     pickSelections();
169     std::cout << "Process finished. Time used: " << timer.elapsed() << " sec.
";
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 }
原文地址:https://www.cnblogs.com/Onlynagesha/p/8584079.html