#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#include <string.h>
#include <errno.h>
#include <stdlib.h>
#include <math.h>

/* This file is C99. Here is some information about compiling this file:
 * $ (source /etc/os-release; echo $PRETTY_NAME); uname -a
 * Debian GNU/Linux trixie/sid
 * Linux redacted 6.11.6-amd64 #1 SMP PREEMPT_DYNAMIC Debian 6.11.6-1 (2024-11-04) x86_64 GNU/Linux
 * $ gcc -ocalc_pi_real calc_pi.c -Wall -Wextra -pedantic -lm -std=c99 -mavx2 -DPOPCOUNT=real_pc -O3
 * $ gcc -ocalc_pi_fake calc_pi.c -Wall -Wextra -pedantic -lm -std=c99 -O3
 * $ diff -s <(./calc_pi_real 1000) <(./calc_pi_fake 1000)
 * Files /dev/fd/63 and /dev/fd/62 are identical
 * $
 */

#ifdef __GNUC__
#define GCC_ALWAYS_INLINE inline __attribute__((always_inline))
#else
#define GCC_ALWAYS_INLINE inline
#endif

/* finds the popcount of a number without using a fancy intrinsic */
GCC_ALWAYS_INLINE unsigned fake_pc(uint64_t pc) {
  pc = (pc & UINT64_C(0x5555555555555555)) + ((pc & UINT64_C(0xAAAAAAAAAAAAAAAA)) >>  1);
  pc = (pc & UINT64_C(0x3333333333333333)) + ((pc & UINT64_C(0xCCCCCCCCCCCCCCCC)) >>  2);
  pc = (pc & UINT64_C(0x0F0F0F0F0F0F0F0F)) + ((pc & UINT64_C(0xF0F0F0F0F0F0F0F0)) >>  4);
  pc = (pc & UINT64_C(0x00FF00FF00FF00FF)) + ((pc & UINT64_C(0xFF00FF00FF00FF00)) >>  8);
  pc = (pc & UINT64_C(0x0000FFFF0000FFFF)) + ((pc & UINT64_C(0xFFFF0000FFFF0000)) >> 16);
  pc = (pc & UINT64_C(0x00000000FFFFFFFF)) + ((pc & UINT64_C(0xFFFFFFFF00000000)) >> 32);
  return (unsigned)pc; /* it will always be small enough to fit into an unsigned int (it won't be more than 63) */
}

GCC_ALWAYS_INLINE unsigned real_pc(uint64_t pc) {
  return (unsigned)(__builtin_popcount((unsigned)(pc & 0xFFFFFFFF))
                  + __builtin_popcount((unsigned)(pc >> 32)));
}

#ifndef POPCOUNT
#define POPCOUNT fake_pc
#endif

#define PRIME_LONG(v, _st) ((_st)[(v) >> 6])
#define PRIME_SHIFT(v)     ((v) & 63)
#define PRIME_NUM(v, _st)  ((PRIME_LONG(v, _st) >> PRIME_SHIFT(v)) & 1)

unsigned long calc_pi(unsigned long in) {
  uint64_t *storage;
  unsigned long anum = in-1;
  unsigned long nlongs = (anum + 63) >> 6;
  unsigned long sum = 0;

  storage = calloc(nlongs, sizeof(uint64_t));
  if (!storage) {
    fprintf(stderr, "Could not allocate space for %zu longs! (about %zu bytes): %s\n",
            nlongs, nlongs * sizeof(uint64_t), strerror(errno));
    abort();
  }

  /* fill storage with 1 bits */
  for (size_t i = 0; i < nlongs - 1; ++i) {
    storage[i] = ~UINT64_C(0);
  }

  storage[nlongs-1] = (UINT64_C(1) << (anum & 63)) - 1;

  /* the sieve */
  for (uint64_t n = 2; n * n <= in; ++n) {
    if (!PRIME_NUM(n - 2, storage)) continue;

    for (uint64_t m = n * 2; m <= in; m += n) {
      PRIME_LONG(m - 2, storage) &= ~(UINT64_C(1) << PRIME_SHIFT(m - 2));
    }
  }

  /* count up the primes */
  for (size_t n = 0; n < nlongs; ++n) {
    sum += POPCOUNT(storage[n]);
  }

  free(storage);

  return sum;
}

unsigned long get_num(const char *str) {
  char *end;
  unsigned long ret = strtoul(str, &end, 10);

  if (!*str || *end) {
    fputs("Invalid number entered :(\n", stderr);
    exit(EXIT_FAILURE);
  }

  if (ret < 2) {
    fputs("Number must not be less than 2\n", stderr);
    exit(EXIT_FAILURE);
  }

  return ret;
}

int main(int argc, char **argv) {
  if (argc < 2) {
    fprintf(stderr, "Usage: %s <num>\nWill print out the number of primes less than or equal to the supplied number.\n", argv[0] ? argv[0] : "(null)");
    return EXIT_FAILURE;
  }

  unsigned long num = get_num(argv[1]);
  unsigned long sum = calc_pi(num);

  double v = (sum * log(num)) / num;

  printf("%lu primes\n", sum);
  printf("%lf num\n", v);
  
  return EXIT_SUCCESS;
}