primepal - a C program to find probable palindromic primes, with many features for searching smoothly undulating palindromic primes, primes with limited numbers of digits, etc.

primepal.c

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <getopt.h>
#include <sys/time.h>
#include <time.h>
#include <gmp.h>

/*

  primepal.c - a c program to find probable palindromic primes

  this requires the GMP bignum lib
  on Mac OS X/clang, depending on gmp location something like this should work:
cc -Wall -L /usr/local/lib -I /usr/local/include -lgmp -o primepal primepal.c -O2

  on Linux/gcc something like this should work:
cc -std=c99 -Wall -lgmp -o primepal primepal.c -O2

  So far I've only compiled on Mac OS X, and a CentOS system.

  There are a few peculiar things going on in here.

  mpz_sizeinbase(n, 10) returns a base 10 len or len+1.  Since all integers
  we're dealing with are odd, we check for even and decrement and dance around
  to deal with this situation.
  The only base10 palindromic prime with an even # of digits is 11.

  the fastest way to check primality of a lot of large numbers is not to factor
  if we can find something's composite cheaply enough, so there's various checks
  mpz_probab_prime_p() does a fast check up to %53 before probabilistic checks
  when there are fast ways to spot %>53 cheaply, these pay off more

  there are a few lurking potential int overflows for really large numbers
  they're not the kind of things we're testing since they're so large they'd
  take very, very long runs

*/

void primepal_cand_inc(mpz_t n);
int interesting(mpz_t n);
void undulating_cand_inc(mpz_t n);
void undularray(int undula_arr[100]);
void second_degree_supp_cand_inc(mpz_t n);
int second_degree_composite_chk(int first_two, int as);
void third_degree_supp_cand_inc(mpz_t n);
int third_degree_composite_chk(int first_two, int as);
void repunit_cand_inc(mpz_t n);
inline bool prime_check(mpz_t n);
void mod_chk(int fac, int len, int as, mpz_t n);
uint64_t get_time();
void prime_check_s(char* s);
int lazy_hop(int lazy);
bool check_supp(int ab, int n);
void divisor_count(mpz_t n);
bool mini_chk(mpz_t n);
void usage();

uint64_t end_time;
uint64_t start_time;
uint64_t candidates = 0;
uint64_t primes = 0;
uint64_t factor_time;

// debug things, globals since we use all over the place
uint64_t mod_cnt[999] = {0};
// 164 elements
const int divisors[] = {
    1,   3,   7,   11,  13,  17,  19,  23,  29,  31,  31,  43,  47,  53,  59,
    61,  67,  71,  73,  79,  83,  97,  101, 103, 107, 109, 113, 127, 131, 137,
    139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
    227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
    311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
    401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487,
    491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
    599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
    683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
    797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
    887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997};

bool lazy = false;
bool verbose;
bool debug = false;
uint_fast8_t uniques;
uint_fast32_t max_len;
uint_fast8_t iterations = 5;

int main(int argc, char* argv[]) {
  mpz_t n;
  mpz_init(n);
  verbose = false;
  uniques = 0;
  max_len = 0;
  int begin_len = 0;
  int increment_mode = 0;
  int init_digit = 0;
  int ab = 0;
  int nn = 0;
  char* start;
  start = (char*)malloc(1);
  start[0] = '\0';

  int c;
  char* p;

  // parse args
  while ((c = getopt(argc, argv, "a:b:n:p:l:du:s:r:m:vf:h")) != -1) {
    switch (c) {
      case 'p':
        prime_check_s(optarg);
        exit(0);
        break;
      case 'v':
        verbose = true;
        break;
      case 'd':
        verbose = true;
        debug = true;
        break;
      case 's':
        start = (char*)realloc(start, strlen(optarg) + 1);
        strcpy(start, optarg);
        printf("Setting start to %s\n", start);
        max_len = strlen(start);
        break;
      case 'm':
        increment_mode = strtoul(optarg, &p, 10);
        if (increment_mode == 0 || increment_mode > 4) {
          printf("unrecognized increment mode, valid modes are 0,1,2,3,4\n");
          usage();
          exit(1);
        }
        break;
      case 'b':
        begin_len = strtoul(optarg, &p, 10);
        if (begin_len == 0) {
          printf("beginning length must be greater than 0\n");
          usage();
          exit(1);
        }
        if (begin_len % 2 == 0) {
          // must be odd
          begin_len++;
        }
        break;
      case 'l':
        max_len = strtoul(optarg, &p, 10);
        if (max_len < 3) max_len = 3;
        if (verbose) printf("Setting length to %d\n", max_len);
        break;
      case 'a':
        ab = strtoul(optarg, &p, 10);
        if ((ab < 10) || (ab > 99)) {
          printf("a must be two digits\n");
          exit(1);
        }
        if (verbose) printf("Setting a to %d\n", ab);
        break;
      case 'n':
        nn = strtoul(optarg, &p, 10);
        if (nn < 1) {
          printf("1 must > 0\n");
          exit(1);
        }
        if (verbose) printf("Setting n to %d\n", nn);
        break;
      case 'f':
        init_digit = (int)optarg[0] - 48;
        if ((init_digit < 1) || (init_digit > 9)) {
          printf("invalid start digit: %s", optarg);
          exit(1);
        }
        break;
      case 'u':
        uniques = labs(strtol(optarg, &p, 10));
        if (uniques == 0) {
          printf("uniques must be greater than 0\n");
          usage();
          exit(1);
        }
        if (verbose) printf("Setting uniques to %d\n", uniques);
        break;
      case 'r':
        iterations = labs(strtol(optarg, &p, 10));
        if (iterations < 5) {
          printf(
              "Setting number of repetitions to %d increases chances of "
              "finding pseudoprimes as length grows.\n",
              iterations);
        }
        if (verbose)
          printf(
              "Setting number of repetitions of probabilistic primality check "
              "to %d\n",
              iterations);
        break;
      case 'h':
        usage();
        exit(0);
    }
  }

  if ((ab > 0) && (nn > 0)) {
    // check (ab*10^len-ba)/99
    int tmp = ab;
    int a = tmp % 10;
    tmp = tmp / 10;
    int b = tmp % 10;
    printf("checking (%d*10^%d-%d%d)/99\n", ab, nn, a, b);
    bool prime = check_supp(ab, nn);
    if (prime) {
      printf("probable prime\n");
    } else {
      printf("not prime\n");
    }
    exit(0);
  }
  if (max_len == 0) {
    usage();
    exit(1);
  } else if (max_len % 2 == 0) {
    max_len++;
    if (verbose) printf("even start length - bumping to %d\n", max_len);
  }
  if (begin_len > max_len) {
    max_len = begin_len;
    printf(
        "beginning len set greater than end len, setting beginning len to %d\n",
        begin_len);
  }
  // set function pointer depending on mode and initial value
  void (*cand_increment)(mpz_t n);
  if (increment_mode == 0) {
    // simple increment
    if (begin_len > 0) {
      printf("begin length ignored\n");
    }
    int start_len = strlen(start);
    cand_increment = primepal_cand_inc;
    if (start_len > 0) {
      // we got a starting number, check and set
      int err = mpz_set_str(n, start, 10);
      if (err) {
        printf("format error with start: %s\n", start);
        exit(1);
      }
      if (start_len % 2 == 0) {
        printf("start digit's even - use integer with odd number of digits.\n");
        exit(1);
      }
    } else {
      // we got a length to start with, set start
      if (init_digit > 0) {
        mpz_ui_pow_ui(n, 10, max_len);
        mpz_sub_ui(n, n, 1);
        mpz_fdiv_q_ui(n, n, 9);
        mpz_mul_ui(n, n, init_digit);
      } else {
        mpz_ui_pow_ui(n, 10, max_len - 1);
        mpz_add_ui(n, n, 1);
      }
    }
  } else if ((increment_mode == 1) || (increment_mode == 5)) {
    // SUPP - abababa
    cand_increment = undulating_cand_inc;
    // we got a length to start with, set start
    if (init_digit > 0) {
      printf("initial digit ignored\n");
    }
    if (begin_len > 0) {
      mpz_ui_pow_ui(n, 10, begin_len);
      mpz_sub_ui(n, n, 1);
      mpz_fdiv_q_ui(n, n, 9);
    } else {
      int err = mpz_set_str(n, "101", 10);
      assert(err == 0);
    }
    if (increment_mode == 1) {
      if (verbose) printf("in SUPP mode\n");
    }
    if (increment_mode == 3) {
      lazy = true;
      if (verbose) printf("in lazy SUPP mode - skipping some lengths\n");
    }
  } else if (increment_mode == 2) {
    // longer period - abbabba
    cand_increment = second_degree_supp_cand_inc;
    if (begin_len > 0) {
      mpz_ui_pow_ui(n, 10, begin_len);
      mpz_sub_ui(n, n, 1);
      mpz_fdiv_q_ui(n, n, 9);
    } else {
      int err = mpz_set_str(n, "1221", 10);
      assert(err == 0);
    }
    if (verbose) printf("in second degree period mode\n");
  } else if (increment_mode == 3) {
    // longer period - abbbabbba
    cand_increment = third_degree_supp_cand_inc;
    if (begin_len > 0) {
      mpz_ui_pow_ui(n, 10, begin_len);
      mpz_sub_ui(n, n, 1);
      mpz_fdiv_q_ui(n, n, 9);
    } else {
      int err = mpz_set_str(n, "16661", 10);
      assert(err == 0);
    }
    if (verbose) printf("in third degree period mode\n");

  } else if (increment_mode == 4) {
    // repunit
    cand_increment = repunit_cand_inc;
    if (begin_len > 0) {
      mpz_ui_pow_ui(n, 10, begin_len);
      mpz_sub_ui(n, n, 1);
      mpz_fdiv_q_ui(n, n, 9);
    } else {
      int err = mpz_set_str(n, "1", 10);
      assert(err == 0);
    }
    if (verbose) printf("in repunit mode\n");
  } else {
    printf("error in getting mode - use 1,2,3\n");
    usage();
    exit(1);
  }

  free(start);

  if (verbose) {
    gmp_printf("starting with: %Zd\n\n", n);
  }
  start_time = get_time();
  factor_time = 0;
  // kick off function to do work, return when done.
  (*cand_increment)(n);

  if (verbose) {
    printf("prime candidates tested: %lld\n", candidates);
    printf("primes found: %lld\n", primes);
    printf("ratio: %f %%\n", (((double)primes) / ((double)candidates)) * 100.0);
  }
  uint64_t end_time = get_time();
  float duration = end_time - start_time;
  float duration_sec = (float)duration / 1000;
  float factor_sec = (float)factor_time / 1000;
  printf("factoring time: %fs\n", factor_sec);
  printf("elapsed time: %fs\n", duration_sec);

  if (debug) {
    // we keep track of some composite factors
    // useful for finding candidates to try to skip
    for (int i = 1; i < 999; i += 2) {
      if (mod_cnt[i] > 0) {
        printf("%d count: %llu\n", i, mod_cnt[i]);
      }
    }
  }
}

