/* Kevin Ryde, January 2024 Usage: ./a.out This is a simple C program using the primesieve library and GNU C on a 64 bit system to generate terms of A146214 = lower prime of the 10^n'th twin prime pair A143585 = prime index of that lower of the 10^n'th pair A146536 = sum of first 10^n twin prime pairs (both of the primes in each pair) Output is free-form like (the 3 sequences respectively), n=3 twin prime 79561 index 7794 sum 69004076 The method is to get all primes with the primesieve library iterator and watch for twin primes, their indices, and sums, up to each target 10^n many twin prime pairs. Requires GNU C and 64 bits due to "__int128" used for the sum (at n=10, A146536(10) is about 76 bits). Terms to n=10 take about 70 minutes on a modest CPU at the time of writing. */ /* Set WANT_ASSERT to 1 for various assert() self-checks. Set WANT_DEBUG to 1 for verbose diagnostic prints. */ #define WANT_ASSERT 0 #define WANT_DEBUG 0 #if ! WANT_ASSERT #define NDEBUG #endif #include #include #include #include #include #include #include #include #include #include #include /* ---------------------------------------------------------------------- */ /* Generic Helpers */ /* GNU C */ typedef unsigned __int128 uint128_t; #if WANT_DEBUG #define DEBUG(expr) do { expr; } while (0) #else #define DEBUG(expr) do { } while (0) #endif #define numberof(array) (sizeof(array)/sizeof((array)[0])) #ifdef __GNUC__ #define LIKELY(cond) __builtin_expect((cond) != 0, 1) #define UNLIKELY(cond) __builtin_expect((cond) != 0, 0) #define ATTRIBUTE_PRINTF __attribute__ ((format (printf, 1, 2))) #else #define LIKELY(cond) (cond) #define UNLIKELY(cond) (cond) #define ATTRIBUTE_PRINTF #endif static noreturn ATTRIBUTE_PRINTF void error (const char *format, ...) { va_list ap; va_start (ap, format); vfprintf (stderr, format, ap); va_end (ap); exit(1); } void clock_gettime_or_die (clockid_t id, struct timespec *t, const char *name) { if (clock_gettime(id, t) != 0) error("Cannot clock_gettime() %s: %s\n", name, strerror(errno)); } double clock_gettime_double_sec_or_die (clockid_t id, const char *name) { struct timespec t; clock_gettime_or_die (id, &t, name); return ((double) t.tv_sec + t.tv_nsec / 1e9); } /* CPU time in seconds used since program start */ double cputime (void) { return (clock_gettime_double_sec_or_die (CLOCK_PROCESS_CPUTIME_ID, "CPUTIME")); } /* Return a string which is the decimal digits of n. Explicit conversion since __int128 is not in printf() etc. */ char * uint128_to_decimal_str (uint128_t n) { /* 3 decimal digits per byte of n, plus '\0' terminator */ static char buf[3*sizeof(n) + 1]; int i = sizeof(buf)-1; do { if (--i < 0) error("oops, uint128_to_decimal_str() buffer too small\n"); int digit = n % 10; n /= 10; buf[i] = digit + '0'; } while (n); return (buf + i); } /* str is a '\0'-terminated character string of decimal digits. Return their value, or error exit if bad string. */ uint128_t uint128_from_decimal_str (const char *str) { uint128_t n = 0; int c; while ((c = *str++)) { if (! isdigit(c)) error("value_from_decimal_str() invalid digit %d \"%c\"\n", c,c); c -= '0'; assert (0 <= c && c <= 9); n = 10*n + c; } return (n); } /* ---------------------------------------------------------------------- */ /* sample data for consistency checks */ const int primes_samples[] = { -1, /* no prime(0) */ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 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 }; int check_prime (uint64_t p, uint64_t num_primes) { if (num_primes < numberof(primes_samples)) { uint64_t want = primes_samples[num_primes]; if (p != want) error("oops, num_primes=%"PRIu64" want p=%"PRIu64" got p=%"PRIu64"\n", num_primes, want, p); } return (1); } static const uint64_t A146214_10pow_lower_twin_prime_samples[] = { 3, /* n=0 */ 107, 3821, 79559, 1260989, 18409199, 252427601, /* n=6 */ 3285916169, 41375648687, 507575862527, 6100479510551 /* n=10 */ }; /* prime index of the 10^n lower twin prime */ static const uint64_t A143585_10pow_lower_twin_prime_index_samples[] = { 2, /* n=0 */ 28, 530, 7793, 97255, 1175775, 13804822, /* n=6 */ 157523559, 1768264167, 19589296501, /* n=9 */ 214795459468 /* n=10 */ }; const char * const A146536_10pow_twin_primes_sum_samples[] = { "8", /* n=0 */ "908", "328184", "69004076", "11556327260", "1707243198956", "237064232862404", /* n=6 */ "31153163750203064", "3947120494191630260", /* n=8 */ "486665774050923191336", /* n=9 */ "58727077924563028184984" /* n=10 (corrected) */ }; /* p is the upper of a twin prime pair, p-2 is the lower. num_primes is how many to p inclusive, so num_primes-1 is index of the lower prime of the pair. */ void check_samples (int n, uint64_t p, uint64_t num_primes, uint128_t sum_twin_primes) { int bad = 0; if (n < numberof(A146214_10pow_lower_twin_prime_samples)) { uint64_t lower = p - 2; /* lower of the twin prime pair */ uint64_t want = A146214_10pow_lower_twin_prime_samples[n]; if (lower != want) { printf("oops, n=%d lower twin prime want %"PRIu64" got %"PRIu64"\n", n, want, lower); bad = 1; } } if (n < numberof(A143585_10pow_lower_twin_prime_index_samples)) { uint64_t got = num_primes - 1; /* index of the lower */ uint64_t want = A143585_10pow_lower_twin_prime_index_samples[n]; if (got != want) { printf("oops, n=%d lower twin prime index want %"PRIu64" got %"PRIu64"\n", n, want, got); bad = 1; } } if (n < numberof(A146536_10pow_twin_primes_sum_samples)) { const char *want_str = A146536_10pow_twin_primes_sum_samples[n]; uint128_t want = uint128_from_decimal_str(want_str); if (sum_twin_primes != want) { printf("oops, n=%d twin primes sum want %s got %s\n", n, want_str, uint128_to_decimal_str(sum_twin_primes)); bad = 1; } } if (bad) exit(1); } /* -------------------------------------------------------------------------- */ static int option_verbose = 0; void output (int target_n, uint64_t p, uint64_t num_primes, uint128_t sum_twin_primes) { printf("n=%d twin prime %"PRIu64" index %"PRIu64" sum %s\n", target_n, p, num_primes, uint128_to_decimal_str(sum_twin_primes)); if (option_verbose) printf(" at CPU time %.2lf minutes\n", cputime()/60.); check_samples (target_n, p, num_primes, sum_twin_primes); } void search (void) { uint64_t num_primes = 0; uint64_t num_twin_prime_pairs = 0; uint128_t sum_twin_primes = 0; /* both primes of each pair */ int target_n = 0; uint64_t target_10_pow_n = 1; /* 10^target_n */ uint64_t prev_prime_plus_2 = 0; primesieve_iterator it; primesieve_init (&it); for (;;) { uint64_t p = primesieve_next_prime (&it); num_primes++; assert (check_prime(p, num_primes)); if (UNLIKELY(p == prev_prime_plus_2)) { sum_twin_primes += 2*p - 2; /* both primes of the pair */ num_twin_prime_pairs++; if (num_twin_prime_pairs == target_10_pow_n) { output (target_n, p, num_primes, sum_twin_primes); target_n++; target_10_pow_n *= 10; } } prev_prime_plus_2 = p + 2; } } /* -------------------------------------------------------------------------- */ int main (int argc, char **argv) { setbuf(stdout,NULL); int i; for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"-v") == 0) { option_verbose++; } else { error("unrecognised command line option: %s\n", arg); } } if (option_verbose) { printf("primesieve version %s\n", primesieve_version()); } search (); return (0); }