Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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 usage();
- inline bool prime_check(mpz_t n);
- void primepal_cand_inc(mpz_t n);
- void undulating_cand_inc(mpz_t n);
- void long_period_cand_inc(mpz_t n);
- int interesting(mpz_t n);
- uint64_t get_time();
- void prime_check_s(char* s);
- void undularray(int undula_arr[100]);
- int lazy_hop(int lazy);
- bool check_supp(int ab, int n);
- void divisor_count(mpz_t n);
- bool mini_chk(mpz_t n);
- int long_period_composite_chk(int first_two, int as);
- void mod_chk(int fac, int len, int as, mpz_t n);
- 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 > 3) {
- printf("unrecognized increment mode, valid modes are 0,1,2,3\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 == 3)) {
- // 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 - abbbabbba
- cand_increment = long_period_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 alternate period 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 long_period_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 = long_period_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 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 long_period_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:
- ret = 6815;
- break;
- case 6815:
- ret = 6959;
- break;
- case 6959:
- 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:
- ret = 12849;
- break;
- case 12849:
- ret = 13221;
- break;
- case 13221:
- ret = 13265;
- break;
- case 13265:
- 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 - large periodic increment mode, currently only three digits "
- "(122212221, etc.)\n"
- " 3 - 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"
- " 4 - impatient periodic mode - same as 2, 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");
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement