使用容斥原理计素数个数

count-primes.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define BOUND 2000000000

int primes_count = 0;
int primes[5000] = {0};
char is_composites[45000] = {0}; // > sqrt(2e9)

void find_primes(int limit) {
  // 筛法
  primes_count = 1;
  for (int n = 3; n <= limit; n += 2) {
    if (!is_composites[n]) {
      primes_count++;
      for (int index = n * n; index <= limit; index += n) {
        is_composites[index] = 1;
      }
    }
  }

  // 收集素数
  primes[0] = 2;
  int prime_index = 1;
  for (int index = 3; index <= limit; index += 2) {
    if (!is_composites[index]) {
      primes[prime_index] = index;
      prime_index++;
    }
  }
}

// 求任意size个素数的积在BOUND内倍数的个数
long long combinate(int size, long long product, int since) {
  if (size > 0) {
    long long total = 0;
    for (int index = since; index <= primes_count - size; index++) {
      long long count = combinate(size - 1, product * primes[index], index + 1);
      if (count > 0) {
        total += count;
      } else {
        break;
      }
    }
    return total;
  } else if (product <= BOUND) {
    return BOUND / product;
  } else {
    return 0;
  }
}

// 根据容斥原理求素数的个数
long long inclusion_exclusion() {
  long long total = BOUND - 1 + primes_count;
  long long sign = -1;
  for (int size = 1; size <= primes_count; size++) {
    long long value = combinate(size, 1, 0);
    if (value > 0) {
      total += sign * value;
      sign = -sign;
    } else {
      break;
    }
  }
  return total;
}

int main(void) {
  find_primes((int)ceil(sqrt(BOUND)));
  printf("%lld\n", inclusion_exclusion());

  return EXIT_SUCCESS;
}