Как рассчитать число взаимно простых подмножеств множества {1,2,3,.., n}

Ответ 1

Хорошо, это товар. После этого программа C получает n = 3000 меньше чем 5 секунд для меня. Мои шляпы у команды (-ов), которые решили эту проблема в конкурентной среде.

Алгоритм основан на идее обработки малых и больших простые по разному. Матрица мала, если ее квадрат не превосходит n. В противном случае, его большой. Обратите внимание, что каждое число, меньшее или равное n имеет не более одного большого простого множителя.

Мы делаем таблицу, индексированную парами. Первый компонент каждой пары указывает количество используемых больших простых чисел. Второй компонент каждая пара задает набор используемых небольших простых чисел. Значение в точке конкретным индексом является количество решений с этой моделью использования, а не содержащее 1 или большое простое число (мы подсчитываем их позже, умножая на соответствующая мощность 2).

Мы итерируем вниз по номерам j без больших простых множителей. На начало каждой итерации, таблица содержит подсчеты для подмножеств of 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;
}

Ответ 2

Здесь ответ, который проходит через первые 200 элементов в последовательность менее чем за секунду, давая правильный ответ 200 → 374855124868136960. С оптимизацией (см. Правление 1) он может вычислить первые 500 записей в возрасте до 90 лет, что быстро, хотя ответ @David Eisenstat, вероятно, будет лучше, если его можно будет разработать. Я думаю, что используется другой подход к алгоритмам, данным до сих пор, включая мой собственный оригинальный ответ, поэтому я отправляю его отдельно.

После оптимизации я понял, что действительно кодирую проблему с графом, поэтому переписал решение как реализацию графика (см. правление 2). Реализация графика позволяет еще несколько оптимизаций, намного более изящная, более чем на порядок быстрее и масштабируется лучше: она вычисляет f(600) в 1,5 с по сравнению с 27.

Основная идея здесь заключается в использовании отношения рекурсии. Для любого множества число подмножеств, удовлетворяющих критерию, представляет собой сумму:

• количество подмножеств с удаленным одним элементом; и
• количество подмножеств с определенным элементом.

Во втором случае, когда элемент определенно включен, любые другие элементы, которые несовместимы с ним, должны быть удалены.

Проблемы с эффективностью:

• Я решил удалить самый большой элемент, чтобы максимизировать вероятность того, что этот элемент уже является общим для всех остальных, и в этом случае нужно сделать только один, а не два рекурсивных вызова.
• Кэширование /memoization помогает.

Код ниже.

#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;
}

Характеристики времени выглядят намного лучше, чем мой первый ответ. Но все равно не будет до 3000 в 5 секунд!

Изменить 1

Есть несколько интересных оптимизаций, которые могут быть сделаны для этого метода. В целом это дает 4-кратное улучшение для более крупных n.

• Все числа в наборе, которые уже являются взаимно простыми, могут быть удалены на одном этапе предварительной обработки: если номер m удален, тогда исходный набор имеет 2 ^m фактор раз больше комбинаций, чем уменьшенный один (потому что для каждого совместного использования вы можете либо иметь его в или из набора независимо от других элементов).
• Самое главное, что можно выбрать элемент для удаления, который находится где угодно в наборе. Оказывается, что удаление самого подключенного элемента работает лучше всего.
• Рекурсивная связь, которая ранее использовалась, может быть обобщена, чтобы удалить более одного элемента, где все удаленные элементы имеют одинаковые основные факторы. Например. для набора {2, 3, 15, 19, 45} числа 15 и 45 имеют одинаковые первичные коэффициенты 3 и 5. Одновременно удаляются 2 числа, поэтому число подмножеств для {2, 3, 15, 19, 45}= в два раза больше комбинаций для 15 или 45 присутствует (для набора {2, 19}, потому что 3 должно отсутствовать, если присутствуют 15 или 45) + число подмножеств для 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;
}

Можно обрабатывать до n=600 примерно через 5 минут. Однако время все еще выглядит экспоненциально, удваивая каждые 50-60 n или около того. График вычисления только одного n показан ниже.

Изменить 2

Это решение гораздо более естественно реализовано в терминах графика. Появились еще две оптимизации:

• Самое главное, если граф G можно разбить на два множества A и B таких, что между A и B нет связей, то взаимно однозначно (G) = взаимно однозначно (A) * взаимно однозначно (B).

• Во-вторых, возможно свернуть все числа для набора простых коэффициентов в один node, поэтому значение для node - это счетчик чисел для его простых факторов.

В приведенном ниже коде класс Graph содержит исходную матрицу смежности и значения node, а класс 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, которое занимает около 2 секунд для n = 200 и замедляется экспоненциально.

{-# 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

Ответ 4

Здесь подход, который получает заданную последовательность до n=62 менее чем за 5 секунд (при оптимизации он работает n=75 в 5 секунд, однако обратите внимание на мою вторую попытку этой проблемы). Я предполагаю, что модульная часть проблемы заключается лишь в том, чтобы избежать числовых ошибок, поскольку функция становится большой, поэтому я игнорирую ее пока.

Подход основан на том, что мы можем выбрать не более одного числа в подмножестве для каждого простого числа.

• Начнем с первого простого, 2. Если мы не будем включать 2, то для этого простого числа мы имеем 1 комбинацию. Если мы включим 2, то у нас будет столько комбинаций, сколько чисел, делящихся на 2.
• Затем мы переходим на второе правое число 3 и решаем, включать или не включать это. Если мы не включим его, у нас есть 1 комбинация для этого простого числа. Если мы включим 2, то у нас будет столько комбинаций, сколько чисел, делящихся на 3.
• ... и т.д.

Взяв пример {1,2,3,4}, сопоставим числа в множестве на простые числа следующим образом. Я включил 1 как простое, поскольку это облегчает экспозицию на этом этапе.

1 → {1}
2 → {2,4}
3 → {3}

У нас есть 2 комбинации для "простого" 1 (не включайте его или 1), 3 комбинации для простого числа 2 (не включайте его или 2 или 4) и 2 комбинации для 3 (не включить его или 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}

но если мы выберем число 6 для простого числа 2, мы не можем выбрать число для простого числа 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, и два пути, где мы делаем: через {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.

Код С++ для рекурсивного запуска через пути приведен ниже. (Это С++ 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;
}

Сроки: вычисляет последовательность до 40 в < 0,1 с, последовательность до 50 в 0,5 с, до 60 в 2,5 с, до 70 в 20 с и до 80 в 157 с.

Несмотря на то, что это, несомненно, выводит правильные цифры, это, как и следовало ожидать, слишком дорого. В частности, требуется, по крайней мере, экспоненциальное время (и, вполне возможно, комбинаторное время).

Очевидно, что этот подход не масштабируется по мере необходимости. Однако здесь может быть что-то, что дает другим людям идеи (или исключает этот подход как провал). Кажется, есть две возможности:

• Этот подход может быть выполнен для работы из-за некоторой комбинации следующего.
  • Есть несколько умных математических оптимизаций, которые я не заметил, которые полностью исключают вычисления.
  • Существуют некоторые оптимизации эффективности, например. используйте bitset, а не set.
  • Caching. Это кажется наиболее перспективным, поскольку можно было бы изменить структуру рекурсивного вызова в древовидную структуру, которая может быть инкрементно обновлена ​​ (где отношение родительский и дочерний указывает на умножение, а отношение родственников указывает на добавление),
    
• Этот подход не может работать.
  • Существует некоторый подход, который в значительной степени не связан с этим.
  • Возможно, что первый подход, который я использовал, можно заставить работать. Это было гораздо более "закрытым" решением, которое эффективно работало для upto n=17 и не удавалось при n=18 и выше, выбывая на небольшое число. Я долгое время писал шаблоны и пытался выяснить, почему он неожиданно потерпел неудачу в отношении n=18, но не смог его увидеть. Я могу вернуться к этому, или я включу его в качестве альтернативного ответа, если кто-то заинтересован.

Изменить. Я сделал некоторые оптимизации, используя несколько трюков, чтобы избежать повторения существующих вычислений там, где это возможно, и код примерно в 10 раз быстрее. Звучит неплохо, но это только постоянное улучшение. Что действительно нужно, так это некоторое понимание этой проблемы - например, можем ли мы основать #subsets(n+1) на #subsets(n)?

Ответ 5

Вот как я это сделаю:

• Найдите простые коэффициенты mod m чисел до n
• Создайте очередь q наборов и добавьте к нему пустой набор и установите счетчик 1
• Пока очередь не пуста, вытащите элемент X из очереди
• Для всех чисел k от max(X) до n, проверьте, являются ли факторы множество пересекается с множителями числа. Если нет, добавьте queue X U k и увеличивать счетчик на 1. В противном случае перейдите к следующему к.
• Возвращаемый счетчик

Следует отметить две важные вещи:

• Вам не нужна факторизация чисел до n, а просто их основные факторы, это означает, что для 12 вам просто нужны 2 и 3. Таким образом, проверка того, что 2 числа взаимно просты, проверяет, пересекает ли пересечение из двух наборов пуст.
• Вы можете отслеживать "набор факторов" каждого нового набора, который вы создаете, таким образом, вам не придется проверять каждый новый номер на каждый другой номер в наборе, а просто пересекать его простые коэффициенты, установленные против один из целого набора.

Ответ 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
    }
}

