如何计算集合 1,2,3,..,n 的互质子集的数量
Posted
技术标签:
【中文标题】如何计算集合 1,2,3,..,n 的互质子集的数量【英文标题】：How to calculate the number of coprime subsets of the set 1,2,3,..,n如何计算集合 1,2,3,..,n 的互质子集的数量 【发布时间】：2013-09-18 09:00:56 【问题描述】：我正在解决this task (problem I)。声明是：
集合1, 2, 3, ..., n 有多少个子集是互质的？如果一组整数的每两个元素互质，则称为互质。如果两个整数的最大公约数等于 1，则它们互质。
输入
输入的第一行包含两个整数n 和m (1 <= n <= 3000, 1 <= m <= 10^9 + 9)
输出
输出1, 2, 3, ..., n模m的互质子集的数量。
例子
输入：4 7 输出：5
1,2,3,4 有 12 个互质子集：、1、2、3、4、1,2、1,3、1,4、@986 987654344@, 1,2,3, 1,3,4.
我认为可以通过使用素数来解决。 （跟踪我们是否使用了每个素数）..但我不确定。
我可以得到一些提示来解决这个任务吗？
你可以在这里找到这个序列：http://oeis.org/A084422【问题讨论】：
如果这个问题太难了，试着解决一个更简单的问题。有多少对小于或等于 n 的数互质？或者更简单：有多少小于 n 的数与 n 互质？ 互质要求立即让我想到了Euler totient。也许this 论文会有所帮助。 这里是相关论文。 math.clemson.edu/~calkin/Papers/calkin_granville.pdf 我相信定理 2 就是你要找的，祝你好运。 @octref，我刚从OEIS(A084422找到这个，而this是序列的表。比如@时可以看到有374855124868136960互质子集987654348@.
@BrettHale 好吧，如果你认为如果不存在 gcd > 1 的两个不同元素的集合是互质的，那么单例和空集都可以。
【参考方案1】：
好的，这里有货。下面的 C 程序在 less 中得到 n=3000 对我来说超过 5 秒。我要感谢解决这个问题的团队 竞争环境中的问题。
该算法是基于处理大小的思想 素数不同。如果素数的平方至多为 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;
【讨论】：
我不明白你的第 4 步。为什么 7, 9, 10 变成 9,2？ @justhalf 将大素数 7 和 5 分开后变成 1, 9, 2，然后是 9, 2。 进展顺利 -- 比起团队合作，您自己完成这项工作给我留下了深刻的印象。 太棒了！我能理解你的解释。感谢您的回复：） 这似乎是一个旧消息 - 尽管这个问题在一个月的大部分时间里一直让我好奇 - 但你能帮助我了解你的程序如何（有效地）计算“解决方案的数量具有特定的使用模式？”我无法理解它。无论如何，阅读它教会了我新的东西。非常聪明。【参考方案2】：这是一个答案，它在不到一秒的时间内遍历sequence 中的前 200 个元素，给出正确答案 200 → 374855124868136960。通过优化（参见编辑 1），它可以计算 90 秒以下的前 500 个条目，这很快 - 尽管@David Eisenstat 的答案可能会更好，如果它可以开发。我认为对目前给出的算法采取了不同的方法，包括我自己的原始答案，所以我单独发布。
优化后，我意识到我实际上是在编写一个图形问题，因此我将解决方案重写为图形实现（参见编辑 2）。图形实现允许进行更多优化，更加优雅，速度快了一个数量级以上，并且扩展性更好：与 27 秒相比，它在 1.5 秒内计算 f(600)。
这里的主要思想是使用递归关系。对于任何集合，满足条件的子集的数量是：
删除一个元素的子集数；和 绝对包含该元素的子集数。在第二种情况下，当元素被确定包含时，任何其他不与它互质的元素都必须被删除。
效率问题：
我选择删除最大的元素，以最大限度地提高该元素与所有其他元素互质的机会，在这种情况下，只需进行一次而不是两次递归调用。 缓存/记忆有助于。代码如下。
#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;
时间特征看起来比my first answer好很多。但是5s内还是不会涨到3000！
编辑 1
可以对此方法进行一些有趣的优化。总体而言，这为较大的n 提供了 4 倍的改进。
m 数字，则原始集合的组合比减少的组合多 2m 倍一个（因为对于每个互质数，您可以将其独立于其他元素进入或退出集合）。
最重要的是，可以选择集合中任意位置的元素来移除。事实证明，删除连接最多的元素效果最好。
以前使用的递归关系可以推广到删除多个元素，其中所有删除的元素都具有相同的素因子。例如。对于集合2, 3, 15, 19, 45，数字 15 和 45 具有相同的质因数 3 和 5。一次删除 2 个数字，因此 2, 3, 15, 19, 45 的子集数 = 两倍存在 15 或 45 的组合数（对于集合 2, 19，因为如果存在 15 或 45，则必须不存在 3）+不存在 15 和 45 的子集数（对于设置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;
最多可以在 5 分钟内处理到 n=600。然而，时间看起来仍然呈指数级增长，每 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 相去甚远！毫无疑问，所需的解决方案是简短而甜蜜的，但这很有趣......
【讨论】：
【参考方案3】：这是 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
【讨论】：
这让我想起了我的第一次尝试。我怀疑对于大的 n 会有太多的素数分区。 @DavidEisenstat 感谢您的检查。haskell: length $ factorsets 3000 => 1823 （一个素数的不同幂算作一个因子集）这不意味着将最多长度为 431 的重叠集相加少于 1823 个吗？
我想我无法从示例中推断出您的算法。 n = 20 的重叠集是什么样的？
等等，那个号码是什么？项数是否大于因子集？
@DavidEisenstat 我想我明白了你的意思...2 3 的因子集在是否可以与5 7 组合在一起方面与2 2 3 相同。这就是我所说的 1823 个因子集的意思。【参考方案4】：
这是一种在 5 秒内将 given sequence 提升到 n=62 的方法（经过优化，它在 5 秒内运行 n=75，但请注意我的 second attempt at this problem 效果更好）。我假设问题的模部分只是为了避免随着函数变大而出现数值错误，所以我现在忽略它。
该方法基于这样一个事实，即我们最多可以在一个子集中为每个素数选择一个数字。
我们从第一个素数 2 开始。如果我们不包括 2，那么我们有 1 个该素数的组合。如果我们确实包含 2，那么我们就有与可被 2 整除的数字一样多的组合。 然后我们进入第二个素数 3，并决定是否包含它。如果我们不包括它，我们有这个素数的 1 个组合。如果我们确实包含 2，那么我们就有尽可能多的组合，因为它们可以被 3 整除。 ...等等。以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 选择一个数字（作为脚注，在我可能会回到的第一种方法中，我将其视为选择当我们选择 6 时，因为 3 被减半，所以我使用 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 的路径，还有两条我们选择的路径：通过2→2,4 和通过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：在
虽然这似乎输出了正确的数字，但正如预期的那样，成本太高了。特别是它至少需要指数级时间（很可能是组合时间）。
显然这种方法不能按要求扩展。但是，这里可能有一些东西可以给其他人一些想法（或将这种方法排除在失败之外）。好像有两种可能：
-
由于以下因素的某种组合，可以使这种方法起作用。
我还没有发现一些巧妙的数学优化，它们完全消除了计算。
有一些效率优化，例如使用
bitset 而不是set。
缓存。这似乎是最有希望的，因为有可能将递归调用结构更改为可以增量更新的树结构（其中父子关系表示相乘，兄弟关系表示相加） .
这种方法无法奏效。
有一些方法与此方法基本无关。
我使用的第一种方法可能会奏效。这更像是一个“封闭形式”的解决方案，它在n=17 上非常有效，但在n=18 及更高版本上却失败了，只有少数人出局。我花了很长时间编写模式并试图弄清楚为什么它突然对n=18 失败，但看不到它。我可能会回到这个，或者如果有人感兴趣，我会将其作为替代答案。
编辑：我使用一些技巧进行了一些优化，尽量避免重做现有计算，并且代码速度提高了大约 10 倍。听起来不错，但这只是不断的改进。真正需要的是对这个问题的一些了解——例如。我们可以将#subsets(n+1) 建立在#subsets(n) 的基础上吗？
【讨论】：
我的想法是这样的，但我计算了subset(n) 和subset(n+1) 之间的差异。我的想法是计算：从之前的subset(n) 子集中可以包含多少个子集n+1？如果n+1 是素数，那么答案显然是2*subset(n)，否则，我们需要按照您在此处所做的计算。我很惊讶您实际上可以在 C++ 中用相当短的代码实现这一点。我在 Python 中使用与您的代码长度大致相同的代码实现了这一点。 （顺便说一句，我的机器速度很快，它可以在 5 秒内计算到 n=160）
@justhalf 听起来不错。我确实实现了2*subset(n) 优化（节省了大约 30%）。我怀疑您可能比我做得更好，尤其是在解释 Python 时，可能值得张贴/解释您所做的事情。我认为“你能在 5 秒内做多少”是对算法的合理判断，特别是如果我们无法获得没有指数/组合爆炸的版本。
我对您尝试提供封闭式解决方案更感兴趣。你能告诉我们更多关于它的信息吗？
@justhalf 它可能会被刚刚发布的公式所取代...如果您仍然感兴趣，请明天告诉我，但它有点脆弱而且非常hacky。我刚才写了另一个答案，它使用基于从集合中删除元素的递归公式。不过，我认为您的方法可能更类似于我的第一个答案。【参考方案5】：
我会这样做：
-
找出直到
n 的数字的质因数mod m
创建一个队列q的集合，并将空集合添加到其中，并将计数器设置为1
当队列不为空时，从队列中弹出一个元素X
对于从max(X)到n的所有数字k，检查是否因数
该集合与数字的因子相交。如果没有，请添加到
排队X U k 并将计数器加1。否则，转到下一个
ķ。
退货计数器
必须指出两件重要的事情：
您不需要对直到n 的数字进行因式分解，而只需对它们的素因数进行分解，这意味着，对于 12，您只需要 2 和 3。这样检查 2 个数字是否互质就变成了检查交集是否两套是空的。
您可以跟踪您创建的每个新集合的“因子集合”，这样您就不必将每个新数字与集合中的每个其他数字进行测试，而只需将他的素因子集合与整套之一。
【讨论】：
你还没有定义k来自什么集合，它的初始值，如何得到下一个k等等。还要指定算法复杂度；如果总共有 s 个互质子集并且需要工作 w 来检查相交，则看起来可能是 O(s·w)。由于 s 类似于 O(2ⁿ)，因此您的方法可能很慢。可能存在 O(2ⁿ mod m) 方法。 @jwpat7 k 不是一个集合，只是一个介于 max(X) 和 n 之间的数字。根据我的计算，成本应该是 O(s*n^3)，其中 s 是 2 个子集相交的成本：事实上，如果你考虑最坏的情况，你必须检查 n 个数字大小为 1、n-1 的子代与大小为 to 的子代，依此类推。 如何确保始终有足够的空间来维护该队列？ 问题中提到的'mod m'是针对最终答案，即mod m的集合数。质因数 mod m 会导致其他结果。【参考方案6】：这是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
但既然这是一个竞争问题，就必须有一个更快更肮脏的解决方案。而且这种方法需要指数级的时间和空间，3000以下的素数有430个，也比较优雅。
【讨论】：
您的代码似乎暗示您想让primeSets 包含2^430 个元素，这已经超过了已知宇宙中的原子数。
我没提过吗？ :p
当然，但考虑到所有因素；我宁愿有一个需要大量时间但能够在普通机器上运行的解决方案。如果由于内存不足而无法运行而不崩溃，则该解决方案不是真正的解决方案。【参考方案7】：
编辑：添加了递归方法。在 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)
【讨论】：
预先计算一个矩阵来检查数字是否被纠正是一个非常好的主意；我希望我能想到这一点。不过，我不确定是否要单独枚举每个集合——我认为要获得一种有效的解决方案，您需要以某种方式将它们分组计数。【参考方案8】：这是我之前提到的不同方法。
它确实比我以前使用的快得多。它可以使用在线解释器（ideone）在不到 5 秒的时间内计算出高达 coprime_subsets(117)。
代码从空集开始构建可能的集合，并将每个数字插入所有可能的子集中。
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
我还有另一个想法，如果有时间我会尝试实现。我相信另一种方法可能会快很多。
【讨论】：
添加了第三个版本，使用了不同的方法。【参考方案9】：看起来提议的答案以及问题的序言所解决的问题与所提出的问题不同。
问题是：
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 and 3, 4 mod 7: 2=(3*3)%7 and 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 对数。
【讨论】：
链接描述中对示例测试用例的解释与您的解释相矛盾。 这真是一个巧妙的观察，我赞成完全是因为环算术的概念。但是，如果您查看原始链接的 pdf，您会看到 exact 表示法是：“输出 1, 2, 3, ..., n 的互质子集的数量，模 m" - 模运算符前有逗号。 我不太明白逗号的意义。至于 Z 或模 m 中的算术，两者都对实际应用有意义。例如，如果问题源于密码学，那么询问与密钥乘积取模的互质是完全合理的。鉴于模 m 和示例之间的矛盾，只有提出问题的人才能说出需要回答哪种解释...... 我认为模数是为了避免舍入误差，以便在竞争环境中给出准确的答案并进行评估。 original question的链接很清楚。以上是关于如何计算集合 1,2,3,..,n 的互质子集的数量的主要内容，如果未能解决你的问题，请参考以下文章