algorithm - tiene - subconjuntos propios
Cómo calcular el número de subconjuntos de coprime del conjunto{1,2,3,.., n} (9)
Estoy resolviendo esta tarea (problema I) . La declaración es:
¿Cuántos subconjuntos del conjunto {1, 2, 3, ..., n} son coprime? Un conjunto de enteros se llama coprime si cada dos de sus elementos es coprime. Dos enteros son coprime si su mayor divisor común es igual a 1.
Entrada
La primera línea de entrada contiene dos enteros n y m ( 1 <= n <= 3000, 1 <= m <= 10^9 + 9 )
Salida
Obtenga el número de subconjuntos de coprime de {1, 2, 3, ..., n} módulo m .
Ejemplo
entrada: 4 7 salida: 5
Hay 12 subconjuntos de coprime de {1,2,3,4} : {} , {1} , {2} , {3} , {4} , {1,2} , {1,3} , {1,4} , {2,3} , {3,4} , {1,2,3} , {1,3,4} .
Creo que se puede resolver utilizando números primos. (mantener un registro de si usamos cada uno de los números primos) .. pero no estoy seguro
¿Puedo obtener algunos consejos para resolver esta tarea?
- Puede encontrar esta secuencia aquí: http://oeis.org/A084422
Aquí hay algo bastante sencillo en Haskell, que tarda unos 2 segundos para n = 200 y se ralentiza exponencialmente.
{-# 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
Aquí hay un enfoque que lleva la this hasta n=62 en menores de 5 años (con optimizaciones se ejecuta n=75 en 5 segundos , sin embargo, tenga en cuenta que mi segundo intento de este problema funciona mejor). Supongo que la parte de módulo del problema solo tiene que ver con evitar los errores numéricos a medida que la función crece, por lo que por ahora lo ignoro.
El enfoque se basa en el hecho de que podemos elegir a lo sumo un número en un subconjunto para cada primo.
- Comenzamos con la primera prima, 2. Si no incluimos 2, entonces tenemos 1 combinación para esta prima. Si incluimos 2, entonces tenemos tantas combinaciones como son números divisibles por 2.
- Luego vamos a la segunda prima, 3, y decidimos si incluir o no eso. Si no lo incluimos, tenemos 1 combinación para este primer. Si incluimos 2, entonces tenemos tantas combinaciones como son números divisibles por 3.
- ... y así.
Tomando el ejemplo {1,2,3,4} , asignamos los números del conjunto a los números primos de la siguiente manera. He incluido 1 como primo, ya que facilita la exposición en esta etapa.
1 → {1}
2 → {2,4}
3 → {3}
Tenemos 2 combinaciones para el "prime" 1 (no lo incluya o 1), 3 combinaciones para el prime 2 (no lo incluya o 2 o 4), y 2 combinaciones para 3 (no lo incluya o 3). Entonces el número de subconjuntos es 2 * 3 * 2 = 12 .
Del mismo modo para {1,2,3,4,5} tenemos
1 → {1}
2 → {2,4}
3 → {3}
5 → {5}
dando 2 * 3 * 2 * 2= 24 .
Pero para {1,2,3,4,5,6} , las cosas no son tan sencillas. Tenemos
1 → {1}
2 → {2,4,6}
3 → {3}
5 → {5}
pero si elegimos el número 6 para el prime 2, no podemos elegir un número para el prime 3 (como nota al pie de página, en mi primer acercamiento, al que puedo volver, traté esto como si las opciones para 3 fueran cut in half when we chose 6, so I used 3.5 rather than 4 for the number of combinations for the prime 2 and 2 * 3.5 * 2 * 2 = 28 gave the right answer. I couldn''t get this approach to work beyond n=17 , however.)
La forma en que lidié con esto es dividir el procesamiento de cada conjunto de factores primos en cada nivel. Así que {2,4}tienen factores primos {2}, mientras que {6}tiene factores primos {2,3}. Omitiendo la entrada espuria para 1 (que no es un número primo), ahora tenemos
2 → {{2}→{2,4}, {2,3}→6}
3 → {{3}→{3}}
5 → {{5}→{5}}
Ahora hay tres rutas para calcular el número de combinaciones, una ruta donde no seleccionamos el Prime 2 y dos rutas donde lo hacemos: a través {2}→{2,4}y por medio {2,3}→6.
- La primera ruta tiene
1 * 2 * 2 = 4combinaciones porque podemos seleccionar 3 o no, y podemos seleccionar 5 o no. - Del mismo modo, el segundo tiene
2 * 2 * 2 = 8combinaciones, teniendo en cuenta que podemos elegir entre 2 o 4. - El tercero tiene
1 * 1 * 2 = 2combinaciones, porque solo tenemos una opción para el Prime 3: no usarlo.
En total, esto nos da 4 + 8 + 2 = 14combinaciones (como una nota de optimización de que la primera y la segunda rutas podrían haberse colapsado juntas para obtener 3 * 2 * 2 = 12). También tenemos la opción de seleccionar 1 o no, por lo que el número total de combinaciones es 2 * 14 = 28.
El código C ++ para ejecutar recursivamente a través de las rutas se encuentra a continuación. (Es C ++ 11, escrito en Visual Studio 2012, pero debería funcionar en otros gcc ya que no he incluido nada específico de la plataforma).
#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;
}
Tiempos: calcula la secuencia hasta 40 en <0.1s, la secuencia hasta 50 en 0.5s, a 60 en 2.5s, a 70 en 20s, y a 80 en 157s.
Aunque esto ciertamente parece dar los números correctos, es, como podría esperarse, demasiado costoso. En particular, toma al menos tiempo exponencial (y posiblemente tiempo combinatorio).
Claramente este enfoque no se escala como se requiere . Sin embargo, puede haber algo aquí que le da ideas a otras personas (o descarta este enfoque como un fracaso). Parece que hay dos posibilidades:
- Se puede hacer que este enfoque funcione, debido a alguna combinación de lo siguiente.
- Hay algunas optimizaciones matemáticas inteligentes que no he detectado que eliminan los cálculos por completo.
- Hay algunas optimizaciones de eficiencia, por ejemplo, uso en
bitsetlugar deset. - Caching Esto parece más prometedor, ya que podría ser posible cambiar la estructura de la llamada recursiva en una estructura de árbol, que podría actualizarse de forma incremental (donde una relación padre-hijo indica multiplicar, y la relación entre hermanos indica agregar).
- Este enfoque no puede hacerse funcionar.
- Hay algún enfoque que no está relacionado en gran medida con este.
- Es posible que el primer enfoque que utilicé se hiciera funcionar. Esto fue mucho más que una solución de "forma cerrada" que funcionó de manera muy eficiente hasta el momento
n=17y fracasó enn=18un principio. Pasé mucho tiempo escribiendo los patrones y tratando de averiguar por qué falló repentinamenten=18, pero no pude verlo. Puedo volver a esto, o lo incluiré como una respuesta alternativa si alguien está interesado.
Edición : He realizado algunas optimizaciones utilizando algunos trucos para evitar rehacer los cálculos existentes cuando sea posible y el código es aproximadamente 10 veces más rápido. Suena bien, pero es solo una mejora constante . Lo que realmente se necesita es una idea de este problema - por ejemplo, podemos basar #subsets(n+1)en #subsets(n)?
Aquí hay una respuesta que recorre los primeros 200 elementos de la this en menos de un segundo, dando la respuesta correcta 200 → 374855124868136960. Con optimizaciones (ver edición 1), puede calcular las primeras 500 entradas en menos de 90, lo cual es rápido: - Aunque la respuesta de @David Eisenstat es probable que sea mejor si se puede desarrollar. Creo que tiene un enfoque diferente a los algoritmos dados hasta ahora, incluida mi propia respuesta original, así que lo publico por separado.
Después de optimizar, me di cuenta de que realmente estaba codificando un problema gráfico, por lo que reescribí la solución como una implementación gráfica (ver edición 2). La implementación del gráfico permite algunas optimizaciones más, es mucho más elegante, más de un orden de magnitud más rápido y se escala mejor: calcula f(600) en 1.5 s, en comparación con 27s.
La idea principal aquí es usar una relación de recursión. Para cualquier conjunto, el número de subconjuntos que cumplen el criterio son la suma de:
- el número de subconjuntos con un elemento eliminado; y
- el número de subconjuntos con ese elemento definitivamente incluido.
En el segundo caso, cuando el elemento se incluye definitivamente, cualquier otro elemento que no sea coprime debe eliminarse.
Cuestiones de eficiencia:
- He elegido eliminar el elemento más grande para maximizar la posibilidad de que ese elemento ya sea coprime a todos los demás, en cuyo caso solo se debe hacer una en lugar de dos llamadas recursivas.
- Caching / memoization ayuda.
Código abajo.
#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;
}
Las características del tiempo se ven mucho mejor que mi primera respuesta . ¡Pero todavía no subirás hasta 3000 en 5s!
Editar 1
Hay algunas optimizaciones interesantes que se pueden hacer a este método. En general, esto da una mejora de 4x para n más grande.
- Todos los números en el conjunto que ya son coprime pueden eliminarse en un solo paso de preprocesamiento: si se elimina
mnúmero, entonces el conjunto original tiene 2 m factor multiplicado por más combinaciones que el reducido (porque para cada coprime, puede tenerlo) dentro o fuera del conjunto independientemente de otros elementos). - Lo más importante es que es posible elegir un elemento para eliminar que esté en cualquier parte del conjunto. Resulta que quitar el elemento más conectado funciona mejor.
- La relación recursiva que se utilizó anteriormente se puede generalizar para eliminar más de un elemento donde todos los elementos eliminados tienen los mismos factores primos. Por ejemplo, para el conjunto
{2, 3, 15, 19, 45}, los números 15 y 45 tienen los mismos factores primos de 3 y 5. Hay 2 números eliminados a la vez, y así el número de subconjuntos para{2, 3, 15, 19, 45}= el doble del número de combinaciones para 15 o 45 presentes (para el conjunto{2, 19}porque 3 deben estar ausentes si hay 15 o 45 presentes) + el número de subconjuntos para 15 y 45 ausentes (para el conjunto{2, 3, 19}) - El uso de
shortpara el tipo de número mejoró el rendimiento en aproximadamente un 10%. - Finalmente, también es posible transformar conjuntos en conjuntos con factores primos equivalentes, con la esperanza de obtener mejores resultados de caché al estandarizar los conjuntos. Por ejemplo,
{ 3, 9, 15}es equivalente (isomorfo) a2, 4, 6. Esta fue la idea más radical, pero probablemente tuvo el menor efecto en el rendimiento.
Probablemente sea mucho más fácil de entender con un ejemplo concreto. He elegido el conjunto {1..12} que es lo suficientemente grande para tener una idea de cómo funciona, pero lo suficientemente pequeño para que sea comprensible.
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
Código abajo
#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;
}
Es posible procesar hasta n=600 en aproximadamente 5 minutos. El tiempo todavía parece exponencial, sin embargo, se duplica cada 50 a 60 no más. La gráfica para calcular solo una n se muestra a continuación.
Editar 2
Esta solución se implementa de forma mucho más natural en términos de una gráfica. Surgieron dos optimizaciones más:
Lo más importante es que si el gráfico G se puede dividir en dos conjuntos A y B de manera que no haya conexiones entre A y B, entonces coprimes (G) = coprimes (A) * coprimes (B).
En segundo lugar, es posible colapsar todos los números para un conjunto de factores primos en un solo nodo, por lo que el valor para el nodo es el conteo de números para sus factores primos.
En el código a continuación, la clase Graph contiene la matriz de adyacencia original y los valores de los nodos, y la clase FilteredGraph contiene la lista actual de nodos restantes como un conjunto de bitset para que, al eliminarse los nodos, la nueva matriz de adyacencia se pueda calcular mediante el enmascaramiento de bits hay relativamente pocos datos para pasar en la recursión).
#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;
}
La gráfica para calcular los coprimes ( n ) se muestra a continuación en rojo (con el enfoque anterior en negro).
Basado en la tasa (exponencial) de aumento observada, la predicción para n=3000 es de 30 horas, asumiendo que el programa no explota. Esto está empezando a parecer computacionalmente factible, especialmente con más optimizaciones, ¡pero no se acerca a los 5s que se requieren! Sin duda, la solución requerida es breve y dulce, pero esto ha sido divertido ...
Así es como lo haría:
- Encuentra los factores primos
mod mde los números hastan - Cree una cola
qde conjuntos, agregue el conjunto vacío y establezca el contador en 1 - Mientras que la cola no está vacía, saca un elemento
Xde la cola - Para todos los números
kdemax(X)an, verifique si los factores del conjunto se intersecan con los factores del número. Si no, agregue a la colaXU ky aumente el contador en 1. De lo contrario, vaya a la siguiente k. - Contador de devoluciones
Hay que señalar dos cosas importantes:
- No necesitas la factorización de los números hasta
n, pero solo sus factores primos, eso significa que para 12 solo necesitas 2 y 3. De esta manera, comprobar si 2 números son coprime se convierte en verificar si la intersección de dos conjuntos está vacía. - Puede realizar un seguimiento del "conjunto de factores" de cada nuevo conjunto que cree, de esta manera no tendrá que probar cada nuevo número contra cada otro número en el conjunto, sino que solo se unirá a sus factores primos frente al de los todo el juego.
Bien, aquí están las mercancías. El programa C que sigue obtiene n = 3000 en menos de 5 segundos para mí. Me quito el sombrero ante el (los) equipo (s) que resolvieron este problema en un entorno competitivo.
El algoritmo se basa en la idea de tratar los primos pequeños y grandes de manera diferente. Un cebado es pequeño si su cuadrado es como máximo n. De lo contrario, es grande . Observe que cada número menor o igual a n tiene a lo sumo un factor primo grande.
Hacemos una tabla indexada por pares. El primer componente de cada par especifica el número de primos grandes en uso. El segundo componente de cada par especifica el conjunto de primos pequeños en uso. El valor en un índice particular es el número de soluciones con ese patrón de uso que no contiene 1 o un número primo grande (los contamos más adelante multiplicándolos por la potencia adecuada de 2).
Recorremos números a la baja sin grandes factores primos. Al comienzo de cada iteración, la tabla contiene los recuentos para los subconjuntos de j..n. Hay dos adiciones en el bucle interno. La primera cuenta para extender los subconjuntos por j mismo, lo que no aumenta el número de primos grandes en uso. La segunda cuenta para extender los subconjuntos en j veces una gran prima, lo que hace. El número de primos grandes adecuados es el número de primos grandes que no son mayores que n / j, menos el número de primos grandes que ya están en uso, ya que la iteración hacia abajo implica que cada cebado grande que ya está en uso no es mayor que n / j.
Al final, sumamos las entradas de la tabla. Cada subconjunto contado en la tabla da lugar a 2 ** k subconjuntos donde k es uno más el número de primos grandes no utilizados, ya que 1 y cada primo grande no utilizado pueden incluirse o excluirse independientemente.
/* 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;
}
Parece que las respuestas propuestas, así como el preámbulo de la pregunta, abordan una pregunta diferente a la formulada.
La pregunta era:
Output the number of coprime subsets of {1, 2, 3, ..., n} modulo m.
Las respuestas propuestas intentan abordar otra:
Output the number of coprime subsets of {1, 2, 3, ..., n}.
Estas preguntas NO son equivalentes. El primero trata con el anillo finito Z_m, y el segundo trata con el anillo de los números enteros Z que tiene una aritmética completamente diferente.
Además, la definición "Dos enteros son coprimo si su mayor divisor común es igual a 1" en el preámbulo de la pregunta no es aplicable a Z_m: los anillos finitos no tienen una comparación aritméticamente estable, por lo que no existe tal cosa como un "mayor" común divisor.
La misma objeción se aplica al ejemplo en la pregunta: 3 y 4 NO son un módulo 7 primo porque los dos son divisibles por 2 módulo 7: 4 = (2 * 2)% 7 y 3 = (2 * 5)% 7.
De hecho, la aritmética Z_m es tan extraña que se puede dar la respuesta en tiempo O (1), al menos para prime m: para cualquier n y prime m NO hay pares de coprime modulo m. Aquí se explica por qué: los elementos que no son cero de Z_m forman un grupo cíclico de orden m-1, lo que implica que para cualquier elemento que no sea cero a y b de Z_m, se puede escribir a = bc para algunas c en Z_m. Esto prueba que no hay pares de coprime en Z_m para prime m.
Del ejemplo de la pregunta: echemos un vistazo a {2, 3} mod 7 y {3, 4} mod 7: 2 = (3 * 3)% 7 y 3 = (4 * 6)% 7. Por lo tanto, {2,3} no son coprime en Z_7 (ambos son divisibles por 3) y {3,4} no son coprime en Z_7 (ambos son divisibles por 4).
Ahora vamos a considerar el caso de m no prime. Escribe ma como un producto de las potencias primarias m = p_1 ^ i_1 * ... * p_k ^ i_k. Si a y b tienen un factor primo común, claramente no son coprime. Si al menos uno de ellos, digamos b, no divide ninguno de los números primos p_1, ..., p_k entonces a y b tienen un factor común por aproximadamente la misma razón que en el caso de prime m: b sería un factor multiplicativo unidad de Z_m, y por lo tanto a sería divisible por b en Z_m.
Entonces, queda por considerar el caso cuando m es compuesto y a y b son divisibles por distintos factores primos de m, digamos que a es divisible por p y b es divisible por q. En este caso sí pueden ser coprimes. Por ejemplo, 2 y 3 módulo 6 son coprimes.
Así que la pregunta para los pares de coprime se reduce a estos pasos:
Encontrar factores primos de m que sean menores que n. Si no hay, entonces no hay pares de coprime.
Enumerar los productos de los poderes de estos factores primos que siguen siendo los factores de m que son menores que n.
Cálculo del número de pares Z-comprime entre estos.
Aquí está el enfoque diferente que mencioné anteriormente.
De hecho, es mucho más rápido que el que he usado antes. Puede calcular hasta coprime_subsets(117)en menos de 5 segundos, utilizando un intérprete en línea (ideone).
El código crea los conjuntos posibles a partir del conjunto vacío e insertando cada número en todos los subconjuntos posibles.
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
Aquí está mi idea y una implementación en python. Parece ser lento, pero no estoy seguro de si es solo la forma en que estaba probando (usando un intérprete en línea) ...
Podría ser que ejecutarlo en una computadora "real" pueda hacer una diferencia, pero puedo '' t prueba que en este momento.
Hay dos partes de este enfoque:
- Pre-generar una lista de factores primos
- Cree una función recursiva en caché para determinar el número de subconjuntos posibles:
- Para cada número, verifique sus factores para ver si se puede agregar al subconjunto
- Si se puede agregar, obtenga el recuento para el caso recursivo y sume al total
Después de eso, supongo que simplemente tomas el módulo ...
Aquí está mi implementación de python (versión mejorada):
# 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
También tengo otra idea, que intentaré implementar si tengo tiempo. Creo que un enfoque diferente podría ser un poco más rápido.
Aquí hay una forma en O (n * 2 ^ p), donde pestá el número de primos debajo n. No hacer uso del módulo.
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
}
}
Pero como este es un problema de competencia, debe haber una solución más rápida y sucia. Y dado que este método necesita tiempo y espacio exponenciales y hay 430 números primos por debajo de 3000, algo más elegante también.
EDITAR: Un enfoque recursivo añadido. Resuelve en 5 segundos hasta 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;
}
Un enfoque ingenuo puede ser (para N = 3000):
Paso 1: Construye una matriz booleana
1. Construye una lista de números primos del 2 al 1500.
2. Para cada número 1 al 3000, crea un conjunto de sus factores primos.
3. Compare cada par de conjuntos y obtenga una matriz booleana [3000] [3000] que indica si los elementos i y j son mutuamente coprime (1) o no (0).
Paso 2: Calcule el número de conjuntos de coprime de longitud k (k = 0 a 3000)
1. Inicialice el conteo = 1 (conjunto vacío). Ahora k = 1. Mantener una lista de conjuntos de longitud k.
2. Construye 3000 conjuntos que contengan solo ese elemento en particular. (incremente la cuenta)
3. Escanee cada elemento de k a 3000 y vea si se puede formar un nuevo conjunto agregándolo a cualquiera de los conjuntos de longitud k existentes. Nota: algunos conjuntos recién formados pueden o no ser idénticos . Si utiliza conjuntos de conjuntos, solo se almacenarán los conjuntos únicos.
4. Eliminar todos los conjuntos que aún tienen una longitud k .
5. Incrementa el conteo por el número actual de conjuntos únicos.
6. k = k + 1 y goto paso 3.
Alternativamente, puede mantener una lista de productos de cada uno de los elementos en un conjunto y verificar si el nuevo elemento es coprime con el producto. En ese caso, no es necesario almacenar la matriz booleana.
La complejidad del algoritmo anterior parece estar entre O (n ^ 2) y O (n ^ 3).
Código completo en C ++: (mejora: condición agregada: ese elemento debe verificarse en un conjunto solo si es> mayor que el valor más grande del conjunto).
#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;
}
El código anterior parece dar un resultado correcto pero consume mucho tiempo: Supongamos que M es lo suficientemente grande:
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)