void primepal_cand_inc(mpz_t n) {
  // for smaller numbers we spend a lot more time building these than factoring
  // precalculate to shave time

  // skip when first digit is 2,4,5,6,8
  static const uint_fast8_t hop[10] = {0, 0, 1, 0, 1, 1, 1, 0, 1, 0};

  mpz_t tmp, ne10, ne11;
  mpz_init(tmp);
  mpz_init(ne10);
  mpz_init(ne11);

  // calc the things we only need to calc once here before we loop
  // sometimes mpz_sizeinbase reports len+1
  uint_fast32_t len = mpz_sizeinbase(n, 10);
  if (len % 2 == 0) {
    len--;
  }
  const uint_fast32_t mid = len / 2;
  const uint_fast8_t uniqs = uniques;

  // ne10 = 10^len-1
  mpz_ui_pow_ui(ne10, 10, len - 1);
  mpz_ui_pow_ui(ne11, 10, len - 1);
  mpz_add_ui(ne11, ne11, 1);

  // precalc n^10,100,1000, etc.
  mpz_t n_powers[mid];
  for (int i = 0; i <= mid; i++) {
    mpz_init(n_powers[i]);
    mpz_ui_pow_ui(n_powers[i], 10, i);
  }
  while (1) {
    /*
     increment based on context:

     if checking uniques, grab first half of digits and increment until we get
     to the next candidate that has the correct number of unique digits, covert
     to string, reverse string, add reversed back in.

     if straight iteration, check if mid digit < 9, if so add 10^len/2 to
     increment, else convert to string and reverse.

     if first digit is 2,4,5,6,8 skip range, since can't be prime.

     once we're past max_len we return
    */
    if (uniqs) {
      // cut len in half, increment as much as needed, then make string
      mpz_tdiv_q(n, n, n_powers[mid]);
      mpz_add_ui(n, n, 1);
      int_fast8_t i;
      while ((i = interesting(n)) != -1) {
        mpz_add(n, n, n_powers[i]);
      }
    } else {
      // 90% of time skip stringy slowness
      // get first digit (basically n0,000 / 10,000)
      mpz_fdiv_q(tmp, n, ne10);
      const uint_fast8_t first = mpz_get_ui(tmp);
      if (hop[first]) {
        // if first digit is 2,4,5,6,8 can't be prime, hop to n0...0n
        mpz_mul_ui(n, ne11, first + 1);
      }
      // get last digit (n/(10^len/2))%10
      mpz_tdiv_q(tmp, n, n_powers[mid]);
      const uint_fast8_t last = mpz_mod_ui(tmp, tmp, 10);
      if (last != 9) {
        // we can just increment the mid-digit n+=(10^len/2) and we're done
        mpz_add(n, n, n_powers[mid]);
        if (mpz_gcd_ui(tmp, n, 22141119) <= 1) {
          prime_check(n);
        }
        continue;
      } else {
        // have to make string - cut len in half and inc by one
        mpz_tdiv_q(n, n, n_powers[mid]);
        mpz_add_ui(n, n, 1);
      }
    }

    // grab first half in str, reverse
    char last_half[mid + 1];
    mpz_get_str(last_half, 10, n);
    // check len of incremented cand.
    // if we've bumped len, last_half[mid + 1] will be a num.
    if (last_half[mid + 1] != '\0') {
      return;
    }
    // check first digit
    const uint_fast8_t first = last_half[0] - 48;
    if (hop[first]) {
      // if first digit is 2,4,5,6,8 can't be prime, hop to n0...0n
      // n = (10^len)*(first+1)+(first+1)
      mpz_mul_ui(n, ne11, first + 1);
    } else {
      // chop last digit and reverse
      last_half[mid] = '\0';
      char *c1, *c2;
      for (c1 = last_half, c2 = last_half + mid - 1; c2 > c1; ++c1, --c2) {
        *c1 ^= *c2;
        *c2 ^= *c1;
        *c1 ^= *c2;
      }
      // last_half from str to bignum
      assert(mpz_set_str(tmp, last_half, 10) == 0);
      // multiply n by 10^mid to get n to full length with latter half zeroed
      // then add first half to reversed half to get next palindrome
      mpz_mul(n, n, n_powers[mid]);
      mpz_add(n, n, tmp);
    }

    // factors that have a higher freq, check if divisible before proceeding
    // to full check
    // best so far:
    // 3*7*11*13*73*101 = 22141119
    if (mpz_gcd_ui(tmp, n, 22141119) <= 1) {
      prime_check(n);
    }
  }
}

int interesting(mpz_t n) {
  // filter for x unique digits in number
  // returns int
  // -1=interesting
  // else return power of ten for next increment

  // convert to str to scan digits
  uint_fast32_t len = (max_len / 2) + 1;
  const uint_fast8_t uniqs = uniques;
  char str[len];
  mpz_get_str(str, 10, n);
  // uint_fast8_t len = strlen(str);
  uint_fast8_t found = 0;
  bool count[10] = {false};
  // walk digits, tracking freq.
  // once we hit > uniqs return place to allow faster increment
  for (int i = 0; i <= len; i++) {
    if (count[str[i] - 48] == false) {
      found++;
      if (found > uniqs) {
        // we've hit the place to increment
        return len - i - 1;
      }
      count[str[i] - 48] = true;
    }
  }
  return -1;
}

void undulating_cand_inc(mpz_t n) {
  // build SUPP xyxyxyx candidates (also 11...11 repunits)
  mpz_t tmp;
  mpz_init(tmp);
  int undula_arr[100];
  undularray(undula_arr);
  bool done = false;
  while (done == false) {
    int len = mpz_sizeinbase(n, 10);
    // get len, if even decrement
    if (len % 2 == 0) {
      len--;
    }
    if (len > max_len) {
      done = true;
      mpz_clear(tmp);
      return;
    }

    // get first two digits
    mpz_ui_pow_ui(tmp, 10, len - 2);
    mpz_tdiv_q(tmp, n, tmp);
    int first_two = mpz_get_ui(tmp);

    // used to check if candidate is %3
    int firsts = len / 2 + 1;
    int seconds = firsts - 1;
    int first, second;

    // increment - walk undula_arr to skip composites
    bool found = false;
    while (found == false) {
      if (first_two == 98) {
        // this gets more complex due to laziness
        if (lazy == true) {
          int next_len = lazy_hop(len);
          if (next_len > 0) {
            len = next_len;
          } else if (next_len == -1) {
            // past the point we have values, shift to slow mode
            printf("\n\n");
            lazy = false;
            len = len + 2;
            firsts++;
            seconds++;
          } else {
            len = len + 2;
            firsts++;
            seconds++;
          }
        } else {
          len = len + 2;
          firsts++;
          seconds++;
        }
        first_two = 10;
      }
      first_two = undula_arr[first_two];

      // hop some composite numbers that are periodic
      if (first_two == 11) {
        // these are %41==0
        if (len % 5 == 0) first_two = 12;
      }
      if (first_two == 12) {
        if (((len + 25) % 30) == 0) first_two = 13;
      } else if (first_two == 14) {
        if (((len + 7) % 30) == 0) first_two = 16;
      } else if (first_two == 16) {
        if ((((len) / 3) % 2) == 1) first_two = 17;
      } else if (first_two == 18) {
        //% 31
        if (((len + 23) % 30) == 0) first_two = 19;
      } else if (first_two == 34) {
        // never prime if len-3%3==0
        //(34*10n-43)/99 - only primes
        if (len % 3 == 0) first_two = 35;
        if (len % 3 == 1) first_two = 35;
      } else if (first_two == 35) {
        // these are %211==0
        if ((len - 9) % 10 == 0) first_two = 37;
      } else if (first_two == 37) {
        if ((((len - 2) / 3) % 2) == 1) first_two = 38;
      } else if (first_two == 71) {
        // biggest hop here
        if ((len - 1) % 3 == 0) {
          first_two = 91;
        }
      } else if (first_two == 75) {
        //% 31
        if (((len + 9) % 30) == 0) first_two = 78;
      } else if (first_two == 78) {
        //% 31
        if (((len + 19) % 30) == 0) first_two = 79;
      } else if (first_two == 92) {
        //% 31
        if (((len + 9) % 30) == 0) first_two = 94;
      } else if (first_two == 94) {
        // 94%3==0 is %13
        if (len % 3 == 0) first_two = 95;
      } else if (first_two == 95) {
        //% 31
        if (((len + 7) % 30) == 0) first_two = 97;
      }

      // check mod3, if clear then check mod11 based on digits before
      // building gmp bignum
      first = first_two;
      second = first % 10;
      first /= 10;
      uint64_t digit_sum = firsts * first + seconds * second;
      if (digit_sum % 3 != 0) {
        uint64_t digit_diff = abs(firsts * first - seconds * second);
        if (digit_diff % 11 != 0) {
          found = true;
        }
      }
    }

    int ab = first * 10 + second;
    int ba = second * 10 + first;

    //(ab*10^len-ba)/99 = aba...aba
    mpz_ui_pow_ui(tmp, 10, len);
    mpz_mul_ui(n, tmp, ab);
    mpz_sub_ui(n, n, ba);
    mpz_tdiv_q_ui(n, n, 99);

    // skip lower digits, mpz_gcd_ui deals with <49
    // 61*71*79*97*239=7932040267
    // 61*71*73*79*97*239=579038939491
    // factors that have a higher freq, check if divisible before proceeding to
    // check
    // probably could use some tuning
    if (mpz_gcd_ui(tmp, n, 579038939491) <= 1) {
      if (prime_check(n) && verbose) {
        printf("(%d*10^%d-%d)/99\n\n", ab, len, ba);
      }
    }
  }
  // done
}

void second_degree_supp_cand_inc(mpz_t n) {
  mpz_t tmp;
  mpz_init(tmp);
  bool done = false;
  while (done == false) {
    int len = mpz_sizeinbase(n, 10);
    // get len, if off by one, decrement
    if ((len - 1) % 3 != 0) {
      len--;
    }

    if (len > max_len) {
      done = true;
      mpz_clear(tmp);
      return;
    }

    // get first two digits
    mpz_ui_pow_ui(tmp, 10, len - 2);
    mpz_tdiv_q(tmp, n, tmp);
    int first_two = mpz_get_ui(tmp);
    int first = first_two / 10;
    int second = first_two % 10;

    // given abbabba...abbabba, as,bs = instances of a and b
    int as = ((len - 1) / 3) + 1;
    int bs = len - as;

    bool found = false;
    while (found == false) {
      if (first_two == 19) {
        first_two = 31;
      } else if (first_two == 39) {
        // 40 = hop
        len = len + 3;
        as = as + 1;
        bs = bs + 2;
        first_two = 12;
        // even len = divisible by 11, skip.
        if (len % 2 == 0) {
          len += 3;
          as = as + 1;
          bs = bs + 2;
        }
      } else {
        first_two++;
      }

      first_two = second_degree_composite_chk(first_two, as);
      // get new first,second
      second = first_two % 10;
      first = first_two / 10;

      // check mod 3 without gmp, if clear then check mod 11 without gmp
      // if both clear, then build the bignum and check primality
      uint64_t digit_sum = as * first + bs * second;
      if (digit_sum % 3 != 0) {
        found = true;
      }
    }
    //(abb*10^n-bba)/999 = abbabba...abbabba
    int abb = first * 100 + second * 10 + second;
    int bba = second * 100 + second * 10 + first;
    mpz_ui_pow_ui(tmp, 10, len);
    mpz_mul_ui(n, tmp, abb);
    mpz_sub_ui(n, n, bba);
    mpz_tdiv_q_ui(n, n, 999);

    // skip lower digits, mpz_gcd_ui deals with <49
    // 7,11,31=2387
    // factors that have a higher freq, check if divisible before proceeding to
    // check
    // probably could use some tuning

    // mod_chk(17, len, as, n);
    if (mpz_gcd_ui(tmp, n, 77) <= 1) {
      if (prime_check(n) && verbose) {
        printf("(%d*10^%d-%d)/999\n\n", abb, len, bba);
      }
    }
  }
  // done
}

