我正在解决这个任务（问题 I）。声明是：
该集合中有多少个子集{1, 2, 3, ..., n}是互质的？如果一组整数的每两个元素互质，则称为互质。如果两个整数的最大公约数等于 1，则它们互质。
输入
输入的第一行包含两个整数n和m( 1 <= n <= 3000, 1 <= m <= 10^9 + 9)
输出
{1, 2, 3, ..., n}输出模的互质子集的数量m。
例子
输入：4 7 输出：5
有 12 个互质子集{1,2,3,4}：{}, {1}, {2}, {3}, {4}, {1,2}, {1,3}, {1,4}, {2,3}, {3,4}, {1,2,3}, {1,3,4}.
我认为可以通过使用素数来解决。（跟踪我们是否使用了每个素数）..但我不确定。
我可以得到一些提示来解决这个任务吗？
好了，货到了。对我来说，接下来的 C 程序在不到 5 秒内得到 n=3000。我要向在竞争环境中解决这个问题的团队致敬。
该算法基于区别对待小素数和大素数的想法。如果素数的平方至多为 n，则素数是小的。否则，它很大。观察每个小于或等于 n 的数至多有一个大素因数。
我们制作了一个按对索引的表。每对的第一个组件指定使用的大素数的数量。每对的第二个组件指定使用的小素数集。特定索引处的值是该使用模式不包含 1 或大素数的解决方案的数量（我们稍后通过乘以 2 的适当幂来计算这些解决方案的数量）。
我们向下迭代数字 j，没有大的素因数。在每次迭代开始时，该表包含 j..n 子集的计数。内循环中有两个添加。第一个解释了 j 本身扩展子集，这不会增加使用的大素数的数量。第二个说明将子集扩展了 j 倍大素数，确实如此。合适的大素数的数量是不大于 n/j 的大素数的数量减去已经使用的大素数的数量，因为向下迭代意味着每个已经使用的大素数不大于 n/j。
最后，我们汇总表条目。表中计数的每个子集产生 2 ** k 个子集，其中 k 是未使用的大素数的数量加一，因为 1 并且每个未使用的大素数可以独立地包含或排除。
/* assumes int, long are 32, 64 bits respectively */
#include <stdio.h>
#include <stdlib.h>
enum {
NMAX = 3000
};
static int n;
static long m;
static unsigned smallfactors[NMAX + 1];
static int prime[NMAX - 1];
static int primecount;
static int smallprimecount;
static int largeprimefactor[NMAX + 1];
static int largeprimecount[NMAX + 1];
static long **table;
static void eratosthenes(void) {
int i;
for (i = 2; i * i <= n; i++) {
int j;
if (smallfactors[i]) continue;
for (j = i; j <= n; j += i) smallfactors[j] |= 1U << primecount;
prime[primecount++] = i;
}
smallprimecount = primecount;
for (; i <= n; i++) {
if (!smallfactors[i]) prime[primecount++] = i;
}
if (0) {
int k;
for (k = 0; k < primecount; k++) printf("%d\n", prime[k]);
}
}
static void makelargeprimefactor(void) {
int i;
for (i = smallprimecount; i < primecount; i++) {
int p = prime[i];
int j;
for (j = p; j <= n; j += p) largeprimefactor[j] = p;
}
}
static void makelargeprimecount(void) {
int i = 1;
int j;
for (j = primecount; j > smallprimecount; j--) {
for (; i <= n / prime[j - 1]; i++) {
largeprimecount[i] = j - smallprimecount;
}
}
if (0) {
for (i = 1; i <= n; i++) printf("%d %d\n", i, largeprimecount[i]);
}
}
static void maketable(void) {
int i;
int j;
table = calloc(smallprimecount + 1, sizeof *table);
for (i = 0; i <= smallprimecount; i++) {
table[i] = calloc(1U << smallprimecount, sizeof *table[i]);
}
table[0][0U] = 1L % m;
for (j = n; j >= 2; j--) {
int lpc = largeprimecount[j];
unsigned sf = smallfactors[j];
if (largeprimefactor[j]) continue;
for (i = 0; i < smallprimecount; i++) {
long *cur = table[i];
long *next = table[i + 1];
unsigned f;
for (f = sf; f < (1U << smallprimecount); f = (f + 1U) | sf) {
cur[f] = (cur[f] + cur[f & ~sf]) % m;
}
if (lpc - i <= 0) continue;
for (f = sf; f < (1U << smallprimecount); f = (f + 1U) | sf) {
next[f] = (next[f] + cur[f & ~sf] * (lpc - i)) % m;
}
}
}
}
static long timesexp2mod(long x, int y) {
long z = 2L % m;
for (; y > 0; y >>= 1) {
if (y & 1) x = (x * z) % m;
z = (z * z) % m;
}
return x;
}
static long computetotal(void) {
long total = 0L;
int i;
for (i = 0; i <= smallprimecount; i++) {
unsigned f;
for (f = 0U; f < (1U << smallprimecount); f++) {
total = (total + timesexp2mod(table[i][f], largeprimecount[1] - i + 1)) % m;
}
}
return total;
}
int main(void) {
scanf("%d%ld", &n, &m);
eratosthenes();
makelargeprimefactor();
makelargeprimecount();
maketable();
if (0) {
int i;
for (i = 0; i < 100; i++) {
printf("%d %ld\n", i, timesexp2mod(1L, i));
}
}
printf("%ld\n", computetotal());
return EXIT_SUCCESS;
}
这是一个答案，它在不到一秒的时间内遍历序列中的前 200 个元素，给出正确答案 200 → 374855124868136960。通过优化（参见编辑 1），它可以在 90 秒内计算前 500 个条目，这很快 - - 尽管@David Eisenstat 的答案可能会更好，如果它可以开发的话。我认为对迄今为止给出的算法采取了不同的方法，包括我自己的原始答案，所以我将它单独发布。
优化后，我意识到我实际上是在编写一个图形问题，所以我将解决方案重写为图形实现（参见编辑 2）。图形实现允许更多优化，更优雅，速度快一个数量级以上，并且扩展性更好：f(600)与 27 秒相比，它的计算时间为 1.5 秒。
这里的主要思想是使用递归关系。对于任何集合，满足标准的子集的数量是以下各项的总和：
在第二种情况下，当确定包含该元素时，必须删除任何其他不与它互质的元素。
效率问题：
代码如下。
#include <cassert>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <iostream>
#include <ctime>
const int PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
{ 2, 3, 5, ...
..., 2969, 2971, 2999 };
const int NPRIMES = sizeof(PRIMES) / sizeof(int);
typedef std::set<int> intset;
typedef std::vector<intset> intsetvec;
const int MAXCALC = 200; // answer at http://oeis.org/A084422/b084422.txt
intsetvec primeFactors(MAXCALC +1);
typedef std::vector<int> intvec;
// Caching / memoization
typedef std::map<intvec, double> intvec2dbl;
intvec2dbl set2NumCoPrimeSets;
double NumCoPrimeSets(const intvec& set)
{
if (set.empty())
return 1;
// Caching / memoization
const intvec2dbl::const_iterator i = set2NumCoPrimeSets.find(set);
if (i != set2NumCoPrimeSets.end())
return i->second;
// Result is the number of coprime sets in:
// setA, the set that definitely has the first element of the input present
// + setB, the set the doesn't have the first element of the input present
// Because setA definitely has the first element, we remove elements it isn't coprime with
// We also remove the first element: as this is definitely present it doesn't make any
// difference to the number of sets
intvec setA(set);
const int firstNum = *setA.begin();
const intset& factors = primeFactors[firstNum];
for(int factor : factors) {
setA.erase(std::remove_if(setA.begin(), setA.end(),
[factor] (int i) { return i % factor == 0; } ), setA.end());
}
// If the first element was already coprime with the rest, then we have setA = setB
// and we can do a single call (m=2). Otherwise we have two recursive calls.
double m = 1;
double c = 0;
assert(set.size() - setA.size() > 0);
if (set.size() - setA.size() > 1) {
intvec setB(set);
setB.erase(setB.begin());
c = NumCoPrimeSets(setB);
}
else {
// first elt coprime with rest
m = 2;
}
const double numCoPrimeSets = m * NumCoPrimeSets(setA) + c;
// Caching / memoization
set2NumCoPrimeSets.insert(intvec2dbl::value_type(set, numCoPrimeSets));
return numCoPrimeSets;
}
int main(int argc, char* argv[])
{
// Calculate prime numbers that factor into each number upto MAXCALC
primeFactors[1].insert(1); // convenient
for(int i=2; i<=MAXCALC; ++i) {
for(int j=0; j<NPRIMES; ++j) {
if (i % PRIMES[j] == 0) {
primeFactors[i].insert(PRIMES[j]);
}
}
}
const clock_t start = clock();
for(int n=1; n<=MAXCALC; ++n) {
intvec v;
for(int i=n; i>0; --i) { // reverse order to reduce recursion
v.push_back(i);
}
const clock_t now = clock();
const clock_t ms = now - start;
const double numCoPrimeSubsets = NumCoPrimeSets(v);
std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
}
return 0;
}
时间特征看起来比我的第一个答案好很多。但仍然不会在 5 秒内达到 3000！
编辑 1
可以对这种方法进行一些有趣的优化。总体而言，这为较大的n.
m数字被删除，则原始集合的组合比减少的组合多 2 m倍（因为对于每个互质，您可以将其放入或独立于其他元素的集合之外）。{2, 3, 15, 19, 45}，数字 15 和 45 具有相同的素因数 3 和 5。一次删除了2 个数字，因此 15 或 45 的子集数{2, 3, 15, 19, 45}=存在 15 或 45 的组合数的两倍（对于set{2, 19}因为如果存在 15 或 45，则 3 必须不存在）+ 15 和 45 不存在的子集数（对于 set {2, 3, 19}）short数字类型将性能提高了约 10%。{ 3, 9, 15}与 等价（同构）2, 4, 6。这是最激进的想法，但可能对性能的影响最小。用一个具体的例子可能更容易理解它。我选择了集合 {1..12}，它足够大，可以感受它的工作原理，但又足够小，以便于理解。
NumCoPrimeSets({ 1 2 3 4 5 6 7 8 9 10 11 12 })
Removed 3 coprimes, giving set { 2 3 4 5 6 8 9 10 12 } multiplication factor now 8
Removing the most connected number 12 with 8 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 3 }
setA = { 5 }
To get setB, remove 2 numbers which have *exactly* the prime factors { 2 3 }
setB = { 2 3 4 5 8 9 10 }
**** Recursing on 2 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)
NumCoPrimeSets({ 5 })
Base case return the multiplier, which is 2
NumCoPrimeSets({ 2 3 4 5 8 9 10 })
Removing the most connected number 10 with 4 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 5 }
setA = { 3 9 }
To get setB, remove 1 numbers which have *exactly* the prime factors { 2 5 }
setB = { 2 3 4 5 8 9 }
**** Recursing on 1 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)
NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Removing the most connected number 4 with 1 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 }
setA = { }
To get setB, remove 2 numbers which have *exactly* the prime factors { 2 }
setB = { }
**** Recursing on 2 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)
NumCoPrimeSets({ })
Base case return the multiplier, which is 1
NumCoPrimeSets({ })
Base case return the multiplier, which is 1
**** Returned from recursing on 2 * NumCoPrimeSets({ }) + NumCoPrimeSets({ })
Caching for{ 2 4 }: 3 = 2 * 1 + 1
Returning for{ 3 9 }: 3 = 1 * 3
NumCoPrimeSets({ 2 3 4 5 8 9 })
Removed 1 coprimes, giving set { 2 3 4 8 9 } multiplication factor now 2
Removing the most connected number 8 with 2 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 }
setA = { 3 9 }
To get setB, remove 3 numbers which have *exactly* the prime factors { 2 }
setB = { 3 9 }
**** Recursing on 3 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)
NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Cache hit, returning 3 = 1 * 3
NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Cache hit, returning 3 = 1 * 3
**** Returned from recursing on 3 * NumCoPrimeSets({ 3 9 }) + NumCoPrimeSets({ 3 9 })
Caching for{ 2 3 4 8 9 }: 12 = 3 * 3 + 3
Returning for{ 2 3 4 5 8 9 }: 24 = 2 * 12
**** Returned from recursing on 1 * NumCoPrimeSets({ 3 9 }) + NumCoPrimeSets({ 2 3 4 5 8 9 })
Caching for{ 2 3 4 5 8 9 10 }: 27 = 1 * 3 + 24
Returning for{ 2 3 4 5 8 9 10 }: 27 = 1 * 27
**** Returned from recursing on 2 * NumCoPrimeSets({ 5 }) + NumCoPrimeSets({ 2 3 4 5 8 9 10 })
Caching for{ 2 3 4 5 6 8 9 10 12 }: 31 = 2 * 2 + 27
Returning for{ 1 2 3 4 5 6 7 8 9 10 11 12 }: 248 = 8 * 31
下面的代码
#include <cassert>
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <ctime>
typedef short numtype;
const numtype PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
...
const numtype NPRIMES = sizeof(PRIMES) / sizeof(numtype);
typedef std::set<numtype> numset;
typedef std::vector<numset> numsetvec;
const numtype MAXCALC = 200; // answer at http://oeis.org/A084422/b084422.txt
numsetvec primeFactors(MAXCALC +1);
typedef std::vector<numtype> numvec;
// Caching / memoization
typedef std::map<numvec, double> numvec2dbl;
numvec2dbl set2NumCoPrimeSets;
double NumCoPrimeSets(const numvec& initialSet)
{
// Preprocessing step: remove numbers which are already coprime
typedef std::unordered_map<numtype, numvec> num2numvec;
num2numvec prime2Elts;
for(numtype num : initialSet) {
const numset& factors = primeFactors[num];
for(numtype factor : factors) {
prime2Elts[factor].push_back(num);
}
}
numset eltsToRemove(initialSet.begin(), initialSet.end());
typedef std::vector<std::pair<numtype,int>> numintvec;
numvec primesRemaining;
for(const num2numvec::value_type& primeElts : prime2Elts) {
if (primeElts.second.size() > 1) {
for (numtype num : primeElts.second) {
eltsToRemove.erase(num);
}
primesRemaining.push_back(primeElts.first);
}
}
double mult = pow(2.0, eltsToRemove.size());
if (eltsToRemove.size() == initialSet.size())
return mult;
// Do the removal by creating a new set
numvec set;
for(numtype num : initialSet) {
if (eltsToRemove.find(num) == eltsToRemove.end()) {
set.push_back(num);
}
}
// Transform to use a smaller set of primes before checking the cache
// (beta code but it seems to work, mostly!)
std::sort(primesRemaining.begin(), primesRemaining.end());
numvec::const_iterator p = primesRemaining.begin();
for(int j=0; p!= primesRemaining.end() && j<NPRIMES; ++p, ++j) {
const numtype primeRemaining = *p;
if (primeRemaining != PRIMES[j]) {
for(numtype& num : set) {
while (num % primeRemaining == 0) {
num = num / primeRemaining * PRIMES[j];
}
}
}
}
// Caching / memoization
const numvec2dbl::const_iterator i = set2NumCoPrimeSets.find(set);
if (i != set2NumCoPrimeSets.end())
return mult * i->second;
// Remove the most connected number
typedef std::unordered_map<numtype, int> num2int;
num2int num2ConnectionCount;
for(numvec::const_iterator srcIt=set.begin(); srcIt!=set.end(); ++srcIt) {
const numtype src = *srcIt;
const numset& srcFactors = primeFactors[src];
for(numvec::const_iterator tgtIt=srcIt +1; tgtIt!=set.end(); ++tgtIt) {
for(numtype factor : srcFactors) {
const numtype tgt = *tgtIt;
if (tgt % factor == 0) {
num2ConnectionCount[src]++;
num2ConnectionCount[tgt]++;
}
}
}
}
num2int::const_iterator connCountIt = num2ConnectionCount.begin();
numtype numToErase = connCountIt->first;
int maxConnCount = connCountIt->second;
for (; connCountIt!=num2ConnectionCount.end(); ++connCountIt) {
if (connCountIt->second > maxConnCount || connCountIt->second == maxConnCount && connCountIt->first > numToErase) {
numToErase = connCountIt->first;
maxConnCount = connCountIt->second;
}
}
// Result is the number of coprime sets in:
// setA, the set that definitely has a chosen element of the input present
// + setB, the set the doesn't have the chosen element(s) of the input present
// Because setA definitely has a chosen element, we remove elements it isn't coprime with
// We also remove the chosen element(s): as they are definitely present it doesn't make any
// difference to the number of sets
numvec setA(set);
const numset& factors = primeFactors[numToErase];
for(numtype factor : factors) {
setA.erase(std::remove_if(setA.begin(), setA.end(),
[factor] (numtype i) { return i % factor == 0; } ), setA.end());
}
// setB: remove all elements which have the same prime factors
numvec setB(set);
setB.erase(std::remove_if(setB.begin(), setB.end(),
[&factors] (numtype i) { return primeFactors[i] == factors; }), setB.end());
const size_t numEltsWithSamePrimeFactors = (set.size() - setB.size());
const double numCoPrimeSets =
numEltsWithSamePrimeFactors * NumCoPrimeSets(setA) + NumCoPrimeSets(setB);
// Caching / memoization
set2NumCoPrimeSets.insert(numvec2dbl::value_type(set, numCoPrimeSets));
return mult * numCoPrimeSets;
}
int main(int argc, char* argv[])
{
// Calculate prime numbers that factor into each number upto MAXCALC
for(numtype i=2; i<=MAXCALC; ++i) {
for(numtype j=0; j<NPRIMES; ++j) {
if (i % PRIMES[j] == 0) {
primeFactors[i].insert(PRIMES[j]);
}
}
}
const clock_t start = clock();
std::ofstream fout("out.txt");
for(numtype n=0; n<=MAXCALC; ++n) {
numvec v;
for(numtype i=1; i<=n; ++i) {
v.push_back(i);
}
const clock_t now = clock();
const clock_t ms = now - start;
const double numCoPrimeSubsets = NumCoPrimeSets(v);
fout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
}
return 0;
}
最多可n=600在 5 分钟内完成处理。然而，时间看起来仍然呈指数增长，每 50 到 60 倍n左右就翻一番。仅计算一个的图表n如下所示。
编辑 2
这个解决方案在图表方面更自然地实现。出现了另外两个优化：
最重要的是，如果图 G 可以划分为两个集合 A 和 B，使得 A 和 B 之间没有连接，那么 coprimes(G) = coprimes(A) * coprimes(B)。
其次，可以将一组素因子的所有数字折叠到一个节点中，因此该节点的值是其素因子的数字计数。
在下面的代码中，Graph该类保存了原始邻接矩阵和节点值，FilteredGraph该类保存了当前剩余节点的列表，bitset以便在删除节点时，可以通过位掩码计算新的邻接矩阵（并且有相对在递归中传递的数据很少）。
#include "Primes.h"
#include <cassert>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <iostream>
#include <ctime>
// Graph declaration
const int MAXGROUPS = 1462; // empirically determined
class Graph
{
typedef std::bitset<MAXGROUPS> bitset;
typedef std::vector<bitset> adjmatrix;
typedef std::vector<int> intvec;
public:
Graph(int numNodes)
: m_nodeValues(numNodes), m_adjMatrix(numNodes) {}
void SetNodeValue(int i, int v) { m_nodeValues[i] = v; }
void SetConnection(int i, int j)
{
m_adjMatrix[i][j] = true;
m_adjMatrix[j][i] = true;
}
int size() const { return m_nodeValues.size(); }
private:
adjmatrix m_adjMatrix;
intvec m_nodeValues;
friend class FilteredGraph;
};
class FilteredGraph
{
typedef Graph::bitset bitset;
public:
FilteredGraph(const Graph* unfiltered);
int FirstNode() const;
int RemoveNode(int node);
void RemoveNodesConnectedTo(int node);
double RemoveDisconnectedNodes();
bool AttemptPartition(FilteredGraph* FilteredGraph);
size_t Hash() const { return std::hash<bitset>()(m_includedNodes); }
bool operator==(const FilteredGraph& x) const
{ return x.m_includedNodes == m_includedNodes && x.m_unfiltered == m_unfiltered; }
private:
bitset RawAdjRow(int i) const {
return m_unfiltered->m_adjMatrix[i];
}
bitset AdjRow(int i) const {
return RawAdjRow(i) & m_includedNodes;
}
int NodeValue(int i) const {
return m_unfiltered->m_nodeValues[i];
}
const Graph* m_unfiltered;
bitset m_includedNodes;
};
// Cache
namespace std {
template<>
class hash<FilteredGraph> {
public:
size_t operator()(const FilteredGraph & x) const { return x.Hash(); }
};
}
typedef std::unordered_map<FilteredGraph, double> graph2double;
graph2double cache;
// MAIN FUNCTION
double NumCoPrimesSubSets(const FilteredGraph& graph)
{
graph2double::const_iterator cacheIt = cache.find(graph);
if (cacheIt != cache.end())
return cacheIt->second;
double rc = 1;
FilteredGraph A(graph);
FilteredGraph B(graph);
if (A.AttemptPartition(&B)) {
rc = NumCoPrimesSubSets(A);
A = B;
}
const int nodeToRemove = A.FirstNode();
if (nodeToRemove < 0) // empty graph
return 1;
// Graph B is the graph with a node removed
B.RemoveNode(nodeToRemove);
// Graph A is the graph with the node present -- and hence connected nodes removed
A.RemoveNodesConnectedTo(nodeToRemove);
// The number of numbers in the node is the number of times it can be reused
const double removedNodeValue = A.RemoveNode(nodeToRemove);
const double A_disconnectedNodesMult = A.RemoveDisconnectedNodes();
const double B_disconnectedNodesMult = B.RemoveDisconnectedNodes();
const double A_coprimes = NumCoPrimesSubSets(A);
const double B_coprimes = NumCoPrimesSubSets(B);
rc *= removedNodeValue * A_disconnectedNodesMult * A_coprimes
+ B_disconnectedNodesMult * B_coprimes;
cache.insert(graph2double::value_type(graph, rc));
return rc;
}
// Program entry point
int Sequence2Graph(Graph** ppGraph, int n);
int main(int argc, char* argv[])
{
const clock_t start = clock();
int n=800; // runs in approx 6s on my machine
Graph* pGraph = nullptr;
const int coPrimesRemoved = Sequence2Graph(&pGraph, n);
const double coPrimesMultiplier = pow(2,coPrimesRemoved);
const FilteredGraph filteredGraph(pGraph);
const double numCoPrimeSubsets = coPrimesMultiplier * NumCoPrimesSubSets(filteredGraph);
delete pGraph;
cache.clear(); // as it stands the cache can't cope with other Graph objects, e.g. for other n
const clock_t now = clock();
const clock_t ms = now - start;
std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
return 0;
}
// Graph implementation
FilteredGraph::FilteredGraph(const Graph* unfiltered)
: m_unfiltered(unfiltered)
{
for(int i=0; i<m_unfiltered->size(); ++i) {
m_includedNodes.set(i);
}
}
int FilteredGraph::FirstNode() const
{
int firstNode=0;
for(; firstNode<m_unfiltered->size() && !m_includedNodes.test(firstNode); ++firstNode) {
}
if (firstNode == m_unfiltered->size())
return -1;
return firstNode;
}
int FilteredGraph::RemoveNode(int node)
{
m_includedNodes.set(node, false);
return NodeValue(node);
}
void FilteredGraph::RemoveNodesConnectedTo(const int node)
{
const bitset notConnected = ~RawAdjRow(node);
m_includedNodes &= notConnected;
}
double FilteredGraph::RemoveDisconnectedNodes()
{
double mult = 1.0;
for(int i=0; i<m_unfiltered->size(); ++i) {
if (m_includedNodes.test(i)) {
const int conn = AdjRow(i).count();
if (conn == 0) {
m_includedNodes.set(i, false);;
mult *= (NodeValue(i) +1);
}
}
}
return mult;
}
bool FilteredGraph::AttemptPartition(FilteredGraph* pOther)
{
typedef std::vector<int> intvec;
intvec includedNodesCache;
includedNodesCache.reserve(m_unfiltered->size());
for(int i=0; i<m_unfiltered->size(); ++i) {
if (m_includedNodes.test(i)) {
includedNodesCache.push_back(i);
}
}
if (includedNodesCache.empty())
return false;
const int startNode= includedNodesCache[0];
bitset currFoundNodes;
currFoundNodes.set(startNode);
bitset foundNodes;
do {
foundNodes |= currFoundNodes;
bitset newFoundNodes;
for(int i : includedNodesCache) {
if (currFoundNodes.test(i)) {
newFoundNodes |= AdjRow(i);
}
}
newFoundNodes &= ~ foundNodes;
currFoundNodes = newFoundNodes;
} while(currFoundNodes.count() > 0);
const size_t foundCount = foundNodes.count();
const size_t thisCount = m_includedNodes.count();
const bool isConnected = foundCount == thisCount;
if (!isConnected) {
if (foundCount < thisCount) {
pOther->m_includedNodes = foundNodes;
m_includedNodes &= ~foundNodes;
}
else {
pOther->m_includedNodes = m_includedNodes;
pOther->m_includedNodes &= ~foundNodes;
m_includedNodes = foundNodes;
}
}
return !isConnected;
}
// Initialization code to convert sequence from 1 to n into graph
typedef short numtype;
typedef std::set<numtype> numset;
bool setIntersect(const numset& setA, const numset& setB)
{
for(int a : setA) {
if (setB.find(a) != setB.end())
return true;
}
return false;
}
int Sequence2Graph(Graph** ppGraph, int n)
{
typedef std::map<numset, int> numset2int;
numset2int factors2count;
int coPrimesRemoved = n>0; // for {1}
// Calculate all sets of prime factors, and how many numbers belong to each set
for(numtype i=2; i<=n; ++i) {
if ((i > n/2) && (std::find(PRIMES, PRIMES+NPRIMES, i) !=PRIMES+NPRIMES)) {
++coPrimesRemoved;
}
else {
numset factors;
for(numtype j=0; j<NPRIMES && PRIMES[j]<n; ++j) {
if (i % PRIMES[j] == 0) {
factors.insert(PRIMES[j]);
}
}
factors2count[factors]++;
}
}
// Create graph
Graph*& pGraph = *ppGraph;
pGraph = new Graph(factors2count.size());
int srcNodeNum = 0;
for(numset2int::const_iterator i = factors2count.begin(); i!=factors2count.end(); ++i) {
pGraph->SetNodeValue(srcNodeNum, i->second);
numset2int::const_iterator j = i;
int tgtNodeNum = srcNodeNum+1;
for(++j; j!=factors2count.end(); ++j) {
if (setIntersect(i->first, j->first)) {
pGraph->SetConnection(srcNodeNum, tgtNodeNum);
}
++tgtNodeNum;
}
++srcNodeNum;
}
return coPrimesRemoved;
}
计算互质数( n) 的图表如下所示为红色（旧方法为黑色）。
根据观察到的（指数）增长率n=3000，假设程序没有崩溃，预测为 30 小时。这开始看起来在计算上是可行的，尤其是在进行更多优化的情况下，但距离所需的 5s 还差得很远！毫无疑问，所需的解决方案是简短而甜蜜的，但这很有趣......
这是 Haskell 中相当简单的事情，对于 n=200 大约需要 2 秒，并且呈指数增长。
{-# OPTIONS_GHC -O2 #-}
f n = 2^(length second + 1) * (g [] first 0) where
second = filter (\x -> isPrime x && x > div n 2) [2..n]
first = filter (flip notElem second) [2..n]
isPrime k =
null [ x | x <- [2..floor . sqrt . fromIntegral $ k], k `mod`x == 0]
g s rrs depth
| null rrs = 2^(length s - depth)
| not $ and (map ((==1) . gcd r) s) = g s rs depth
+ g s' rs' (depth + 1)
| otherwise = g (r:s) rs depth
where r:rs = rrs
s' = r : filter ((==1) . gcd r) s
rs' = filter ((==1) . gcd r) rs
这是一种方法，可以在 5 秒内获得给定的序列（经过n=62优化，它n=75在 5 秒内运行，但请注意我在这个问题上的第二次尝试效果更好）。我假设问题的模部分只是为了避免随着函数变大而出现数值错误，所以我现在忽略它。
该方法基于这样一个事实，即我们最多可以在一个子集中为每个素数选择一个数字。
举个例子{1,2,3,4}，我们将集合中的数字映射到素数上，如下所示。我已将 1 作为质数包括在内，因为它使现阶段的说明更容易。
1 → {1}
2 → {2,4}
3 → {3}
“素数” 1 有 2 种组合（不包括它或 1），素数 2 有 3 种组合（不包括它或 2 或 4），3 有 2 种组合（不包括它或3）。所以子集的个数是2 * 3 * 2 = 12。
同样，{1,2,3,4,5}我们有
1 → {1}
2 → {2,4}
3 → {3}
5 → {5}
给予2 * 3 * 2 * 2= 24.
但是对于{1,2,3,4,5,6}，事情并不是那么简单。我们有
1 → {1}
2 → {2,4,6}
3 → {3}
5 → {5}
但是如果我们为素数 2 选择数字 6，我们就不能为素数 3 选择一个数字（作为脚注，在我可能会回到的第一种方法中，我将其视为 3 的选择是当我们选择 6 时减少了一半，所以我使用 3.5 而不是 4 作为素数 2 的组合数量并2 * 3.5 * 2 * 2 = 28给出了正确的答案。但是，我无法让这种方法超越n=17。）
我处理这个问题的方法是在每个级别对每组主要因素进行处理。所以{2,4}有素因数{2}，而{6}有素因数{2,3}。省略 1 的虚假条目（它不是素数），我们现在有
2 → {{2}→{2,4}, {2,3}→6}
3 → {{3}→{3}}
5 → {{5}→{5}}
现在有三种路径来计算组合的数量，一条我们不选择素数 2 的路径，还有两条我们选择的路径： through{2}→{2,4}和 through {2,3}→6。
1 * 2 * 2 = 4组合，因为我们可以选择 3 或不选择，也可以选择 5 或不选择。2 * 2 * 2 = 8组合，注意我们可以选择 2 或 4。1 * 1 * 2 = 2组合，因为我们对素数 3 只有一个选择——不要使用它。总的来说，这给了我们4 + 8 + 2 = 14组合（作为优化说明，第一条和第二条路径可以折叠在一起得到3 * 2 * 2 = 12）。我们也可以选择是否选择 1，所以组合的总数是2 * 14 = 28。
下面是递归运行路径的 C++ 代码。（它是 C++11，在 Visual Studio 2012 上编写，但应该可以在其他 gcc 上工作，因为我没有包含任何特定于平台的内容）。
#include <cassert>
#include <vector>
#include <set>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <ctime>
const int PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
{ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
103, 107, 109, 113, 127, 131, 137, 139, 149,
151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199 };
const int NPRIMES = sizeof(PRIMES) / sizeof(int);
typedef std::vector<int> intvec;
typedef std::set<int> intset;
typedef std::vector<std::set<int>> intsetvec;
struct FactorSetNumbers
{
intset factorSet;
intvec numbers; // we only need to store numbers.size(), but nice to see the vec itself
FactorSetNumbers() {}
FactorSetNumbers(const intset& factorSet_, int n)
: factorSet(factorSet_)
{
numbers.push_back(n);
}
};
typedef std::vector<FactorSetNumbers> factorset2numbers;
typedef std::vector<factorset2numbers> factorset2numbersArray;
double NumCoPrimeSubsets(
const factorset2numbersArray& factorSet2Numbers4FirstPrime,
int primeIndex, const intset& excludedPrimes)
{
const factorset2numbers& factorSet2Numbers = factorSet2Numbers4FirstPrime[primeIndex];
if (factorSet2Numbers.empty())
return 1;
// Firstly, we may choose not to use this prime number at all
double numCoPrimeSubSets = NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
primeIndex + 1, excludedPrimes);
// Optimization: if we're not excluding anything, then we can collapse
// the above call and the first call in the loop below together
factorset2numbers::const_iterator i = factorSet2Numbers.begin();
if (excludedPrimes.empty()) {
const FactorSetNumbers& factorSetNumbers = *i;
assert(factorSetNumbers.factorSet.size() == 1);
numCoPrimeSubSets *= (1 + factorSetNumbers.numbers.size());
++i;
}
// We are using this prime number. The number of subsets for this prime number is the sum of
// the number of subsets for each set of integers whose factors don't include an excluded factor
for(; i!=factorSet2Numbers.end(); ++i) {
const FactorSetNumbers& factorSetNumbers = *i;
intset intersect;
std::set_intersection(excludedPrimes.begin(),excludedPrimes.end(),
factorSetNumbers.factorSet.begin(),factorSetNumbers.factorSet.end(),
std::inserter(intersect,intersect.begin()));
if (intersect.empty()) {
intset unionExcludedPrimes;
std::set_union(excludedPrimes.begin(),excludedPrimes.end(),
factorSetNumbers.factorSet.begin(),factorSetNumbers.factorSet.end(),
std::inserter(unionExcludedPrimes,unionExcludedPrimes.begin()));
// Optimization: don't exclude on current first prime,
// because can't possibly occur later on
unionExcludedPrimes.erase(unionExcludedPrimes.begin());
numCoPrimeSubSets += factorSetNumbers.numbers.size() *
NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
primeIndex + 1, unionExcludedPrimes);
}
}
return numCoPrimeSubSets;
}
int main(int argc, char* argv[])
{
const int MAXCALC = 80;
intsetvec primeFactors(MAXCALC +1);
// Calculate prime numbers that factor into each number upto MAXCALC
for(int i=2; i<=MAXCALC; ++i) {
for(int j=0; j<NPRIMES; ++j) {
if (i % PRIMES[j] == 0) {
primeFactors[i].insert(PRIMES[j]);
}
}
}
const clock_t start = clock();
factorset2numbersArray factorSet2Numbers4FirstPrime(NPRIMES);
for(int n=2; n<=MAXCALC; ++n) {
{
// For each prime, store all the numbers whose first prime factor is that prime
// E.g. for the prime 2, for n<=20, we store
// {2}, { 2, 4, 8, 16 }
// {2, 3}, { 6, 12, 18 }
// {2, 5}, { 5, 10, 20 }
// {2, 7}, { 14 }
const int firstPrime = *primeFactors[n].begin();
const int firstPrimeIndex = std::find(PRIMES, PRIMES + NPRIMES, firstPrime) - PRIMES;
factorset2numbers& factorSet2Numbers = factorSet2Numbers4FirstPrime[firstPrimeIndex];
const factorset2numbers::iterator findFactorSet = std::find_if(factorSet2Numbers.begin(), factorSet2Numbers.end(),
[&](const FactorSetNumbers& x) { return x.factorSet == primeFactors[n]; });
if (findFactorSet == factorSet2Numbers.end()) {
factorSet2Numbers.push_back(FactorSetNumbers(primeFactors[n], n));
}
else {
findFactorSet->numbers.push_back(n);
}
// The number of coprime subsets is the number of coprime subsets for the first prime number,
// starting with an empty exclusion list
const double numCoPrimeSubSetsForNEquals1 = 2;
const double numCoPrimeSubsets = numCoPrimeSubSetsForNEquals1 *
NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
0, // primeIndex
intset()); // excludedPrimes
const clock_t now = clock();
const clock_t ms = now - start;
std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
}
}
return 0;
}
Timings：在 <0.1s 内计算最多 40 个序列，在 0.5s 内计算最多 50 个序列，在 2.5s 内计算到 60 个，在 20s 内计算到 70 个，在 157s 内计算到 80 个。
尽管这似乎输出了正确的数字，但正如预期的那样，成本太高了。特别是它至少需要指数时间（很可能是组合时间）。
显然，这种方法不能按要求扩展。但是，这里可能有一些东西可以给其他人一些想法（或将这种方法排除在失败之外）。好像有两种可能：
bitset而不是set.n=17和以上的情况下非常有效地工作，并且失败了n=18，被少数人淘汰。我花了很长时间编写模式并试图弄清楚为什么它突然失败了n=18，但看不到它。我可能会回到这个，或者如果有人感兴趣，我会将其作为替代答案。编辑：我使用一些技巧进行了一些优化，尽量避免在可能的情况下重做现有计算，并且代码速度提高了大约 10 倍。听起来不错，但这只是不断的改进。真正需要的是对这个问题的一些了解——例如我们可以#subsets(n+1)基于#subsets(n)吗？
我会这样做：
mod m找出最多数的质因数nq，并将空集合添加到其中，并将计数器设置为 1X从队列中弹出一个元素k从max(X)到的所有数字n，检查集合的因子是否与数字的因子相交。如果不是，则添加到队列中X U k并将计数器加 1。否则，转到下一个 k。 必须指出两个重要的事情：
n，而只需要它们的素因数，这意味着，对于 12，您只需要 2 和 3。这样检查 2 个数字是否互质变成了检查两组的交集是否为空。这是 O(n*2^p) 中的一种方式，其中p是 下的素数数n。不使用模数。
class FailureCoprimeSubsetCounter{
int[] primes;//list of primes under n
PrimeSet[] primeSets;//all 2^primes.length
//A set of primes under n. And a count which goes with it.
class PrimeSet{
BitSet id;//flag x is 1 iff prime[x] is a member of this PrimeSet
long tally;//number of coprime sets that do not have a factor among these primes and do among all the other primes
//that is, we count the number of coprime sets whose maximal coprime subset of primes[] is described by this object
PrimeSet(int np){...}
}
int coprimeSubsets(int n){
//... initialization ...
for(int k=1; k<=n; k++){
PrimeSet p = listToPrimeSet(PrimeFactorizer.factorize(k));
for(int i=0; i<Math.pow(2,primes.length); i++){
//if p AND primes[i] is empty
//add primes[i].tally to PrimeSet[ p OR primes[i] ]
}
}
//return sum of all the tallies
}
}
但既然这是一个竞争问题，就必须有一个更快更肮脏的解决方案。而且由于这种方法需要指数时间和空间，并且有 430 个低于 3000 的素数，所以也更优雅。
编辑：添加了递归方法。在 5 秒内求解，直到 n = 50。
#include <iostream>
#include <vector>
using namespace std;
int coPrime[3001][3001] = {0};
int n, m;
// function that checks whether a new integer is coprime with all
//elements in the set S.
bool areCoprime ( int p, vector<int>& v ) {
for ( int i = 0; i < v.size(); i++ ) {
if ( !coPrime[v[i]][p] )
return false;
}
return true;
}
// implementation of Euclid's GCD between a and b
bool isCoprimeNumbers( int a, int b ) {
for ( ; ; ) {
if (!(a %= b)) return b == 1 ;
if (!(b %= a)) return a == 1 ;
}
}
int subsets( vector<int>& coprimeList, int index ) {
int count = 0;
for ( int i = index+1; i <= n; i++ ) {
if ( areCoprime( i, coprimeList ) ) {
count = ( count + 1 ) % m;
vector<int> newVec( coprimeList );
newVec.push_back( i );
count = ( count + subsets( newVec, i ) ) % m;
}
}
return count;
}
int main() {
cin >> n >> m;
int count = 1; // empty set
count += n; // sets with 1 element each.
// build coPrime matrix
for ( int i = 1; i <= 3000; i++ )
for ( int j = i+1; j <= 3000; j++ )
if ( isCoprimeNumbers( i, j ) )
coPrime[i][j] = 1;
// find sets beginning with i
for ( int i = 1; i <= n; i++ ) {
vector<int> empty;
empty.push_back( i );
count = ( count + subsets( empty, i ) ) % m;
}
cout << count << endl;
return 0;
}
一种天真的方法可以是（对于 N = 3000）：
第一步：建立布尔矩阵
1. 建立一个从 2 到 1500 的素数列表。
2. 对于 1 到 3000 的每个数字，建立其素数因子的集合。
3. 比较每对集合，得到一个布尔矩阵[3000][3000]，它表明元素 i 和 j 是否互质 (1) 或不互质 (0)。
步骤 2：计算长度为 k 的互质集的数量（k = 0 到 3000）
1. 初始化 count = 1（空集）。现在 k = 1。维护一个长度为 k 的集合列表。
2. 构建仅包含该特定元素的 3000 个集合。（增加计数） 3. 从k 到 3000
扫描每个元素，看看是否可以通过将其添加到任何现有的长度为 k 的集合中来形成一个新集合。注意：一些新形成的集合可能相同也可能不同。如果您使用集合集，则仅应存储唯一集合。
4.删除所有长度仍为 k 的集合。
5. 按当前唯一集数增加计数。
6. k = k + 1 并转到步骤 3。
或者，您可以维护集合中每个元素的产品列表，并检查新元素是否与产品互质。在这种情况下，您不需要存储布尔矩阵。
上述算法的复杂度似乎介于 O(n^2) 和 O(n^3) 之间。
C++ 中的完整代码：（改进：添加的条件是，仅当元素大于集合中的最大值时，才应在集合中检查元素）。
#include <iostream>
#include <vector>
#include <set>
using namespace std;
int coPrime[3001][3001] = {0};
// function that checks whether a new integer is coprime with all
//elements in the set S.
bool areCoprime ( int p, set<int> S ) {
set<int>::iterator it_set;
for ( it_set = S.begin(); it_set != S.end(); it_set++ ) {
if ( !coPrime[p][*it_set] )
return false;
}
return true;
}
// implementation of Euclid's GCD between a and b
bool isCoprimeNumbers( int a, int b ) {
for ( ; ; ) {
if (!(a %= b)) return b == 1 ;
if (!(b %= a)) return a == 1 ;
}
}
int main() {
int n, m;
cin >> n >> m;
int count = 1; // empty set
set< set<int> > setOfSets;
set< set<int> >::iterator it_setOfSets;
// build coPrime matrix
for ( int i = 1; i <= 3000; i++ )
for ( int j = 1; j <= 3000; j++ )
if ( i != j && isCoprimeNumbers( i, j ) )
coPrime[i][j] = 1;
// build set of sets containing 1 element.
for ( int i = 1; i <= n; i++ ) {
set<int> newSet;
newSet.insert( i );
setOfSets.insert( newSet );
count = (count + 1) % m;
}
// Make sets of length k
for ( int k = 2; k <= n; k++ ) {
// Scane each element from k to n
set< set<int> > newSetOfSets;
for ( int i = k; i <= n; i++ ) {
//Scan each existing set.
it_setOfSets = setOfSets.begin();
for ( ; it_setOfSets != setOfSets.end(); it_setOfSets++ ) {
if ( i > *(( *it_setOfSets ).rbegin()) && areCoprime( i, *it_setOfSets ) ) {
set<int> newSet( *it_setOfSets );
newSet.insert( i );
newSetOfSets.insert( newSet );
}
}
}
count = ( count + newSetOfSets.size() ) % m;
setOfSets = newSetOfSets;
}
cout << count << endl;
return 0;
}
上面的代码似乎给出了正确的结果，但消耗了大量时间：假设 M 足够大：
For N = 4, count = 12. (almost instantaneous)
For N = 20, count = 3232. (2-3 seconds)
For N = 25, count = 11168. (2-3 seconds)
For N = 30, count = 31232 (4 seconds)
For N = 40, count = 214272 (30 seconds)
这是我之前提到的不同方法。
它确实比我以前使用的快得多。使用在线解释器（ideone），它可以coprime_subsets(117)在不到 5 秒的时间内完成计算。
该代码从空集开始构建可能的集合，并将每个数字插入所有可能的子集中。
primes_to_3000 = set([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999])
# primes up to sqrt(3000), used for factoring numbers
primes = set([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53])
factors = [set() for _ in xrange(3001)]
for p in primes:
for n in xrange(p, 3001, p):
factors[n].add(p)
def coprime_subsets(highest):
count = 1
used = {frozenset(): 1}
for n in xrange(1, highest+1):
if n in primes_to_3000:
# insert the primes into all sets
count <<= 1
if n < 54:
used.update({k.union({n}): v for k, v in used.iteritems()})
else:
for k in used:
used[k] *= 2
else:
for k in used:
# only insert into subsets that don't share any prime factors
if not factors[n].intersection(k):
count += used[k]
used[k.union(factors[n])] += used[k]
return count
这是我的想法和python中的实现。它似乎很慢，但我不确定这是否只是我测试的方式（使用在线解释器）......
可能是在“真实”计算机上运行它可能会有所作为，但我可以'现在测试一下。
这种方法有两个部分：
在那之后，我猜你只是取模......
这是我的python实现（改进版）：
# primes up to 1500
primes = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499
factors = [set() for _ in xrange(3001)]
for p in primes:
for n in xrange(p, 3001, p):
factors[n].add(p)
def coprime_subsets(highest, current=1, factors_used=frozenset(), cache={}):
"""
Determine the number of possible coprime subsets of numbers,
using numbers starting at index current.
factor_product is used for determining if a number can be added
to the current subset.
"""
if (current, factors_used) in cache:
return cache[current, factors_used]
count = 1
for n in xrange(current, highest+1):
if factors_used.intersection(factors[n]):
continue
count += coprime_subsets(highest, n+1, factors_used.union(factors[n]))
cache[current, factors_used] = count
return count
我还有另一个想法，如果有时间我会尝试实施。我相信另一种方法可能会更快一些。
看起来提议的答案以及问题的序言解决了一个与所提出的问题不同的问题。
问题是：
Output the number of coprime subsets of {1, 2, 3, ..., n} modulo m.
建议的答案试图解决另一个问题：
Output the number of coprime subsets of {1, 2, 3, ..., n}.
这些问题并不等同。第一个处理有限环 Z_m，第二个处理具有完全不同算术的整数环 Z。
此外，问题前言中的定义“如果最大公约数等于 1，则两个整数互质”不适用于 Z_m：有限环没有算术上稳定的比较，因此不存在“最大”公共除数。
同样的反对意见也适用于问题中的示例：3 和 4 不是模 7 的互质数，因为两者都可以被模 7 的 2 整除：4=(2*2)%7 和 3=(2*5)%7。
事实上，Z_m 算术非常奇怪，以至于可以在 O(1) 时间内给出答案，至少对于素数 m 来说是这样：对于任何 n 和素数 m，没有模 m 的互素对。原因如下：Z_m 的非零元素形成一个 m-1 阶循环群，这意味着对于 Z_m 中的任何非零元素 a 和 b，对于 Z_m 中的某个 c，可以写 a=bc。这证明在 Z_m 中对于素数 m 没有互素对。
从问题中的示例：让我们看一下 {2, 3} mod 7 和 {3, 4} mod 7：2=(3*3)%7 和 3=(4*6)%7。因此 {2,3} 在 Z_7 中不是互质的（都可以被 3 整除）并且 {3,4} 在 Z_7 中不是互质的（都可以被 4 整除）。
现在让我们考虑非素数 m 的情况。将 ma 写为素数幂 m=p_1^i_1*...*p_k^i_k 的乘积。如果 a 和 b 有一个共同的质因子，那么它们显然不是互质的。如果其中至少一个，比如 b，不整除任何素数 p_1,...,p_k，那么 a 和 b 有一个公因数，原因与素数 m 的情况大致相同：b 将是一个乘法Z_m 的单位，因此在 Z_m 中 a 可以被 b 整除。
因此，仍然需要考虑当 m 是合数并且 a 和 b 可以被 m 的不同素因数整除的情况，假设 a 可以被 p 整除，b 可以被 q 整除。在这种情况下，它们确实可以互质。例如，2 和 3 模 6 是互质数。
因此，互质对的问题归结为以下步骤：
找到小于 n 的 m 的质因子。如果没有，则没有互质对。
枚举这些素因子的幂乘积，这些素因子仍然是小于 n 的 m 的因子。
计算这些之间的 Z-comprime 对的数量。