GRFreire/mathomatic
1
0
Fork 0
mirror of https://github.com/mfillpot/mathomatic.git synced 2026-01-09 04:59:37 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity
master
mathomatic/primes/matho-primes.c

597 lines
16 KiB
C
Raw Permalink Blame History

/*
* 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) 1987-2012 George Gesslein II.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
The chief copyright holder can be contacted at gesslein@mathomatic.org, or
George Gesslein II, P.O. Box 224, Lansing, NY 14882-0224 USA.
*/
/*
* Changes:
*
* 11/22/05 - converted everything to long doubles. Now uses C99 long double functions.
* 3/25/06 - made primes buffer variable size.
* 3/30/08 - Allow long double to be aliased to double when long double isn't supported.
* 2/11/09 - Cleanup calculation of number of decimal digits and max_integer.
* 9/12/10 - General cleanup and added error message for when the requested number of primes to display is not reached.
* 10/25/10 - Added -c and -h options.
* 10/26/10 - Using usage2() instead of usage() most of the time, now.
* 10/28/10 - Fixed to work with 16-bit integers.
* 1/16/11 - Allow to use doubles instead of long doubles with USE_DOUBLES define, for systems that don't fully support long doubles.
* 5/10/11 - Added -u (unbuffered output) option, for real-time output.
* 11/13/11 - Use EXIT_SUCCESS and EXIT_FAILURE macros.
* 11/30/11 - Added -m option.
* 12/31/11 - Added a compile-time warning when the number of digits of precision of long doubles is less than 18.
* 07/11/12 - Added -v option to display program name with version number and then exit successfully.
* 08/09/12 - Allow "matho-primes all" for endless output of consecutive primes.
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include <errno.h>
#include <math.h>
#include <float.h>
#include <assert.h>
#if !NO_GETOPT_H
#include <getopt.h>
#endif
#define VERSION "1.4" /* The current version number of this primes program. */
#define true 1
#define false 0
#ifndef max
#define max(a, b) (((a) > (b)) ? (a) : (b)) /* return the maximum of two values */
#endif
#ifndef min
#define min(a, b) (((a) < (b)) ? (a) : (b)) /* return the minimum of two values */
#endif
#define ARR_CNT(a) ((int) (sizeof(a)/sizeof(a[0]))) /* returns the number of elements in an array */
#if !USE_DOUBLES
#ifndef LDBL_DIG
#define USE_DOUBLES 1
#warning "long doubles not supported by this C compiler."
#elif LDBL_DIG < 18
#warning "long double number of digits precision is less than 18, probably meaning matho-primes is double precision only."
#warning "Due to the above floating point warning, more testing or compiling with -DUSE_DOUBLES is recommended."
#endif
#endif
#if USE_DOUBLES
typedef double double_type;
#warning "Using only double floats instead of long double floats, as requested."
#else
typedef long double double_type;
#endif
/* Maximum memory usage in bytes; can be set to any size. */
#ifndef BUFFER_SIZE
#define BUFFER_SIZE 2000000
#endif
#if BUFFER_SIZE >= (INT_MAX / 2) || BUFFER_SIZE < 100
#warning BUFFER_SIZE out of range, using default.
#undef BUFFER_SIZE
#define BUFFER_SIZE min((INT_MAX / 2), 2000000)
#endif
void generate_primes(void);
int test_pal(double_type d, double_type base);
void usage(int ev);
void usage2(int ev);
int get_double_type_int(char *cp, double_type *dp);
#if !USE_DOUBLES
long double powl(), ceill(), sqrtl(), fmodl(), strtold();
#endif
double_type max_integer; /* largest value of a double_type integral value */
double_type start_value; /* where to start finding primes */
double_type number; /* number of prime lines to display */
int count_requested; /* true if the number of primes to display is set by "number" above */
double_type default_number = 20; /* default number of primes to display */
double_type end_value; /* where to stop finding primes */
double_type 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) */
int pal_flag, twin_flag;
double_type pal_base = 10; /* The palindrome base, if displaying palindromic primes. */
char *prime; /* The boolean sieve array (buffer) for the current window of numbers being tested for primality; */
/* each char (byte) contains true or false, true if prime. */
int buffer_size; /* Number of bytes for variable size prime[] buffer above. */
char *prog_name = "matho-primes";
int
main(int argc, char *argv[])
{
#if NO_GETOPT_H /* if no getopt.h is available */
extern char *optarg; /* set by getopt(3) */
extern int optind;
#endif
int i;
char buf[1000];
char *cp = NULL;
double new_size = 0;
buffer_size = BUFFER_SIZE;
/* set the highest number this program will work with: */
#if USE_DOUBLES
max_integer = pow(10.0, (double) (DBL_DIG));
#else
max_integer = powl(10.0L, (long double) (LDBL_DIG));
#endif
while (max_integer == max_integer + 1.0) {
#if USE_DOUBLES
fprintf(stderr, "Warning: max_integer (%g) is too large; size of double = %u bytes.\n", max_integer, (unsigned) sizeof(double));
max_integer /= 10.0;
#else
fprintf(stderr, "Warning: max_integer (%Lg) is too large; size of long double = %u bytes.\n", max_integer, (unsigned) sizeof(long double));
max_integer /= 10.0L;
#endif
}
#if MINGW || __APPLE__
do_again:
#endif
start_value = -1.0;
end_value = max_integer;
number = 0;
count_requested = false;
/* process command line options: */
while ((i = getopt(argc, argv, "c:thuvp:m:")) >= 0) {
switch (i) {
case 'c':
count_requested = true;
if (optarg && !get_double_type_int(optarg, &number)) {
usage2(EXIT_FAILURE);
}
break;
case 'h':
usage2(EXIT_SUCCESS);
break;
case 't':
twin_flag = true;
break;
case 'p':
pal_flag = true;
if (optarg && !get_double_type_int(optarg, &pal_base)) {
usage2(EXIT_FAILURE);
}
break;
case 'm':
if (optarg)
new_size = strtod(optarg, &cp) * BUFFER_SIZE;
if (optarg == NULL || cp == NULL || *cp || new_size < 100 || new_size >= (INT_MAX / 2)) {
fprintf(stderr, "%s: Invalid memory size multiplier specified.\n", prog_name);
exit(EXIT_FAILURE);
}
buffer_size = (int) new_size;
fprintf(stderr, "%s: Window size = %d bytes.\n", prog_name, buffer_size);
break;
case 'u':
setbuf(stdout, NULL); /* make output unbuffered */
setbuf(stderr, NULL);
break;
case 'v':
printf("%s version %s\n", prog_name, VERSION);
exit(0);
default:
usage2(EXIT_FAILURE);
}
}
for (; argc > optind; optind++) {
if (strcasecmp(argv[optind], "all") == 0) {
if (start_value < 0)
start_value = 0;
count_requested = true;
number = max_integer;
continue;
}
if (strcasecmp(argv[optind], "twin") == 0) {
twin_flag = true;
continue;
}
break;
}
if (argc > optind && isdigit(argv[optind][0])) {
if (get_double_type_int(argv[optind], &start_value)) {
optind++;
} else {
usage2(EXIT_FAILURE);
}
if (argc > optind && isdigit(argv[optind][0])) {
if (get_double_type_int(argv[optind], &end_value)) {
if (end_value < start_value) {
fprintf(stderr, "End value is less than start value.\n");
usage2(EXIT_FAILURE);
}
optind++;
if (number == 0)
number = max_integer;
} else {
usage2(EXIT_FAILURE);
}
}
}
for (; argc > optind; optind++) {
if (strcasecmp(argv[optind], "all") == 0) {
if (start_value < 0)
start_value = 0;
count_requested = true;
number = max_integer;
continue;
}
if (strcasecmp(argv[optind], "twin") == 0) {
twin_flag = true;
continue;
}
break;
}
if (argc > optind) {
if (strncasecmp(argv[optind], "pal", 3) == 0) {
pal_flag = true;
optind++;
} else {
fprintf(stderr, "Unrecognized argument: \"%s\".\n", argv[optind]);
usage(EXIT_FAILURE);
}
if (argc > optind && isdigit(argv[optind][0])) {
if (!get_double_type_int(argv[optind], &pal_base)) {
usage(EXIT_FAILURE);
}
optind++;
}
}
for (; argc > optind; optind++) {
if (strcasecmp(argv[optind], "all") == 0) {
if (start_value < 0)
start_value = 0;
count_requested = true;
number = max_integer;
continue;
}
if (strcasecmp(argv[optind], "twin") == 0) {
twin_flag = true;
continue;
}
break;
}
if (argc > optind) {
fprintf(stderr, "Unrecognized argument: \"%s\".\n", argv[optind]);
usage(EXIT_FAILURE);
}
if (pal_base < 2 || pal_base >= INT_MAX) {
fprintf(stderr, "Palindrome number base must be >= 2.\n");
usage(EXIT_FAILURE);
}
if (start_value < 0.0) {
get1:
fprintf(stderr, "Enter number to start finding consecutive primes at (0): ");
fflush(NULL);
if (fgets(buf, sizeof(buf), stdin) == NULL)
exit(EXIT_SUCCESS);
switch (buf[0]) {
case '\0':
case '\n':
case '\r':
start_value = 0;
break;
default:
if (!get_double_type_int(buf, &start_value)) {
goto get1;
}
}
}
if (number == 0) {
get2:
fprintf(stderr, "Enter number of%s%s primes to output (0 to end)",
pal_flag ? " palindromic" : " consecutive", twin_flag ? " twin" : "");
#if USE_DOUBLES
fprintf(stderr, " (%g): ", default_number);
#else
fprintf(stderr, " (%Lg): ", default_number);
#endif
fflush(NULL);
if (fgets(buf, sizeof(buf), stdin) == NULL)
exit(EXIT_SUCCESS);
switch (buf[0]) {
case '\0':
case '\n':
case '\r':
number = default_number;
break;
default:
if (!get_double_type_int(buf, &number)) {
goto get2;
}
}
count_requested = true;
}
if (prime == NULL) {
/* allocate the prime[] buffer: */
prime = (char *) malloc(buffer_size);
}
if (prime == NULL) {
fprintf(stderr, "%s: Not enough memory for buffer_size = %d.\n", prog_name, buffer_size);
exit(EXIT_FAILURE);
}
generate_primes();
#if MINGW || __APPLE__
fflush(NULL);
if (argc <= 1 && number > 0) {
goto do_again;
}
#endif
exit(EXIT_SUCCESS);
}
/*
* Eliminate all multiples of "arg" from the sieve array by setting their entries to false.
* The sieve array "prime[]" is the prime truth values of a consecutive
* batch of numbers starting at "start_value".
*
* When all multiplies of all primes up to the square root
* of the highest value represented in sieve array are registered,
* the sieve array will tell whether the numbers in it are prime or not.
* This is how the Sieve of Eratosthenes works.
*/
void
elim_factor(double_type arg)
{
double_type d;
int i, j;
#if USE_DOUBLES
d = ceil(start_value / arg);
#else
d = ceill(start_value / arg);
#endif
if (d < 2.0)
d = 2.0;
d *= arg;
d -= start_value;
if (d >= buffer_size)
return;
i = (int) d;
if (i >= buffer_size)
return;
assert(i >= 0);
if (arg >= buffer_size) {
prime[i] = false;
} else {
j = (int) arg;
assert(j > 0);
for (; i < buffer_size; i += j) {
prime[i] = false;
}
}
}
/*
* Generate and display at most "number" consecutive prime numbers,
* between "start_value" and "end_value".
*/
void
generate_primes(void)
{
int n, j;
double_type count, d, ii, sqrt_value, last_prime = -3.0;
for (count = 0; count < number; start_value += buffer_size) {
if (start_value > end_value) {
goto check_return;
}
/* generate a batch of consecutive primes with the prime number sieve */
memset(prime, true, buffer_size); /* set the prime array to all true */
#if USE_DOUBLES
sqrt_value = 1.0 + sqrt(min(start_value + (double) buffer_size, end_value));
#else
sqrt_value = 1.0 + sqrtl(min(start_value + (long double) buffer_size, end_value));
#endif
elim_factor((double_type) 2.0);
elim_factor((double_type) 3.0);
elim_factor((double_type) 5.0);
elim_factor((double_type) 7.0);
d = 1.0;
while (d <= sqrt_value) {
for (j = 0; j < ARR_CNT(skip_multiples); j++) {
d += skip_multiples[j];
elim_factor(d);
}
}
/* display the requested part of the batch of generated prime numbers */
for (n = 0; n < buffer_size && count < number; n++) {
if (prime[n]) { /* if prime number */
ii = start_value + n;
if (ii > end_value) {
goto check_return;
}
if (ii > 1.0) {
if (pal_flag && !test_pal(ii, pal_base))
continue;
if (twin_flag) {
if ((last_prime + 2.0) == ii) {
#if USE_DOUBLES
printf("%.0f %.0f\n", last_prime, ii);
#else
printf("%.0Lf %.0Lf\n", last_prime, ii);
#endif
count++;
}
} else {
#if USE_DOUBLES
printf("%.0f\n", ii);
#else
printf("%.0Lf\n", ii);
#endif
count++;
}
last_prime = ii;
}
}
}
}
check_return:
if (count_requested && count < number) {
fprintf(stderr, "%s: Number of primes requested not reached.\n", prog_name);
exit(EXIT_FAILURE);
}
}
/*
* Parse a space or null terminated ASCII number in the string pointed to by cp.
*
* Return true with a floating point double_type value in *dp
* if a valid positive integer or zero,
* otherwise display an error message and return false.
*/
int
get_double_type_int(char *cp, double_type *dp)
{
char *cp1;
if (cp == NULL || dp == NULL) {
return false;
}
errno = 0;
#if USE_DOUBLES
*dp = strtod(cp, &cp1);
#else
*dp = strtold(cp, &cp1);
#endif
if (errno) {
perror(NULL);
return false;
}
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_integer) {
#if USE_DOUBLES
fprintf(stderr, "Number is too large, maximum is %g.\n", max_integer);
#else
fprintf(stderr, "Number is too large, maximum is %Lg.\n", max_integer);
#endif
return false;
}
#if USE_DOUBLES
if (*dp < 0.0 || fmod(*dp, 1.0) != 0.0) {
#else
if (*dp < 0.0 || fmodl(*dp, 1.0L) != 0.0) {
#endif
fprintf(stderr, "Number must be a positive integer or zero.\n");
return false;
}
return true;
}
/*
* Return true if "d" is a palindromic number, base "base".
*/
int
test_pal(double_type d, double_type 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);
#if USE_DOUBLES
digits[i] = (int) (fmod(d, base));
#else
digits[i] = (int) (fmodl(d, base));
#endif
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(int ev)
{
printf("Prime number generator version %s\n", VERSION);
printf("Usage: %s [start [stop] or \"all\"] [\"twin\"] [\"pal\" [base]]\n\n", prog_name);
#if USE_DOUBLES
printf("Generate consecutive prime numbers from start to stop, up to %g.\n", max_integer);
#else
printf("Generate consecutive prime numbers from start to stop, up to %Lg.\n", max_integer);
#endif
printf("If \"twin\" is specified, output only twin primes.\n");
printf("If \"pal\" is specified, output only palindromic primes.\n");
printf("The palindrome number base may be specified, the default is base 10.\n");
exit(ev);
}
void
usage2(int ev)
{
printf("Prime number generator version %s\n", VERSION);
printf("Usage: %s [options] [start [stop]]\n\n", prog_name);
#if USE_DOUBLES
printf("Generate consecutive prime numbers from start to stop, up to %g.\n", max_integer);
#else
printf("Generate consecutive prime numbers from start to stop, up to %Lg.\n", max_integer);
#endif
printf("Options:\n");
printf(" -c count Count lines of primes, stop when count reached.\n");
printf(" -h Display this help and exit.\n");
printf(" -m number Specify a memory size multiplier.\n");
printf(" -p base Output only palindromic primes.\n");
printf(" -t Output only twin primes.\n");
printf(" -u Set all output to be unbuffered.\n");
printf(" -v Display version number, then exit successfully.\n");
exit(ev);
}
Reference in New Issue View Git Blame Copy Permalink