/*
* Generate batches of consecutive prime numbers using a modified sieve of Eratosthenes
* algorithm that doesn't use much memory by using a windowing sieve buffer.
*
* Copyright (C) 2007 George Gesslein II.
*/
/*
Usage: matho-primes [start [stop]] ["twin"] ["pal" [base]]
Generate batches of consecutive prime numbers up to 19 digits.
If "twin" is specified, output only twin primes.
If "pal" is specified, output only palindromic primes.
The palindromic base may be specified, the default is base 10.
or
Usage: matho-primes [options] [start [stop]]
Generate batches of consecutive prime numbers up to 19 digits.
Options:
-t Output only twin primes.
-p base Output only palindromic primes.
*/
/*
* Changes:
*
* 11/22/05 - converted everything to long doubles.
* 3/25/06 - made primes buffer variable size.
*/
#include <stdio.h>
#include <stdlib.h>
#include <libgen.h>
#include <unistd.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include <errno.h>
#include <math.h>
#include <assert.h>
#define true 1
#define false 0
#ifndef max
#define max(a, b) (((a) > (b)) ? (a) : (b))
#endif
#ifndef min
#define min(a, b) (((a) < (b)) ? (a) : (b))
#endif
#define ARR_CNT(a) (sizeof(a)/sizeof(a[0])) /* returns the number of elements in an array */
#define MAX_K_INTEGER 1.0e19L /* 19 digits for long doubles */
/* Maximum memory usage in bytes; can be set to any size; the larger, the faster: */
#define BUFFER_SIZE 2000000
#if BUFFER_SIZE >= (INT_MAX / 2) || BUFFER_SIZE <= 0
#error "BUFFER_SIZE too large."
#endif
void generate_primes(void);
int test_pal(long double d, long double base);
void usage(void);
void usage2(void);
int get_long_double_int(char *cp, long double *dp);
long double strtold(const char *nptr, char **endptr);
long double start_value, number, end_value;
long double skip_multiples[] = { /* Additive array that skips over multiples of 2, 3, 5, and 7. */
10, 2, 4, 2, 4, 6, 2, 6,
4, 2, 4, 6, 6, 2, 6, 4,
2, 6, 4, 6, 8, 4, 2, 4,
2, 4, 8, 6, 4, 6, 2, 4,
6, 2, 6, 6, 4, 2, 4, 6,
2, 6, 4, 2, 4, 2,10, 2
}; /* sum of all numbers = 210 = (2*3*5*7) */
long double last_prime = -3.0;
int pal_flag, twin_flag;
long double pal_base = 10;
int buffer_size = BUFFER_SIZE; /* for variable size prime buffer */
char *prime; /* the boolean sieve array for the current window of numbers */
char *prog_name = "matho-primes";
int
main(int argc, char *argv[])
{
extern char *optarg;
extern int optind;
int i;
char buf[1000];
prog_name = strdup(basename(argv[0]));
start_value = -1.0;
end_value = MAX_K_INTEGER;
if (end_value == end_value + 1.0) {
fprintf(stderr, "MAX_K_INTEGER is too large!\n");
}
number = 0;
/* process command line options */
while ((i = getopt(argc, argv, "tp:")) >= 0) {
switch (i) {
case 't':
twin_flag = true;
break;
case 'p':
pal_flag = true;
if (!get_long_double_int(optarg, &pal_base)) {
usage2();
}
break;
default:
usage2();
}
}
if (argc > optind) {
if (strncmp(argv[optind], "twin", 4) == 0) {
twin_flag = true;
optind++;
}
}
if (argc > optind && isdigit(argv[optind][0])) {
if (get_long_double_int(argv[optind], &start_value)) {
optind++;
} else {
usage();
}
if (argc > optind && isdigit(argv[optind][0])) {
if (get_long_double_int(argv[optind], &end_value)) {
if (end_value < start_value) {
fprintf(stderr, "End value is less than start value.\n");
usage();
}
optind++;
number = MAX_K_INTEGER;
} else {
usage();
}
}
}
if (argc > optind) {
if (strncmp(argv[optind], "twin", 4) == 0) {
twin_flag = true;
optind++;
}
}
if (argc > optind) {
if (strncmp(argv[optind], "pal", 3) == 0) {
pal_flag = true;
optind++;
} else {
usage();
}
if (argc > optind) {
if (!get_long_double_int(argv[optind], &pal_base)) {
usage();
}
optind++;
}
}
if (argc > optind) {
usage();
}
if (pal_base < 2 || pal_base >= INT_MAX) {
fprintf(stderr, "Palindromic base must be >= 2.\n");
usage();
}
if (start_value < 0.0) {
get1:
fprintf(stderr, "Enter number (up to 19 digits) to start finding primes at (0): ");
if (fgets(buf, sizeof(buf), stdin) == NULL)
exit(1);
switch (buf[0]) {
case '\0':
case '\n':
case '\r':
start_value = 0;
break;
default:
if (!get_long_double_int(buf, &start_value)) {
goto get1;
}
}
}
if (number == 0) {
get2:
fprintf(stderr, "Enter number of primes to output (40): ");
if (fgets(buf, sizeof(buf), stdin) == NULL)
exit(1);
switch (buf[0]) {
case '\0':
case '\n':
case '\r':
number = 40;
break;
default:
if (!get_long_double_int(buf, &number)) {
goto get2;
}
}
}
buffer_size = (int) min(BUFFER_SIZE, (end_value - start_value) + 1.0L);
prime = (char *) malloc(buffer_size);
if (prime == NULL) {
fprintf(stderr, "%s: Not enough memory.\n", prog_name);
exit(2);
}
generate_primes();
exit(0);
}
/*
* Eliminate all multiples of "arg" from the sieve array.
*/
inline void
elim_factor(long double arg)
{
long double d;
int i, j;
d = ceill(start_value / arg);
if (d < 2.0)
d = 2.0;
d *= arg;
d -= start_value;
if (d >= buffer_size)
return;
i = (int) d;
assert(i >= 0);
if (arg >= buffer_size) {
prime[i] = false;
} else {
j = (int) arg;
for (; i < buffer_size; i += j) {
prime[i] = false;
}
}
}
/*
* Generate and display consecutive prime numbers.
*/
void
generate_primes()
{
int n, j;
long double count, ii, vv;
for (count = 0; count < number; start_value += buffer_size) {
if (start_value > end_value) {
return;
}
memset(prime, true, buffer_size); /* set the prime array to all true */
vv = 1.0 + sqrtl(min(start_value + (long double) buffer_size, end_value));
elim_factor(2.0L);
elim_factor(3.0L);
elim_factor(5.0L);
elim_factor(7.0L);
ii = 1.0;
for (;;) {
for (j = 0; j < ARR_CNT(skip_multiples); j++) {
ii += skip_multiples[j];
elim_factor(ii);
}
if (ii > vv)
break;
}
/* display the batch of generated prime numbers */
for (n = 0; n < buffer_size && count < number; n++) {
if (prime[n]) {
ii = start_value + n;
if (ii > end_value) {
return;
}
if (ii > 1.0) {
if (pal_flag && !test_pal(ii, pal_base))
continue;
if (twin_flag) {
if ((last_prime + 2.0L) == ii) {
printf("%.0Lf %.0Lf\n", last_prime, ii);
count += 2;
}
} else {
printf("%.0Lf\n", ii);
count++;
}
last_prime = ii;
}
}
}
}
}
/*
* Parse a space or null terminated ASCII number in the string pointed to by "cp" and
* return true with a floating point long double value in "*dp" if a valid integer,
* otherwise display an error message and return false.
*/
int
get_long_double_int(char *cp, long double *dp)
{
char *cp1;
*dp = strtold(cp, &cp1);
if (cp == cp1) {
fprintf(stderr, "Invalid number.\n");
return false;
}
switch (*cp1) {
case '\0':
case '\r':
case '\n':
break;
default:
if (isspace(*cp1)) {
break;
}
fprintf(stderr, "Invalid number.\n");
return false;
}
if (*dp > MAX_K_INTEGER) {
fprintf(stderr, "Number is too large.\n");
return false;
}
if (*dp < 0.0 || fmodl(*dp, 1.0L) != 0.0) {
fprintf(stderr, "Number must be a positive integer.\n");
return false;
}
return true;
}
/*
* Return true if "d" is a palindrome base "base".
*/
int
test_pal(long double d, long double base)
{
#define MAX_DIGITS 1000
int i, j;
int digits[MAX_DIGITS];
/* build the array of digits[] */
for (i = 0; d >= 1.0; i++) {
assert(i < MAX_DIGITS);
digits[i] = (int) fmodl(d, base);
d /= base;
}
/* compare the array of digits[] end to end */
for (j = 0, i--; j < i; j++, i--) {
if (digits[i] != digits[j])
return false;
}
return true;
}
void
usage()
{
fprintf(stderr, "\nUsage: %s [start [stop]] [\"twin\"] [\"pal\" [base]]\n\n", prog_name);
fprintf(stderr, "Generate batches of consecutive prime numbers up to 19 digits.\n");
fprintf(stderr, "If \"twin\" is specified, output only twin primes.\n");
fprintf(stderr, "If \"pal\" is specified, output only palindromic primes.\n");
fprintf(stderr, "The palindromic base may be specified, the default is base 10.\n");
exit(1);
}
void
usage2()
{
fprintf(stderr, "\nUsage: %s [options] [start [stop]]\n\n", prog_name);
fprintf(stderr, "Generate batches of consecutive prime numbers up to 19 digits.\n");
fprintf(stderr, "Options:\n");
fprintf(stderr, " -t Output only twin primes.\n");
fprintf(stderr, " -p base Output only palindromic primes.\n");
exit(1);
}
syntax highlighted by Code2HTML, v. 0.9.1