void third_degree_supp_cand_inc(mpz_t n) {
  mpz_t tmp;
  mpz_init(tmp);
  bool done = false;
  while (done == false) {
    int len = mpz_sizeinbase(n, 10);
    // get len, if even decrement
    if (len % 2 == 0) {
      len--;
    }
    if (len > max_len) {
      done = true;
      mpz_clear(tmp);
      return;
    }

    // get first two digits
    mpz_ui_pow_ui(tmp, 10, len - 2);
    mpz_tdiv_q(tmp, n, tmp);
    int first_two = mpz_get_ui(tmp);
    int first = first_two / 10;
    int second = first_two % 10;

    // given abbbabbba...abbbabbba, as,bs = instances of a and b
    int as = ((len - 1) / 4) + 1;
    int bs = len - as;

    bool found = false;
    while (found == false) {
      if (first_two == 19) {
        // 19..70==prime wasteland
        first_two = 71;
      } else {
        first_two++;
      }

      if (first_two == 79) {
        len = len + 4;
        as = as + 1;
        bs = bs + 3;
        first_two = 11;
      }
      first_two = third_degree_composite_chk(first_two, as);
      // get new first,second
      second = first_two % 10;
      first = first_two / 10;

      // check mod 3 without gmp, if clear then check mod 11 without gmp
      // if both clear, then build the bignum and check primality
      uint64_t digit_sum = as * first + bs * second;
      if (digit_sum % 3 != 0) {
        // check mod 11 - add up odd digits & even digits
        // subtract evens from odds and check for divisibility by 11
        int odds = (first * as) + (as - 1) * second;
        int evens = (len / 2) * second;
        if (abs(odds - evens) % 11 != 0) {
          found = true;
        }
      }
    }
    //(abbb*10^n-bbba)/9999 = abbbabbba...abbbabbba
    int abbb = first * 1000 + second * 100 + second * 10 + second;
    int bbba = second * 1000 + second * 100 + second * 10 + first;
    mpz_ui_pow_ui(tmp, 10, len);
    mpz_mul_ui(n, tmp, abbb);
    mpz_sub_ui(n, n, bbba);
    mpz_tdiv_q_ui(n, n, 9999);

    // 3*7*11*61*101=1423191
    // 3*7*11*19*61*101=27040629
    // 3*7*11*61*101*449=639012759
    // 7*61*101*239*449=4628001497

    // skip lower digits, mpz_gcd_ui deals with <49
    // 61*101*239*449=661143071
    // factors that have a higher freq, check if divisible before proceeding to
    // check
    // probably could use some tuning
    if (mpz_gcd_ui(tmp, n, 661143071) <= 1) {
      if (prime_check(n) && verbose) {
        printf("(%d*10^%d-%d)/9999\n\n", abbb, len, bbba);
      }
    }
  }
  // done
}

void repunit_cand_inc(mpz_t n) {
  // build repunit cands
  // length is always prime, so hop along primes for lengths to check
  // they are 2, 19, 23, 317, 1031, 49081, 86453, 109297, 270343...
  // so after 1031 we're grinding a really long time.
  // use -b to start somewhere interesting

  bool done = false;
  mpz_t tmp;
  mpz_init(tmp);
  int len = mpz_sizeinbase(n, 10);
  // get len, if even decrement
  if (len % 2 == 0) {
    len--;
  }
  while (done == false) {
    mpz_set_ui(tmp, len);
    bool next_len_found = false;
    // mpz_nextprime() skips some primes, so increment and test
    while (next_len_found == false) {
      uint_fast8_t p = mpz_probab_prime_p(tmp, 1);
      if (p > 0) {
        next_len_found = true;
      } else {
        mpz_add_ui(tmp, tmp, 1);
      }
    }
    len = mpz_get_ui(tmp);
    if (len > max_len) {
      done = true;
      mpz_clear(tmp);
      return;
    }
    //((10^len)-1)/9 = aba...aba
    mpz_ui_pow_ui(n, 10, len);
    mpz_sub_ui(n, n, 1);
    mpz_fdiv_q_ui(n, n, 9);

    if (prime_check(n) && verbose) {
      printf("((10^%d)-1)/9\n\n", len);
    }
    // len is odd so we can inc by 2 and skip 3s
    len += 2;
    if (len % 3 == 0) {
      len += 2;
    }
  }
  // done
}

void mod_chk(int fac, int len, int as, mpz_t n) {
  // log whether fac is a factor
  // useful for finding composites in undulating candidates to skip
  if (mpz_divisible_ui_p(n, fac) > 0) {
    gmp_printf("mod %d - %Zd - %d \nas - %d\n", fac, n, len, as);
  }
}

int second_degree_composite_chk(int first_two, int as) {
  // TODO
  if (first_two == 14) {
    //% 31
    if ((as + 1) % 10 == 0) {
      first_two = 15;
    }
  } else if (first_two == 19) {
    //% 19
    if ((as + 1) % 6 == 0) {
      first_two = 31;
    }
  }
  if (first_two == 31) {
    //% 31
    if ((as + 7) % 10 == 0) {
      first_two = 32;
    }
  }
  if (first_two == 38) {
    //% 19
    if ((as + 1) % 6 == 0) {
      first_two = 39;
    }
  }
  return first_two;
}

int third_degree_composite_chk(int first_two, int as) {
  /*
    Should have documented these more as they filled in...
    there are various patterns that show up with some sets of
    first_two at various lengths that are always composite.
    -- based on tests, not proofs
    -- there may be cases where this skips primes
  */

  if (first_two == 10) {
    // don't have proof, but these have been composite up to 20,000 digits
    first_two = 11;
  }
  if (first_two == 11) {
    if ((as + 15) % 13 == 0) {
      //%19==0
      first_two = 12;
    } else if ((as - 17) % 13 == 0) {
      //%19==0
      first_two = 12;
    } else if ((as + 22) % 13 == 0) {
      //%53==0
      first_two = 12;
    } else if ((as - 2) % 15 == 0) {
      //%71==0
      first_two = 12;
    }
  }
  if (first_two == 12) {
    if ((as + 17) % 21 == 0) {
      first_two = 13;
    } else if ((as + 10) % 7 == 0) {
      first_two = 13;
    } else if ((as + 6) % 13 == 0) {
      first_two = 13;
    }
  }
  if (first_two == 13) {
    if ((as + 65) % 69 == 0) {
      first_two = 14;
    } else if ((as + 15) % 23 == 0) {
      first_two = 14;
    } else if ((as + 4) % 9 == 0) {
      first_two = 14;
    } else if ((as + 13) % 17 == 0) {
      first_two = 14;
    } else if ((as + 21) % 53 == 0) {
      first_two = 14;
    }
  }
  if (first_two == 14) {
    if ((as + 38) % 69 == 0) {
      first_two = 15;
    } else if ((as + 65) % 69 == 0) {
      first_two = 15;
    } else if ((as + 7) % 15 == 0) {
      first_two = 15;
    } else if ((as + 16) % 13 == 0) {
      first_two = 15;
    } else if ((as + 1) % 21 == 0) {
      first_two = 15;
    }
  }
  if (first_two == 15) {
    if ((as + 6) % 11 == 0) {
      //%23==0
      first_two = 16;
    } else if ((as + 34) % 45 == 0) {
      first_two = 16;
    } else if ((as + 4) % 45 == 0) {
      first_two = 16;
    } else if ((as + 4) % 15 == 0) {
      //%31
      first_two = 16;
    } else if ((as + 56) % 69 == 0) {
      first_two = 16;
    }
  }
  if (first_two == 16) {
    if (as % 3 == 0) {
      //%19==0
      first_two = 17;
    } else if ((as + 1) % 9 == 0) {
      first_two = 17;
    } else if ((as + 3) % 23 == 0) {
      first_two = 17;
    } else if ((as + 5) % 21 == 0) {
      first_two = 17;
    } else if ((as + 24) % 29 == 0) {
      first_two = 17;
    } else if ((as + 49) % 53 == 0) {
      // mod107
      first_two = 17;
    }
  }
  if (first_two == 17) {
    if ((as + 4) % 11 == 0) {
      first_two = 18;
    } else if ((as + 1) % 3 == 0) {
      //%13
      first_two = 18;
    } else if ((as + 11) % 21 == 0) {
      first_two = 18;
    } else if ((as + 2) % 39 == 0) {
      first_two = 18;
    } else if ((as + 53) % 69 == 0) {
      first_two = 18;
    } else if ((as + 14) % 29 == 0) {
      first_two = 18;
    } else if ((as + 4) % 53 == 0) {
      //%107
      first_two = 18;
    }
  }
  if (first_two == 18) {
    if ((as - 3) % 4 == 0) {
      first_two = 19;
    } else if ((as + 8) % 21 == 0) {
      first_two = 19;
    } else if ((as + 1) % 21 == 0) {
      first_two = 19;
    } else if ((as + 24) % 13 == 0) {
      first_two = 19;
    } else if ((as + 43) % 69 == 0) {
      first_two = 19;
    } else if ((as + 20) % 69 == 0) {
      first_two = 19;
    }
  }
  if (first_two == 19) {
    if ((as + 5) % 11 == 0) {
      first_two = 70;
    } else if ((as + 1) % 13 == 0) {
      first_two = 70;
    } else if ((as + 10) % 21 == 0) {
      first_two = 70;
    } else if ((as + 8) % 23 == 0) {
      first_two = 70;
    } else if ((as + 11) % 29 == 0) {
      first_two = 70;
    } else if ((as + 4) % 5 == 0) {
      first_two = 70;
    }
  }
  if (first_two == 70) {
    first_two = 71;
  }
  if (first_two == 71) {
    if ((as + 7) % 9 == 0) {
      first_two = 72;
    } else if ((as + 7) % 13 == 0) {
      first_two = 72;
    } else if ((as + 2) % 23 == 0) {
      first_two = 72;
    } else if ((as + 83) % 105 == 0) {
      first_two = 72;
    } else if ((as + 38) % 53 == 0) {
      //%107
      first_two = 72;
    }
  }
  if (first_two == 72) {
    if ((as + 4) % 39 == 0) {
      first_two = 73;
    } else if ((as + 7) % 69 == 0) {
      first_two = 73;
    } else if ((as + 4) % 9 == 0) {
      first_two = 73;
    } else if ((as + 8) % 15 == 0) {
      //%31
      first_two = 73;
    } else if ((as + 88) % 105 == 0) {
      first_two = 73;
    } else if ((as + 15) % 53 == 0) {
      //%107
      first_two = 73;
    }
  }
  if (first_two == 73) {
    if ((as + 4) % 11 == 0) {
      first_two = 74;
    } else if ((as + 18) % 23 == 0) {
      first_two = 74;
    } else if ((as + 4) % 17 == 0) {
      first_two = 74;
    } else if ((as + 27) % 29 == 0) {
      first_two = 74;
    }
  }
  if (first_two == 74) {
    if (as % 3 == 0) {
      //%19==0
      first_two = 75;
    } else if ((as + 2) % 11 == 0) {
      first_two = 75;
    } else if ((as + 4) % 7 == 0) {
      first_two = 75;
    } else if ((as + 16) % 7 == 0) {
      first_two = 75;
    } else if ((as + 4) % 15 == 0) {
      first_two = 75;
    } else if ((as + 1) % 9 == 0) {
      first_two = 75;
    } else if ((as + 1) % 23 == 0) {
      first_two = 75;
    } else if ((as + 11) % 29 == 0) {
      first_two = 75;
    } else if ((as + 17) % 21 == 0) {
      first_two = 75;
    }
  }
  if (first_two == 75) {
    if ((as + 1) % 4 == 0) {
      //%19==0
      first_two = 76;
    } else if ((as + 2) % 9 == 0) {
      //%19==0
      first_two = 76;
    } else if ((as + 25) % 29 == 0) {
      //%19==0
      first_two = 76;
    } else if ((as + 10) % 15 == 0) {
      // mod 71, needs review
      first_two = 76;
    } else if ((as + 5) % 15 == 0) {
      // mod 71, needs review
      first_two = 76;
    }
  }
  if (first_two == 76) {
    if ((as + 1) % 15 == 0) {
      //%19==0
      first_two = 78;
    } else if ((as + 1) % 21 == 0) {
      first_two = 78;
    } else if ((as + 11) % 33 == 0) {
      //%67
      first_two = 78;
    } else if ((as + 7) % 41 == 0) {
      //%83
      first_two = 78;
    } else if ((as + 37) % 53 == 0) {
      //%107
      first_two = 78;
    }
  }
  if (first_two == 77) {
    first_two = 78;
  }
  // for 78 we just have to let it go
  return first_two;
}

