Advertisement
Guest User

Untitled

a guest
Feb 7th, 2016
49
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 33.73 KB | None | 0 0
  1. #include <assert.h>
  2. #include <math.h>
  3. #include <stdbool.h>
  4. #include <stdint.h>
  5. #include <string.h>
  6. #include <stdio.h>
  7. #include <stdlib.h>
  8. #include <unistd.h>
  9. #include <getopt.h>
  10. #include <sys/time.h>
  11. #include <time.h>
  12. #include <gmp.h>
  13.  
  14. /*
  15.  
  16. primepal.c - a c program to find probable palindromic primes
  17.  
  18. this requires the GMP bignum lib
  19. on Mac OS X/clang, depending on gmp location something like this should work:
  20. cc -Wall -L /usr/local/lib -I /usr/local/include -lgmp -o primepal primepal.c -O2
  21.  
  22. on Linux/gcc something like this should work:
  23. cc -std=c99 -Wall -lgmp -o primepal primepal.c -O2
  24.  
  25. So far I've only compiled on Mac OS X, and a CentOS system.
  26.  
  27. There are a few peculiar things going on in here.
  28.  
  29. mpz_sizeinbase(n, 10) returns a base 10 len or len+1. Since all integers
  30. we're dealing with are odd, we check for even and decrement and dance around
  31. to deal with this situation.
  32. The only base10 palindromic prime with an even # of digits is 11.
  33.  
  34. the fastest way to check primality of a lot of large numbers is not to factor
  35. if we can find something's composite cheaply enough, so there's various checks
  36. mpz_probab_prime_p() does a fast check up to %53 before probabilistic checks
  37. when there are fast ways to spot %>53 cheaply, these pay off more
  38.  
  39. there are a few lurking potential int overflows for really large numbers
  40. they're not the kind of things we're testing since they're so large they'd
  41. take very, very long runs
  42.  
  43. */
  44.  
  45. void usage();
  46. inline bool prime_check(mpz_t n);
  47. void primepal_cand_inc(mpz_t n);
  48. void undulating_cand_inc(mpz_t n);
  49. void long_period_cand_inc(mpz_t n);
  50. int interesting(mpz_t n);
  51. uint64_t get_time();
  52. void prime_check_s(char* s);
  53. void undularray(int undula_arr[100]);
  54. int lazy_hop(int lazy);
  55. bool check_supp(int ab, int n);
  56. void divisor_count(mpz_t n);
  57. bool mini_chk(mpz_t n);
  58. int long_period_composite_chk(int first_two, int as);
  59. void mod_chk(int fac, int len, int as, mpz_t n);
  60.  
  61. uint64_t end_time;
  62. uint64_t start_time;
  63. uint64_t candidates = 0;
  64. uint64_t primes = 0;
  65. uint64_t factor_time;
  66.  
  67.  
  68. // debug things, globals since we use all over the place
  69. uint64_t mod_cnt[999] = {0};
  70. // 164 elements
  71. const int divisors[] = {
  72. 1, 3, 7, 11, 13, 17, 19, 23, 29, 31, 31, 43, 47, 53, 59,
  73. 61, 67, 71, 73, 79, 83, 97, 101, 103, 107, 109, 113, 127, 131, 137,
  74. 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
  75. 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
  76. 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
  77. 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487,
  78. 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
  79. 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
  80. 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
  81. 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
  82. 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997};
  83.  
  84. bool lazy = false;
  85. bool verbose;
  86. bool debug = false;
  87. uint_fast8_t uniques;
  88. uint_fast32_t max_len;
  89. uint_fast8_t iterations = 5;
  90.  
  91. int main(int argc, char* argv[]) {
  92. mpz_t n;
  93. mpz_init(n);
  94. verbose = false;
  95. uniques = 0;
  96. max_len = 0;
  97. int begin_len = 0;
  98. int increment_mode = 0;
  99. int init_digit = 0;
  100. int ab = 0;
  101. int nn = 0;
  102. char* start;
  103. start = (char*)malloc(1);
  104. start[0] = '\0';
  105.  
  106. int c;
  107. char* p;
  108.  
  109. // parse args
  110. while ((c = getopt(argc, argv, "a:b:n:p:l:du:s:r:m:vf:h")) != -1) {
  111. switch (c) {
  112. case 'p':
  113. prime_check_s(optarg);
  114. exit(0);
  115. break;
  116. case 'v':
  117. verbose = true;
  118. break;
  119. case 'd':
  120. verbose = true;
  121. debug = true;
  122. break;
  123. case 's':
  124. start = (char*)realloc(start, strlen(optarg) + 1);
  125. strcpy(start, optarg);
  126. printf("Setting start to %s\n", start);
  127. max_len = strlen(start);
  128. break;
  129. case 'm':
  130. increment_mode = strtoul(optarg, &p, 10);
  131. if (increment_mode == 0 || increment_mode > 3) {
  132. printf("unrecognized increment mode, valid modes are 0,1,2,3\n");
  133. usage();
  134. exit(1);
  135. }
  136. break;
  137. case 'b':
  138. begin_len = strtoul(optarg, &p, 10);
  139. if (begin_len == 0) {
  140. printf("beginning length must be greater than 0\n");
  141. usage();
  142. exit(1);
  143. }
  144. if (begin_len % 2 == 0) {
  145. // must be odd
  146. begin_len++;
  147. }
  148. break;
  149. case 'l':
  150. max_len = strtoul(optarg, &p, 10);
  151. if (max_len < 3) max_len = 3;
  152. if (verbose) printf("Setting length to %d\n", max_len);
  153. break;
  154. case 'a':
  155. ab = strtoul(optarg, &p, 10);
  156. if ((ab < 10) || (ab > 99)) {
  157. printf("a must be two digits\n");
  158. exit(1);
  159. }
  160. if (verbose) printf("Setting a to %d\n", ab);
  161. break;
  162. case 'n':
  163. nn = strtoul(optarg, &p, 10);
  164. if (nn < 1) {
  165. printf("1 must > 0\n");
  166. exit(1);
  167. }
  168. if (verbose) printf("Setting n to %d\n", nn);
  169. break;
  170. case 'f':
  171. init_digit = (int)optarg[0] - 48;
  172. if ((init_digit < 1) || (init_digit > 9)) {
  173. printf("invalid start digit: %s", optarg);
  174. exit(1);
  175. }
  176. break;
  177. case 'u':
  178. uniques = labs(strtol(optarg, &p, 10));
  179. if (uniques == 0) {
  180. printf("uniques must be greater than 0\n");
  181. usage();
  182. exit(1);
  183. }
  184. if (verbose) printf("Setting uniques to %d\n", uniques);
  185. break;
  186. case 'r':
  187. iterations = labs(strtol(optarg, &p, 10));
  188. if (iterations < 5) {
  189. printf(
  190. "Setting number of repetitions to %d increases chances of "
  191. "finding pseudoprimes as length grows.\n",
  192. iterations);
  193. }
  194. if (verbose)
  195. printf(
  196. "Setting number of repetitions of probabilistic primality check "
  197. "to %d\n",
  198. iterations);
  199. break;
  200. case 'h':
  201. usage();
  202. exit(0);
  203. }
  204. }
  205.  
  206. if ((ab > 0) && (nn > 0)) {
  207. // check (ab*10^len-ba)/99
  208. int tmp = ab;
  209. int a = tmp % 10;
  210. tmp = tmp / 10;
  211. int b = tmp % 10;
  212. printf("checking (%d*10^%d-%d%d)/99\n", ab, nn, a, b);
  213. bool prime = check_supp(ab, nn);
  214. if (prime) {
  215. printf("probable prime\n");
  216. } else {
  217. printf("not prime\n");
  218. }
  219. exit(0);
  220. }
  221. if (max_len == 0) {
  222. usage();
  223. exit(1);
  224. } else if (max_len % 2 == 0) {
  225. max_len++;
  226. if (verbose) printf("even start length - bumping to %d\n", max_len);
  227. }
  228. if (begin_len > max_len) {
  229. max_len = begin_len;
  230. printf(
  231. "beginning len set greater than end len, setting beginning len to %d\n",
  232. begin_len);
  233. }
  234. // set function pointer depending on mode and initial value
  235. void (*cand_increment)(mpz_t n);
  236. if (increment_mode == 0) {
  237. // simple increment
  238. if (begin_len > 0) {
  239. printf("begin length ignored\n");
  240. }
  241. int start_len = strlen(start);
  242. cand_increment = primepal_cand_inc;
  243. if (start_len > 0) {
  244. // we got a starting number, check and set
  245. int err = mpz_set_str(n, start, 10);
  246. if (err) {
  247. printf("format error with start: %s\n", start);
  248. exit(1);
  249. }
  250. if (start_len % 2 == 0) {
  251. printf("start digit's even - use integer with odd number of digits.\n");
  252. exit(1);
  253. }
  254. } else {
  255. // we got a length to start with, set start
  256. if (init_digit > 0) {
  257. mpz_ui_pow_ui(n, 10, max_len);
  258. mpz_sub_ui(n, n, 1);
  259. mpz_fdiv_q_ui(n, n, 9);
  260. mpz_mul_ui(n, n, init_digit);
  261. } else {
  262. mpz_ui_pow_ui(n, 10, max_len - 1);
  263. mpz_add_ui(n, n, 1);
  264. }
  265. }
  266. } else if ((increment_mode == 1) || (increment_mode == 3)) {
  267. // SUPP - abababa
  268. cand_increment = undulating_cand_inc;
  269. // we got a length to start with, set start
  270. if (init_digit > 0) {
  271. printf("initial digit ignored\n");
  272. }
  273. if (begin_len > 0) {
  274. mpz_ui_pow_ui(n, 10, begin_len);
  275. mpz_sub_ui(n, n, 1);
  276. mpz_fdiv_q_ui(n, n, 9);
  277. } else {
  278. int err = mpz_set_str(n, "101", 10);
  279. assert(err == 0);
  280. }
  281. if (increment_mode == 1) {
  282. if (verbose) printf("in SUPP mode\n");
  283. }
  284. if (increment_mode == 3) {
  285. lazy = true;
  286. if (verbose) printf("in lazy SUPP mode - skipping some lengths\n");
  287. }
  288.  
  289. } else if (increment_mode == 2) {
  290. // longer period - abbbabbba
  291. cand_increment = long_period_cand_inc;
  292. if (begin_len > 0) {
  293. mpz_ui_pow_ui(n, 10, begin_len);
  294. mpz_sub_ui(n, n, 1);
  295. mpz_fdiv_q_ui(n, n, 9);
  296. } else {
  297. int err = mpz_set_str(n, "16661", 10);
  298. assert(err == 0);
  299. }
  300. if (verbose) printf("in alternate period mode\n");
  301. } else {
  302. printf("error in getting mode - use 1,2,3\n");
  303. usage();
  304. exit(1);
  305. }
  306. free(start);
  307.  
  308. if (verbose) {
  309. gmp_printf("starting with: %Zd\n\n", n);
  310. }
  311. start_time = get_time();
  312. factor_time = 0;
  313. // kick off function to do work, return when done.
  314. (*cand_increment)(n);
  315.  
  316. if (verbose) {
  317. printf("prime candidates tested: %lld\n", candidates);
  318. printf("primes found: %lld\n", primes);
  319. printf("ratio: %f %%\n", (((double)primes) / ((double)candidates)) * 100.0);
  320. }
  321. uint64_t end_time = get_time();
  322. float duration = end_time - start_time;
  323. float duration_sec = (float)duration / 1000;
  324. float factor_sec = (float)factor_time / 1000;
  325. printf("factoring time: %fs\n", factor_sec);
  326. printf("elapsed time: %fs\n", duration_sec);
  327.  
  328. /*
  329. if (debug) {
  330. // we keep track of some composite factors
  331. // useful for finding candidates to try to skip
  332. for (int i = 1; i < 999; i += 2) {
  333. if (mod_cnt[i] > 0) {
  334. printf("%d count: %llu\n", i, mod_cnt[i]);
  335. }
  336. }
  337. }*/
  338. }
  339.  
  340. void primepal_cand_inc(mpz_t n) {
  341. // for smaller numbers we spend a lot more time building these than factoring
  342. // precalculate to shave time
  343.  
  344. // skip when first digit is 2,4,5,6,8
  345. static const uint_fast8_t hop[10] = {0, 0, 1, 0, 1, 1, 1, 0, 1, 0};
  346.  
  347. mpz_t tmp, ne10, ne11;
  348. mpz_init(tmp);
  349. mpz_init(ne10);
  350. mpz_init(ne11);
  351.  
  352. // calc the things we only need to calc once here before we loop
  353. // sometimes mpz_sizeinbase reports len+1
  354. uint_fast32_t len = mpz_sizeinbase(n, 10);
  355. if (len % 2 == 0) {
  356. len--;
  357. }
  358. const uint_fast32_t mid = len / 2;
  359. const uint_fast8_t uniqs = uniques;
  360.  
  361. // ne10 = 10^len-1
  362. mpz_ui_pow_ui(ne10, 10, len - 1);
  363. mpz_ui_pow_ui(ne11, 10, len - 1);
  364. mpz_add_ui(ne11, ne11, 1);
  365.  
  366. //precalc n^10,100,1000, etc.
  367. mpz_t n_powers[mid];
  368. for (int i = 0; i<=mid; i++){
  369. mpz_init(n_powers[i]);
  370. mpz_ui_pow_ui(n_powers[i], 10, i);
  371. }
  372. while (1) {
  373. /*
  374. increment based on context:
  375.  
  376. if checking uniques, grab first half of digits and increment until we get
  377. to the next candidate that has the correct number of unique digits, covert
  378. to string, reverse string, add reversed back in.
  379.  
  380. if straight iteration, check if mid digit < 9, if so add 10^len/2 to
  381. increment, else convert to string and reverse.
  382.  
  383. if first digit is 2,4,5,6,8 skip range, since can't be prime.
  384.  
  385. once we're past max_len we return
  386. */
  387. if (uniqs) {
  388. // cut len in half, increment as much as needed, then make string
  389. mpz_tdiv_q(n, n, n_powers[mid]);
  390. mpz_add_ui(n, n, 1);
  391. int_fast8_t i;
  392. while ((i = interesting(n)) != -1) {
  393. mpz_add(n, n, n_powers[i]);
  394. }
  395. } else {
  396. // 90% of time skip stringy slowness
  397. // get first digit (basically n0,000 / 10,000)
  398. mpz_fdiv_q(tmp, n, ne10);
  399. const uint_fast8_t first = mpz_get_ui(tmp);
  400. if (hop[first]) {
  401. // if first digit is 2,4,5,6,8 can't be prime, hop to n0...0n
  402. mpz_mul_ui(n, ne11, first + 1);
  403. }
  404. // get last digit (n/(10^len/2))%10
  405. mpz_tdiv_q(tmp, n, n_powers[mid]);
  406. const uint_fast8_t last = mpz_mod_ui(tmp, tmp, 10);
  407. if (last != 9) {
  408. // we can just increment the mid-digit n+=(10^len/2) and we're done
  409. mpz_add(n, n, n_powers[mid]);
  410. if (mpz_gcd_ui(tmp, n, 22141119) <= 1) {
  411. prime_check(n);
  412. }
  413. continue;
  414. } else {
  415. // have to make string - cut len in half and inc by one
  416. mpz_tdiv_q(n, n, n_powers[mid]);
  417. mpz_add_ui(n, n, 1);
  418. }
  419. }
  420.  
  421. // grab first half in str, reverse
  422. char last_half[mid + 1];
  423. mpz_get_str(last_half, 10, n);
  424. // check len of incremented cand.
  425. // if we've bumped len, last_half[mid + 1] will be a num.
  426. if (last_half[mid + 1] != '\0') {
  427. return;
  428. }
  429. // check first digit
  430. const uint_fast8_t first = last_half[0] - 48;
  431. if (hop[first]) {
  432. // if first digit is 2,4,5,6,8 can't be prime, hop to n0...0n
  433. // n = (10^len)*(first+1)+(first+1)
  434. mpz_mul_ui(n, ne11, first + 1);
  435. } else {
  436. // chop last digit and reverse
  437. last_half[mid] = '\0';
  438. char *c1, *c2;
  439. for (c1 = last_half, c2 = last_half + mid - 1; c2 > c1; ++c1, --c2) {
  440. *c1 ^= *c2;
  441. *c2 ^= *c1;
  442. *c1 ^= *c2;
  443. }
  444. // last_half from str to bignum
  445. assert(mpz_set_str(tmp, last_half, 10) == 0);
  446. // multiply n by 10^mid to get n to full length with latter half zeroed
  447. // then add first half to reversed half to get next palindrome
  448. mpz_mul(n, n, n_powers[mid]);
  449. mpz_add(n, n, tmp);
  450. }
  451.  
  452. // factors that have a higher freq, check if divisible before proceeding
  453. // to full check
  454. // best so far:
  455. // 3*7*11*13*73*101 = 22141119
  456. if (mpz_gcd_ui(tmp, n, 22141119) <= 1) {
  457. prime_check(n);
  458. }
  459. }
  460. }
  461.  
  462. int interesting(mpz_t n) {
  463. // filter for x unique digits in number
  464. // returns int
  465. // -1=interesting
  466. // else return power of ten for next increment
  467.  
  468. // convert to str to scan digits
  469. uint_fast32_t len = (max_len/2)+1;
  470. const uint_fast8_t uniqs = uniques;
  471. char str[len];
  472. mpz_get_str(str, 10, n);
  473. //uint_fast8_t len = strlen(str);
  474. uint_fast8_t found = 0;
  475. bool count[10] = {false};
  476. // walk digits, tracking freq.
  477. // once we hit > uniqs return place to allow faster increment
  478. for (int i = 0; i <= len; i++) {
  479. if (count[str[i] - 48] == false) {
  480. found++;
  481. if (found > uniqs) {
  482. // we've hit the place to increment
  483. return len - i - 1;
  484. }
  485. count[str[i] - 48] = true;
  486. }
  487. }
  488. return -1;
  489. }
  490.  
  491. void undulating_cand_inc(mpz_t n) {
  492. // build SUPP xyxyxyx candidates (also 11...11 repunits)
  493. mpz_t tmp;
  494. mpz_init(tmp);
  495. int undula_arr[100];
  496. undularray(undula_arr);
  497. bool done = false;
  498. while (done == false) {
  499. int len = mpz_sizeinbase(n, 10);
  500. // get len, if even decrement
  501. if (len % 2 == 0) {
  502. len--;
  503. }
  504. if (len > max_len) {
  505. done = true;
  506. mpz_clear(tmp);
  507. return;
  508. }
  509.  
  510. // get first two digits
  511. mpz_ui_pow_ui(tmp, 10, len - 2);
  512. mpz_tdiv_q(tmp, n, tmp);
  513. int first_two = mpz_get_ui(tmp);
  514.  
  515. // used to check if candidate is %3
  516. int firsts = len / 2 + 1;
  517. int seconds = firsts - 1;
  518. int first, second;
  519.  
  520. // increment - walk undula_arr to skip composites
  521. bool found = false;
  522. while (found == false) {
  523. if (first_two == 98) {
  524. // this gets more complex due to laziness
  525. if (lazy == true) {
  526. int next_len = lazy_hop(len);
  527. if (next_len > 0) {
  528. len = next_len;
  529. } else if (next_len == -1) {
  530. // past the point we have values, shift to slow mode
  531. printf("\n\n");
  532. lazy = false;
  533. len = len + 2;
  534. firsts++;
  535. seconds++;
  536. } else {
  537. len = len + 2;
  538. firsts++;
  539. seconds++;
  540. }
  541. } else {
  542. len = len + 2;
  543. firsts++;
  544. seconds++;
  545. }
  546. first_two = 10;
  547. }
  548. first_two = undula_arr[first_two];
  549.  
  550. // hop some composite numbers that are periodic
  551. if (first_two == 11) {
  552. // these are %41==0
  553. if (len % 5 == 0) first_two = 12;
  554. }
  555. if (first_two == 12) {
  556. if (((len + 25) % 30) == 0) first_two = 13;
  557. } else if (first_two == 14) {
  558. if (((len + 7) % 30) == 0) first_two = 16;
  559. } else if (first_two == 16) {
  560. if ((((len) / 3) % 2) == 1) first_two = 17;
  561. } else if (first_two == 18) {
  562. //% 31
  563. if (((len + 23) % 30) == 0) first_two = 19;
  564. } else if (first_two == 34) {
  565. // never prime if len-3%3==0
  566. //(34*10n-43)/99 - only primes
  567. if (len % 3 == 0) first_two = 35;
  568. if (len % 3 == 1) first_two = 35;
  569. } else if (first_two == 35) {
  570. // these are %211==0
  571. if ((len - 9) % 10 == 0) first_two = 37;
  572. } else if (first_two == 37) {
  573. if ((((len - 2) / 3) % 2) == 1) first_two = 38;
  574. } else if (first_two == 71) {
  575. // biggest hop here
  576. if ((len - 1) % 3 == 0) {
  577. first_two = 91;
  578. }
  579. } else if (first_two == 75) {
  580. //% 31
  581. if (((len + 9) % 30) == 0) first_two = 78;
  582. } else if (first_two == 78) {
  583. //% 31
  584. if (((len + 19) % 30) == 0) first_two = 79;
  585. } else if (first_two == 92) {
  586. //% 31
  587. if (((len + 9) % 30) == 0) first_two = 94;
  588. } else if (first_two == 94) {
  589. // 94%3==0 is %13
  590. if (len % 3 == 0) first_two = 95;
  591. } else if (first_two == 95) {
  592. //% 31
  593. if (((len + 7) % 30) == 0) first_two = 97;
  594. }
  595.  
  596. // check mod3, if clear then check mod11 based on digits before
  597. // building gmp bignum
  598. first = first_two;
  599. second = first % 10;
  600. first /= 10;
  601. uint64_t digit_sum = firsts * first + seconds * second;
  602. if (digit_sum % 3 != 0) {
  603. uint64_t digit_diff = abs(firsts * first - seconds * second);
  604. if (digit_diff % 11 != 0) {
  605. found = true;
  606. }
  607. }
  608. }
  609.  
  610. int ab = first * 10 + second;
  611. int ba = second * 10 + first;
  612.  
  613. //(ab*10^len-ba)/99 = aba...aba
  614. mpz_ui_pow_ui(tmp, 10, len);
  615. mpz_mul_ui(n, tmp, ab);
  616. mpz_sub_ui(n, n, ba);
  617. mpz_tdiv_q_ui(n, n, 99);
  618.  
  619. // skip lower digits, mpz_gcd_ui deals with <49
  620. // 61*71*79*97*239=7932040267
  621. // 61*71*73*79*97*239=579038939491
  622. // factors that have a higher freq, check if divisible before proceeding to
  623. // check
  624. // probably could use some tuning
  625. if (mpz_gcd_ui(tmp, n, 579038939491) <= 1) {
  626. if (prime_check(n) && verbose) {
  627. printf("(%d*10^%d-%d)/99\n\n", ab, len, ba);
  628. }
  629. }
  630. }
  631. // done
  632. }
  633.  
  634. void long_period_cand_inc(mpz_t n) {
  635. mpz_t tmp;
  636. mpz_init(tmp);
  637. bool done = false;
  638. while (done == false) {
  639. int len = mpz_sizeinbase(n, 10);
  640. // get len, if even decrement
  641. if (len % 2 == 0) {
  642. len--;
  643. }
  644. if (len > max_len) {
  645. done = true;
  646. mpz_clear(tmp);
  647. return;
  648. }
  649.  
  650. // get first two digits
  651. mpz_ui_pow_ui(tmp, 10, len - 2);
  652. mpz_tdiv_q(tmp, n, tmp);
  653. int first_two = mpz_get_ui(tmp);
  654. int first = first_two / 10;
  655. int second = first_two % 10;
  656.  
  657. // given abbbabbba...abbbabbba, as,bs = instances of a and b
  658. int as = ((len - 1) / 4) + 1;
  659. int bs = len - as;
  660.  
  661. bool found = false;
  662. while (found == false) {
  663. if (first_two == 19) {
  664. // 19..70==prime wasteland
  665. first_two = 71;
  666. } else {
  667. first_two++;
  668. }
  669.  
  670. if (first_two == 79) {
  671. len = len + 4;
  672. as = as + 1;
  673. bs = bs + 3;
  674. first_two = 11;
  675. }
  676. first_two = long_period_composite_chk(first_two, as);
  677. // get new first,second
  678. second = first_two % 10;
  679. first = first_two / 10;
  680.  
  681. // check mod 3 without gmp, if clear then check mod 11 without gmp
  682. // if both clear, then build the bignum and check primality
  683. uint64_t digit_sum = as * first + bs * second;
  684. if (digit_sum % 3 != 0) {
  685. // check mod 11 - add up odd digits & even digits
  686. // subtract evens from odds and check for divisibility by 11
  687. int odds = (first * as) + (as - 1) * second;
  688. int evens = (len / 2) * second;
  689. if (abs(odds - evens) % 11 != 0) {
  690. found = true;
  691. }
  692. }
  693. }
  694. //(abbb*10^n-bbba)/9999 = abbbabbba...abbbabbba
  695. int abbb = first * 1000 + second * 100 + second * 10 + second;
  696. int bbba = second * 1000 + second * 100 + second * 10 + first;
  697. mpz_ui_pow_ui(tmp, 10, len);
  698. mpz_mul_ui(n, tmp, abbb);
  699. mpz_sub_ui(n, n, bbba);
  700. mpz_tdiv_q_ui(n, n, 9999);
  701.  
  702. // 3*7*11*61*101=1423191
  703. // 3*7*11*19*61*101=27040629
  704. // 3*7*11*61*101*449=639012759
  705. // 7*61*101*239*449=4628001497
  706.  
  707. // skip lower digits, mpz_gcd_ui deals with <49
  708. // 61*101*239*449=661143071
  709. // factors that have a higher freq, check if divisible before proceeding to
  710. // check
  711. // probably could use some tuning
  712. if (mpz_gcd_ui(tmp, n, 661143071) <= 1) {
  713. if (prime_check(n) && verbose) {
  714. printf("(%d*10^%d-%d)/9999\n\n", abbb, len, bbba);
  715. }
  716. }
  717. }
  718. // done
  719. }
  720.  
  721. void mod_chk(int fac, int len, int as, mpz_t n) {
  722. // log whether fac is a factor
  723. // useful for finding composites in undulating candidates to skip
  724. if (mpz_divisible_ui_p(n, fac) > 0) {
  725. gmp_printf("mod %d - %Zd - %d \nas - %d\n", fac, n, len, as);
  726. }
  727. }
  728.  
  729. int long_period_composite_chk(int first_two, int as) {
  730. /*
  731. Should have documented these more as they filled in...
  732. there are various patterns that show up with some sets of
  733. first_two at various lengths that are always composite.
  734. -- based on tests, not proofs
  735. -- there may be cases where this skips primes
  736. */
  737.  
  738. if (first_two == 10) {
  739. // don't have proof, but these have been composite up to 20,000 digits
  740. first_two = 11;
  741. }
  742. if (first_two == 11) {
  743. if ((as + 15) % 13 == 0) {
  744. //%19==0
  745. first_two = 12;
  746. } else if ((as - 17) % 13 == 0) {
  747. //%19==0
  748. first_two = 12;
  749. } else if ((as + 22) % 13 == 0) {
  750. //%53==0
  751. first_two = 12;
  752. } else if ((as - 2) % 15 == 0) {
  753. //%71==0
  754. first_two = 12;
  755. }
  756. }
  757. if (first_two == 12) {
  758. if ((as + 17) % 21 == 0) {
  759. first_two = 13;
  760. } else if ((as + 10) % 7 == 0) {
  761. first_two = 13;
  762. } else if ((as + 6) % 13 == 0) {
  763. first_two = 13;
  764. }
  765. }
  766. if (first_two == 13) {
  767. if ((as + 65) % 69 == 0) {
  768. first_two = 14;
  769. } else if ((as + 15) % 23 == 0) {
  770. first_two = 14;
  771. } else if ((as + 4) % 9 == 0) {
  772. first_two = 14;
  773. } else if ((as + 13) % 17 == 0) {
  774. first_two = 14;
  775. } else if ((as + 21) % 53 == 0) {
  776. first_two = 14;
  777. }
  778. }
  779. if (first_two == 14) {
  780. if ((as + 38) % 69 == 0) {
  781. first_two = 15;
  782. } else if ((as + 65) % 69 == 0) {
  783. first_two = 15;
  784. } else if ((as + 7) % 15 == 0) {
  785. first_two = 15;
  786. } else if ((as + 16) % 13 == 0) {
  787. first_two = 15;
  788. } else if ((as + 1) % 21 == 0) {
  789. first_two = 15;
  790. }
  791. }
  792. if (first_two == 15) {
  793. if ((as + 6) % 11 == 0) {
  794. //%23==0
  795. first_two = 16;
  796. } else if ((as + 34) % 45 == 0) {
  797. first_two = 16;
  798. } else if ((as + 4) % 45 == 0) {
  799. first_two = 16;
  800. } else if ((as + 4) % 15 == 0) {
  801. //%31
  802. first_two = 16;
  803. } else if ((as + 56) % 69 == 0) {
  804. first_two = 16;
  805. }
  806. }
  807. if (first_two == 16) {
  808. if (as % 3 == 0) {
  809. //%19==0
  810. first_two = 17;
  811. } else if ((as + 1) % 9 == 0) {
  812. first_two = 17;
  813. } else if ((as + 3) % 23 == 0) {
  814. first_two = 17;
  815. } else if ((as + 5) % 21 == 0) {
  816. first_two = 17;
  817. } else if ((as + 24) % 29 == 0) {
  818. first_two = 17;
  819. } else if ((as + 49) % 53 == 0) {
  820. // mod107
  821. first_two = 17;
  822. }
  823. }
  824. if (first_two == 17) {
  825. if ((as + 4) % 11 == 0) {
  826. first_two = 18;
  827. } else if ((as + 1) % 3 == 0) {
  828. //%13
  829. first_two = 18;
  830. } else if ((as + 11) % 21 == 0) {
  831. first_two = 18;
  832. } else if ((as + 2) % 39 == 0) {
  833. first_two = 18;
  834. } else if ((as + 53) % 69 == 0) {
  835. first_two = 18;
  836. } else if ((as + 14) % 29 == 0) {
  837. first_two = 18;
  838. } else if ((as + 4) % 53 == 0) {
  839. //%107
  840. first_two = 18;
  841. }
  842. }
  843. if (first_two == 18) {
  844. if ((as - 3) % 4 == 0) {
  845. first_two = 19;
  846. } else if ((as + 8) % 21 == 0) {
  847. first_two = 19;
  848. } else if ((as + 1) % 21 == 0) {
  849. first_two = 19;
  850. } else if ((as + 24) % 13 == 0) {
  851. first_two = 19;
  852. } else if ((as + 43) % 69 == 0) {
  853. first_two = 19;
  854. } else if ((as + 20) % 69 == 0) {
  855. first_two = 19;
  856. }
  857. }
  858. if (first_two == 19) {
  859. if ((as + 5) % 11 == 0) {
  860. first_two = 70;
  861. } else if ((as + 1) % 13 == 0) {
  862. first_two = 70;
  863. } else if ((as + 10) % 21 == 0) {
  864. first_two = 70;
  865. } else if ((as + 8) % 23 == 0) {
  866. first_two = 70;
  867. } else if ((as + 11) % 29 == 0) {
  868. first_two = 70;
  869. } else if ((as + 4) % 5 == 0) {
  870. first_two = 70;
  871. }
  872. }
  873. if (first_two == 70) {
  874. first_two = 71;
  875. }
  876. if (first_two == 71) {
  877. if ((as + 7) % 9 == 0) {
  878. first_two = 72;
  879. } else if ((as + 7) % 13 == 0) {
  880. first_two = 72;
  881. } else if ((as + 2) % 23 == 0) {
  882. first_two = 72;
  883. } else if ((as + 83) % 105 == 0) {
  884. first_two = 72;
  885. } else if ((as + 38) % 53 == 0) {
  886. //%107
  887. first_two = 72;
  888. }
  889. }
  890. if (first_two == 72) {
  891. if ((as + 4) % 39 == 0) {
  892. first_two = 73;
  893. } else if ((as + 7) % 69 == 0) {
  894. first_two = 73;
  895. } else if ((as + 4) % 9 == 0) {
  896. first_two = 73;
  897. } else if ((as + 8) % 15 == 0) {
  898. //%31
  899. first_two = 73;
  900. } else if ((as + 88) % 105 == 0) {
  901. first_two = 73;
  902. } else if ((as + 15) % 53 == 0) {
  903. //%107
  904. first_two = 73;
  905. }
  906. }
  907. if (first_two == 73) {
  908. if ((as + 4) % 11 == 0) {
  909. first_two = 74;
  910. } else if ((as + 18) % 23 == 0) {
  911. first_two = 74;
  912. } else if ((as + 4) % 17 == 0) {
  913. first_two = 74;
  914. } else if ((as + 27) % 29 == 0) {
  915. first_two = 74;
  916. }
  917. }
  918. if (first_two == 74) {
  919. if (as % 3 == 0) {
  920. //%19==0
  921. first_two = 75;
  922. } else if ((as + 2) % 11 == 0) {
  923. first_two = 75;
  924. } else if ((as + 4) % 7 == 0) {
  925. first_two = 75;
  926. } else if ((as + 16) % 7 == 0) {
  927. first_two = 75;
  928. } else if ((as + 4) % 15 == 0) {
  929. first_two = 75;
  930. } else if ((as + 1) % 9 == 0) {
  931. first_two = 75;
  932. } else if ((as + 1) % 23 == 0) {
  933. first_two = 75;
  934. } else if ((as + 11) % 29 == 0) {
  935. first_two = 75;
  936. } else if ((as + 17) % 21 == 0) {
  937. first_two = 75;
  938. }
  939. }
  940. if (first_two == 75) {
  941. if ((as + 1) % 4 == 0) {
  942. //%19==0
  943. first_two = 76;
  944. } else if ((as + 2) % 9 == 0) {
  945. //%19==0
  946. first_two = 76;
  947. } else if ((as + 25) % 29 == 0) {
  948. //%19==0
  949. first_two = 76;
  950. } else if ((as + 10) % 15 == 0) {
  951. // mod 71, needs review
  952. first_two = 76;
  953. } else if ((as + 5) % 15 == 0) {
  954. // mod 71, needs review
  955. first_two = 76;
  956. }
  957. }
  958. if (first_two == 76) {
  959. if ((as + 1) % 15 == 0) {
  960. //%19==0
  961. first_two = 78;
  962. } else if ((as + 1) % 21 == 0) {
  963. first_two = 78;
  964. } else if ((as + 11) % 33 == 0) {
  965. //%67
  966. first_two = 78;
  967. } else if ((as + 7) % 41 == 0) {
  968. //%83
  969. first_two = 78;
  970. } else if ((as + 37) % 53 == 0) {
  971. //%107
  972. first_two = 78;
  973. }
  974. }
  975. if (first_two == 77) {
  976. first_two = 78;
  977. }
  978. // for 78 we just have to let it go
  979. return first_two;
  980. }
  981.  
  982. bool prime_check(mpz_t n) {
  983. candidates++;
  984. uint64_t before_chk = get_time();
  985. uint_fast8_t p = mpz_probab_prime_p(n, iterations);
  986. factor_time += get_time() - before_chk;
  987. if (p > 0) {
  988. gmp_printf("%Zd\n", n);
  989. primes++;
  990. return true;
  991. }/* else if (debug) {
  992. divisor_count(n);
  993. }*/
  994. return false;
  995. }
  996.  
  997. uint64_t get_time() {
  998. struct timeval tv;
  999. uint64_t ret;
  1000. gettimeofday(&tv, NULL);
  1001. // get microseconds, convert to milliseconds
  1002. ret = tv.tv_usec;
  1003. ret /= 1000;
  1004. // get seconds, convert to milliseconds
  1005. ret += (tv.tv_sec * 1000);
  1006. return ret;
  1007. }
  1008.  
  1009. bool mini_chk(mpz_t n) {
  1010. // check smaller primes, return true if we aren't divisible by them
  1011. // return immediately if we find we're composite
  1012. // this is similar to what mpz_probab_prime_p() does but we go a
  1013. // little past %53
  1014.  
  1015. // poor perf. once numbers are >1000 digits
  1016. for (int i = 1; i < 27; i++) {
  1017. if (mpz_divisible_ui_p(n, divisors[i]) > 0) {
  1018. return false;
  1019. }
  1020. }
  1021. return true;
  1022. }
  1023.  
  1024. void divisor_count(mpz_t n) {
  1025. // count divisors to find any that are interesting
  1026. for (int i = 1; i < 164; i++) {
  1027. if (mpz_divisible_ui_p(n, divisors[i]) > 0) {
  1028. mod_cnt[divisors[i]]++;
  1029. }
  1030. }
  1031. }
  1032.  
  1033. void undularray(int undula_arr[100]) {
  1034. // it works, stop judging me
  1035.  
  1036. // build array to loop for next candidate
  1037. // skip matching digits if not 1
  1038. // skip leading 20s,40s,50s,60s,80s
  1039. // skip 76, always factorable by either 3, 7, or 13
  1040. // skip 77, always factorable
  1041. undula_arr[10] = 11;
  1042. undula_arr[11] = 12;
  1043. undula_arr[12] = 13;
  1044. undula_arr[13] = 14;
  1045. undula_arr[14] = 15;
  1046. undula_arr[15] = 16;
  1047. undula_arr[16] = 17;
  1048. undula_arr[17] = 18;
  1049. undula_arr[18] = 19;
  1050. undula_arr[19] = 31;
  1051. undula_arr[31] = 32;
  1052. undula_arr[32] = 34;
  1053. undula_arr[34] = 35;
  1054. undula_arr[35] = 37;
  1055. undula_arr[37] = 38;
  1056. undula_arr[38] = 71;
  1057. undula_arr[71] = 72;
  1058. undula_arr[72] = 73;
  1059. undula_arr[73] = 74;
  1060. undula_arr[74] = 75;
  1061. undula_arr[75] = 78;
  1062. undula_arr[78] = 79;
  1063. undula_arr[79] = 91;
  1064. undula_arr[91] = 92;
  1065. undula_arr[92] = 94;
  1066. undula_arr[94] = 95;
  1067. undula_arr[95] = 97;
  1068. undula_arr[97] = 98;
  1069. undula_arr[98] = 10;
  1070. }
  1071.  
  1072. int lazy_hop(int lazy) {
  1073. // skip lengths with no SUPPs
  1074. // yes, this is silly in many ways, but it helps to skip past the shorter
  1075. // range to grind through to the bigger ones
  1076. int ret = 0;
  1077. switch (lazy) {
  1078. case 627:
  1079. ret = 849;
  1080. break;
  1081. case 849:
  1082. ret = 857;
  1083. break;
  1084. case 857: // 857 is prime
  1085. ret = 961;
  1086. break;
  1087. case 961:
  1088. ret = 1031; // 1031 is prime
  1089. break;
  1090. case 1031:
  1091. ret = 1079;
  1092. break;
  1093. case 1079:
  1094. ret = 1295;
  1095. break;
  1096. case 1295:
  1097. ret = 1399;
  1098. break;
  1099. case 1399: // 1399 is prime
  1100. ret = 1415;
  1101. break;
  1102. case 1415:
  1103. ret = 1479;
  1104. break;
  1105. case 1479:
  1106. ret = 1597;
  1107. break;
  1108. case 1597: // 1597 is prime
  1109. ret = 1631;
  1110. break;
  1111. case 1631:
  1112. ret = 1791;
  1113. break;
  1114. case 1791:
  1115. ret = 1893;
  1116. break;
  1117. case 1893:
  1118. ret = 1969;
  1119. break;
  1120. case 1969:
  1121. ret = 1979;
  1122. break;
  1123. case 1979: // 1979 is prime
  1124. ret = 2069;
  1125. break;
  1126. case 2069: // 2069 is prime
  1127. ret = 2075;
  1128. break;
  1129. case 2075:
  1130. ret = 2165;
  1131. break;
  1132. case 2165:
  1133. ret = 2177;
  1134. break;
  1135. case 2177:
  1136. ret = 2335;
  1137. break;
  1138. case 2335:
  1139. ret = 2529;
  1140. break;
  1141. case 2529:
  1142. ret = 2883;
  1143. break;
  1144. case 2883:
  1145. ret = 3015;
  1146. break;
  1147. case 3015:
  1148. ret = 3047;
  1149. break;
  1150. case 3047:
  1151. ret = 3057;
  1152. break;
  1153. case 3057:
  1154. ret = 3147;
  1155. break;
  1156. case 3147:
  1157. ret = 3407;
  1158. break;
  1159. case 3407: // 3407 is prime
  1160. ret = 3503;
  1161. break;
  1162. case 3503:
  1163. ret = 3657;
  1164. break;
  1165. case 3657:
  1166. ret = 4157;
  1167. break;
  1168. case 4157: // 4157 is prime
  1169. ret = 4233;
  1170. break;
  1171. case 4233:
  1172. ret = 4335;
  1173. break;
  1174. case 4335:
  1175. ret = 4573;
  1176. break;
  1177. case 4573:
  1178. ret = 4859;
  1179. break;
  1180. case 4859:
  1181. ret = 4885;
  1182. break;
  1183. case 4885:
  1184. ret = 6249;
  1185. break;
  1186. case 6249:
  1187. ret = 6343;
  1188. break;
  1189. case 6343:
  1190. ret = 6815;
  1191. break;
  1192. case 6815:
  1193. ret = 6959;
  1194. break;
  1195. case 6959:
  1196. ret = 7233;
  1197. break;
  1198. case 7233:
  1199. ret = 7387;
  1200. break;
  1201. case 7387:
  1202. ret = 7809;
  1203. break;
  1204. case 7809:
  1205. ret = 8315;
  1206. break;
  1207. case 8315:
  1208. ret = 9273;
  1209. break;
  1210. case 9273:
  1211. ret = 9591;
  1212. break;
  1213. case 9591:
  1214. ret = 9599;
  1215. break;
  1216. case 9599:
  1217. ret = 10419;
  1218. break;
  1219. case 10419:
  1220. ret = 11393;
  1221. break;
  1222. case 11393: // 11393 is prime
  1223. ret = 11399;
  1224. break;
  1225. case 11399:
  1226. ret = 12849;
  1227. break;
  1228. case 12849:
  1229. ret = 13221;
  1230. break;
  1231. case 13221:
  1232. ret = 13265;
  1233. break;
  1234. case 13265:
  1235. ret = -1;
  1236. break;
  1237. }
  1238. return ret;
  1239. }
  1240.  
  1241. bool check_supp(int ab, int nn) {
  1242. // just for -a, -b flags
  1243. mpz_t n;
  1244. mpz_init(n);
  1245. mpz_t tmp;
  1246. mpz_init(tmp);
  1247. //(ab*10^len-ba)/99 = aba...aba
  1248. mpz_ui_pow_ui(tmp, 10, nn);
  1249. mpz_mul_ui(n, tmp, ab);
  1250. mpz_sub_ui(n, n, ab);
  1251. mpz_cdiv_q_ui(n, n, 99);
  1252. if (verbose) {
  1253. gmp_printf("%Zd\n", n);
  1254. }
  1255. int pp = mpz_probab_prime_p(n, 20);
  1256. if (pp != 0) {
  1257. return true;
  1258. }
  1259. return false;
  1260. }
  1261.  
  1262. void prime_check_s(char* s) {
  1263. // just for -p flag
  1264. mpz_t n;
  1265. mpz_init(n);
  1266. int err = mpz_set_str(n, s, 10);
  1267. if (err) {
  1268. printf("format error integer arg: %s\n", s);
  1269. exit(1);
  1270. }
  1271. int pp = mpz_probab_prime_p(n, 25);
  1272. if (pp == 0) {
  1273. puts("not prime.\n");
  1274. } else if (pp == 1) {
  1275. puts("probably prime.\n");
  1276. } else {
  1277. puts("prime.\n");
  1278. }
  1279. }
  1280.  
  1281. void usage() {
  1282. printf(
  1283. "usage: primepal (-p [number to test]) or (-s [starting number for "
  1284. "search])\n"
  1285. " -p [number] - check if individual number is prime\n"
  1286. " -s [start range]\n"
  1287. " [-m] 0 - (default) increment to next palindrome, if n, increment to "
  1288. "next palindrome with that number of digits\n"
  1289. " 1 - undulating mode (101, 18181, etc.), starting with 101, "
  1290. "ending when max_len set with -l is hit\n"
  1291. " 2 - large periodic increment mode, currently only three digits "
  1292. "(122212221, etc.)\n"
  1293. " 3 - impatient undulating mode - same as 1, but skips grinding "
  1294. "on some ranges in longer lengths in the 1000-10000 digit range where "
  1295. "there are no primes to be found\n"
  1296. " 4 - impatient periodic mode - same as 2, but skips grinding on "
  1297. "some ranges in longer lengths in the 1000-10000 digit range where "
  1298. "there are no primes to be found\n"
  1299. " [-d] limit search to primes with n unique digits (mode 0) \n"
  1300. " [-f] set all starting digits to num (mode 0)\n"
  1301. " [-r] set repetitions of primality test (default 15)\n"
  1302. " [-v] verbose\n"
  1303. " [-h] show this usage message\n");
  1304. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement