Search code examples
algorithmintegerprimes

How to calculate the number of coprime subsets of the set {1,2,3,..,n}


I am solving this task (problem I). The statement is:

How many subsets of the set {1, 2, 3, ..., n} are coprime? A set of integers is called coprime if every two of its elements are coprime. Two integers are coprime if their greatest common divisor equals 1.

Input

First line of input contains two integers n and m (1 <= n <= 3000, 1 <= m <= 10^9 + 9)

Output

Output the number of coprime subsets of {1, 2, 3, ..., n} modulo m.

Example

input: 4 7 output: 5

There are 12 coprime subsets of {1,2,3,4}: {}, {1}, {2}, {3}, {4}, {1,2}, {1,3}, {1,4}, {2,3}, {3,4}, {1,2,3}, {1,3,4}.


I think it can be solved by using prime numbers. (keeping track of if we used each prime numbers) ..but I'm not sure.

Can I get some hints to solve this task?


Solution

  • Okay, here’s the goods. The C program that follows gets n=3000 in less than 5 seconds for me. My hat’s off to the team(s) that solved this problem in a competitive setting.

    The algorithm is based on the idea of treating the small and large primes differently. A prime is small if its square is at most n. Otherwise, it’s large. Observe that each number less than or equal to n has at most one large prime factor.

    We make a table indexed by pairs. The first component of each pair specifies the number of large primes in use. The second component of each pair specifies the set of small primes in use. The value at a particular index is the number of solutions with that usage pattern not containing 1 or a large prime (we count those later by multiplying by the appropriate power of 2).

    We iterate downward over numbers j with no large prime factors. At the beginning of each iteration, the table contains the counts for subsets of j..n. There are two additions in the inner loop. The first accounts for extending subsets by j itself, which does not increase the number of large primes in use. The second accounts for extending subsets by j times a large prime, which does. The number of suitable large primes is the number of large primes not greater than n/j, minus the number of large primes already in use, since the downward iteration implies that each large prime already in use is not greater than n/j.

    At the end, we sum the table entries. Each subset counted in the table gives rise to 2 ** k subsets where k is one plus the number of unused large primes, as 1 and each unused large prime can be included or excluded independently.

    /* 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;
    }