add code

1 files changed, 110 insertions, 0 deletions

diff --git a/calc_pi.c b/calc_pi.c
new file mode 100644
index 0000000..c3aa4da
--- /dev/null
+++ b/calc_pi.c
@@ -0,0 +1,110 @@
+#include <stdio.h>
+#include <stdint.h>
+#include <inttypes.h>
+#include <string.h>
+#include <errno.h>
+#include <stdlib.h>
+#include <math.h>
+
+#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;
+}