Но поскольку это проблема конкуренции, должно быть более быстрое и более грязное решение. И поскольку этот метод требует экспоненциального времени и пространства, и есть 430 простых чисел до 3000, что-то более элегантное тоже.

Ответ 7

EDIT: добавлен рекурсивный подход. Решает за 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 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: построение булевой матрицы
1. Составьте список простых чисел от 2 до 1500.
2. Для каждого номера от 1 до 3000 постройте набор его основных факторов.
3. Сравните каждую пару множеств и получите булевую матрицу [3000] [3000], в которой указывается, являются ли элементы я и j взаимно взаимными (1) или нет (0).

Шаг 2: Рассчитайте количество взаимно простых множеств длины k (k = 0 до 3000)
1. Инициализировать count = 1 (пустой набор). Теперь k = 1. Ведем список множеств длины k.
2. Создайте 3000 наборов, содержащих только этот конкретный элемент. (увеличение счетчика)
3. Сканируйте каждый элемент с до 3000 и посмотрите, можно ли создать новый набор, добавив его к любому из существующих наборов длины k. Примечание: некоторые вновь сформированные наборы могут быть или не быть идентичными. Если вы используете набор наборов, то сохраняются только уникальные множества.
4. Удалить все наборы, все еще имеющие длину k.
5. Инкремент рассчитывается по текущему количеству уникальных наборов.
6. k = k + 1 и перейти к шагу 3.

В качестве альтернативы вы можете поддерживать список продуктов каждого из элементов в наборе и проверять, является ли новый элемент совпадающим с продуктом. В этом случае вам не нужно хранить логическую матрицу.

Сложность для вышеуказанного алгоритма кажется где-то между O (n ^ 2) и O (n ^ 3).

Полный код в С++: (улучшение: условие добавлено, что элемент должен быть проверен в наборе, только если он больше чем наибольшее значение в наборе).

#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 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

Вот другой подход, который я упомянул ранее.
Это действительно намного быстрее, чем тот, который я использовал раньше. Он может рассчитывать до coprime_subsets(117) менее чем за 5 секунд, используя онлайн-интерпретатор (ideone).

В коде строятся возможные множества, начиная с пустого набора, и вставляя каждое число во все возможные подмножества.

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, потому что оба они делятся на 2 по модулю 7: 4 = (2 * 2)% 7 и 3 = (2 * 5)% 7.

На самом деле, арифметика Z_m настолько странна, что можно дать ответ в O (1) раз, по крайней мере для простого m: для любого n и простого m не существует NO взаимно простых пар по модулю m. Вот почему: ненулевые элементы Z_m образуют циклическую группу порядка m-1, откуда следует, что для любых ненулевых элементов a и b из Z_m можно записать a = bc для некоторого c из Z_m. Это доказывает отсутствие взаимных пар в 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, и поэтому a будет делиться на b в Z_m.

Таким образом, остается рассмотреть случай, когда m является составным и a и b делятся на различные простые множители m, скажем, a делится на p и b делится на q. В этом случае они действительно могут быть взаимными. Например, 2 и 3 по модулю 6 являются взаимными.

Таким образом, вопрос о взаимных парах сводится к следующим шагам:

• Нахождение простых множителей m, меньших n. Если нет, то нет взаимно простых пар.

• Перечисление произведений степеней этих простых множителей, которые остаются факторами m, меньшими n.

• Вычисление числа пар Z-содержащихся среди thoses.