如何计算集合{1,2,3,..,n}的互质子集数量
Lov*_*per 25 algorithm primes integer
我正在解决这个问题(问题一).声明是:
该集合中{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}.
我认为可以通过使用素数来解决.(跟踪我们是否使用了每个素数)..但我不确定.
我可以得到一些解决这个任务的提示吗？
- 您可以在此处找到此序列:http://oeis.org/A084422
Dav*_*tat 16
好的,这是货.随后的C程序在不到5秒的时间内获得n = 3000.我的帽子是在竞争激烈的环境中解决这个问题的团队.
该算法基于以不同方式处理小和大素数的想法.如果素数最多为n,则素数很小.否则,它很大.观察到小于或等于n的每个数字最多具有一个大的素因子.
我们创建一个由成对索引的表.每对的第一个组件指定使用的大质数的数量.每对的第二个组件指定使用的小素数集.特定索引处的值是该使用模式不包含1或大素数的解的数量(我们稍后通过乘以适当的2的幂来计算它们).
我们向下迭代数j,没有大的素因子.在每次迭代开始时,该表包含j..n子集的计数.内循环中有两个添加项.第一个考虑通过j本身扩展子集,这不会增加使用中的大质数的数量.第二个用于将子集扩展j次大素数,这样做.合适的大素数的数量是不大于n/j的大素数的数量减去已经使用的大素数的数量,因为向下迭代意味着已经使用的每个大素数不大于n/j.
最后,我们对表条目求和.在表中计数的每个子集产生2**k个子集,其中k是1加上未使用的大质数的数量,为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？通过优化(参见编辑1),它可以计算90岁以下的前500个条目,这很快 - 尽管如果可以开发@David Eisenstat的答案可能会更好.我认为到目前为止给出的算法采用了不同的方法,包括我自己的原始答案,所以我将它单独发布.
在优化之后,我意识到我真的在编写图形问题,所以我将解决方案重写为图形实现(参见编辑2).图形实现允许一些更优化,更优雅,更快一个数量级,并且更好地扩展:它计算f(600)在1.5s,而27s.
The main idea here is to use a recursion relationship. For any set, the number of subsets meeting the criterion are the sum of:
- the number of subsets with one element removed; and
- the number of subsets with that element definitely included.
In the second case, when the element is definitely included, any other elements which aren't coprime with it must be removed.
Efficiency issues:
- I've chosen to remove the largest element to maximize the chance of that element already being coprime to all the others, in which case only one rather than two recursive calls need to be made.
- Caching/memoization helps.
Code below.
#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;
}
Time characteristics look a lot better than my first answer. But still won't go upto 3000 in 5s!
Edit 1
There are some interesting optimizations that can be made to this method. Overall this gives a 4x improvement for larger n.
- All numbers in the set that are already coprime can be removed in one single preprocessing step: if
mnumber are removed, then the original set has 2m factor times more combinations than the reduced one (because for each coprime, you can either have it in or out of the set independently of other elements). - Most importantly, it is possible to choose an element to remove that is anywhere in the set. It turns out that removing the most connected element works best.
- The recursive relationship that was previously used can be generalized to remove more than one element where all elements removed have the same prime factors. E.g. for the set
{2, 3, 15, 19, 45}, the numbers 15 and 45 have the same prime factors of 3 and 5. There are 2 numbers removed at once, and so the number of subsets for{2, 3, 15, 19, 45}= twice the number of combinations for either 15 or 45 present (for the set{2, 19}because 3 has to be absent if either 15 or 45 are present) + the number of subsets for 15 and 45 absent (for the set{2, 3, 19}) - Using
shortfor the number type improved performance by about 10%. - Finally, it is also possible to transform sets into sets with equivalent prime factors, in the hope of getting better cache hits by standardizing the sets. E.g.,
{ 3, 9, 15}is equivalent (isomorphic) to2, 4, 6. This was the most radical idea but probably had the least effect on performance.
It's probably a lot easier to understand it with a concrete example. I've chosen the set {1..12} which is large enough to get a feel for how it works but small enough so that it's comprehensible.
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
Code below
#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;
}
It's possible to process upto n=600 in about 5 minutes. Time still looks exponential however, doubling every 50 to 60 n or so. The graph for calculating just one n is shown below.
Edit 2
This solution is much more naturally implemented in terms of a graph. Two more optimizations arose:
最重要的是,如果图G可以被分成两组A和B,使得A和B之间没有连接,那么互质(G)=互质(A)*互质(B).
其次,可以将一组素因子的所有数字折叠成单个节点,因此节点的值是其素因子的数量计数.
在下面的代码中,Graph类保存原始邻接矩阵和节点值,并且FilteredGraph类将剩余节点的当前列表保存为a,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;
}
计算coprimes(n)的图表如下所示为红色(旧方法为黑色).
基于观察到的(指数)增长率n=3000,假设程序没有爆炸,预测为30小时.这开始看起来在计算上是可行的,特别是在进行更多优化的情况下,但是远远不能满足要求的5s!毫无疑问,所需的解决方案简短而甜蜜,但这很有趣......