bool prime_check(mpz_t n) {
  candidates++;
  uint64_t before_chk = get_time();
  uint_fast8_t p = mpz_probab_prime_p(n, iterations);
  factor_time += get_time() - before_chk;
  if (p > 0) {
    gmp_printf("%Zd\n", n);
    primes++;
    return true;
  } else if (debug) {
    divisor_count(n);
  }
  return false;
}

uint64_t get_time() {
  struct timeval tv;
  uint64_t ret;
  gettimeofday(&tv, NULL);
  // get microseconds, convert to milliseconds
  ret = tv.tv_usec;
  ret /= 1000;
  // get seconds, convert to milliseconds
  ret += (tv.tv_sec * 1000);
  return ret;
}

bool mini_chk(mpz_t n) {
  // check smaller primes, return true if we aren't divisible by them
  // return immediately if we find we're composite
  // this is similar to what mpz_probab_prime_p() does but we go a
  // little past %53

  // poor perf. once numbers are >1000 digits
  for (int i = 1; i < 27; i++) {
    if (mpz_divisible_ui_p(n, divisors[i]) > 0) {
      return false;
    }
  }
  return true;
}

void divisor_count(mpz_t n) {
  // count divisors to find any that are interesting
  for (int i = 1; i < 164; i++) {
    if (mpz_divisible_ui_p(n, divisors[i]) > 0) {
      mod_cnt[divisors[i]]++;
    }
  }
}

void undularray(int undula_arr[100]) {
  // it works, stop judging me

  // build array to loop for next candidate
  // skip matching digits if not 1
  // skip leading 20s,40s,50s,60s,80s
  // skip 76, always factorable by either 3, 7, or 13
  // skip 77, always factorable
  undula_arr[10] = 11;
  undula_arr[11] = 12;
  undula_arr[12] = 13;
  undula_arr[13] = 14;
  undula_arr[14] = 15;
  undula_arr[15] = 16;
  undula_arr[16] = 17;
  undula_arr[17] = 18;
  undula_arr[18] = 19;
  undula_arr[19] = 31;
  undula_arr[31] = 32;
  undula_arr[32] = 34;
  undula_arr[34] = 35;
  undula_arr[35] = 37;
  undula_arr[37] = 38;
  undula_arr[38] = 71;
  undula_arr[71] = 72;
  undula_arr[72] = 73;
  undula_arr[73] = 74;
  undula_arr[74] = 75;
  undula_arr[75] = 78;
  undula_arr[78] = 79;
  undula_arr[79] = 91;
  undula_arr[91] = 92;
  undula_arr[92] = 94;
  undula_arr[94] = 95;
  undula_arr[95] = 97;
  undula_arr[97] = 98;
  undula_arr[98] = 10;
}

int lazy_hop(int lazy) {
  // skip lengths with no SUPPs
  // yes, this is silly in many ways, but it helps to skip past the shorter
  // range to grind through to the bigger ones
  int ret = 0;
  switch (lazy) {
    case 627:
      ret = 849;
      break;
    case 849:
      ret = 857;
      break;
    case 857:  // 857 is prime
      ret = 961;
      break;
    case 961:
      ret = 1031;  // 1031 is prime
      break;
    case 1031:
      ret = 1079;
      break;
    case 1079:
      ret = 1295;
      break;
    case 1295:
      ret = 1399;
      break;
    case 1399:  // 1399 is prime
      ret = 1415;
      break;
    case 1415:
      ret = 1479;
      break;
    case 1479:
      ret = 1597;
      break;
    case 1597:  // 1597 is prime
      ret = 1631;
      break;
    case 1631:
      ret = 1791;
      break;
    case 1791:
      ret = 1893;
      break;
    case 1893:
      ret = 1969;
      break;
    case 1969:
      ret = 1979;
      break;
    case 1979:  // 1979 is prime
      ret = 2069;
      break;
    case 2069:  // 2069 is prime
      ret = 2075;
      break;
    case 2075:
      ret = 2165;
      break;
    case 2165:
      ret = 2177;
      break;
    case 2177:
      ret = 2335;
      break;
    case 2335:
      ret = 2529;
      break;
    case 2529:
      ret = 2883;
      break;
    case 2883:
      ret = 3015;
      break;
    case 3015:
      ret = 3047;
      break;
    case 3047:
      ret = 3057;
      break;
    case 3057:
      ret = 3147;
      break;
    case 3147:
      ret = 3407;
      break;
    case 3407:  // 3407 is prime
      ret = 3503;
      break;
    case 3503:
      ret = 3657;
      break;
    case 3657:
      ret = 4157;
      break;
    case 4157:  // 4157 is prime
      ret = 4233;
      break;
    case 4233:
      ret = 4335;
      break;
    case 4335:
      ret = 4573;
      break;
    case 4573:
      ret = 4859;
      break;
    case 4859:
      ret = 4885;
      break;
    case 4885:
      ret = 6249;
      break;
    case 6249:
      ret = 6343;
      break;
    case 6343:  // 6343 is prime
      ret = 6815;
      break;
    case 6815:
      ret = 6959;
      break;
    case 6959:  // 6959 is prime
      ret = 7233;
      break;
    case 7233:
      ret = 7387;
      break;
    case 7387:
      ret = 7809;
      break;
    case 7809:
      ret = 8315;
      break;
    case 8315:
      ret = 9273;
      break;
    case 9273:
      ret = 9591;
      break;
    case 9591:
      ret = 9599;
      break;
    case 9599:
      ret = 10419;
      break;
    case 10419:
      ret = 11393;
      break;
    case 11393:  // 11393 is prime
      ret = 11399;
      break;
    case 11399:  // prime is prime
      ret = 12849;
      break;
    case 12849:
      ret = 13221;
      break;
    case 13221:
      ret = 13265;
      break;
    case 13265:
      ret = 14059;
      break;
    case 14059:
      ret = 14337;
      break;
    case 14337:
      ret = 14579;
      break;
    case 14579:
      ret = 14689;
      break;
    case 14689:
      ret = 14691;
      break;
    case 14691:
      ret = 15293;
      break;
    case 15293:
      ret = -1;
      break;
  }
  return ret;
}

bool check_supp(int ab, int nn) {
  // just for -a, -b flags
  mpz_t n;
  mpz_init(n);
  mpz_t tmp;
  mpz_init(tmp);
  //(ab*10^len-ba)/99 = aba...aba
  mpz_ui_pow_ui(tmp, 10, nn);
  mpz_mul_ui(n, tmp, ab);
  mpz_sub_ui(n, n, ab);
  mpz_cdiv_q_ui(n, n, 99);
  if (verbose) {
    gmp_printf("%Zd\n", n);
  }
  int pp = mpz_probab_prime_p(n, 20);
  if (pp != 0) {
    return true;
  }
  return false;
}

void prime_check_s(char* s) {
  // just for -p flag
  mpz_t n;
  mpz_init(n);
  int err = mpz_set_str(n, s, 10);
  if (err) {
    printf("format error integer arg: %s\n", s);
    exit(1);
  }
  int pp = mpz_probab_prime_p(n, 25);
  if (pp == 0) {
    puts("not prime.\n");
  } else if (pp == 1) {
    puts("probably prime.\n");
  } else {
    puts("prime.\n");
  }
}

void usage() {
  printf(
      "usage: primepal (-p [number to test]) or (-s [starting number for "
      "search])\n"
      " -p [number] - check if individual number is prime\n"
      " -s [start range]\n"
      " [-m] 0 - (default) increment to next palindrome, if n, increment to "
      "next palindrome with that number of digits\n"
      "      1 - undulating mode (101, 18181, etc.), starting with 101, "
      "ending when max_len set with -l is hit\n"
      "      2 - second degree periodic mode, e.g. abbabbabba \n"
      "      3 - third degree periodic mode, e.g. abbbabbbabbba \n"
      "      4 - repunit mode, e.g. 11, 1111111111111111111 \n"
      "      5 - impatient undulating mode - same as 1, but skips grinding "
      "on some ranges in longer lengths in the 1000-10000 digit range where "
      "there are no primes to be found\n"
      " [-d] limit search to primes with n unique digits (mode 0) \n"
      " [-f] set all starting digits to num (mode 0)\n"
      " [-r] set repetitions of primality test (default 15)\n"
      " [-v] verbose\n"
      " [-h] show this usage message\n");
}