Logo Search packages:      
Sourcecode: mathomatic version File versions  Download package

matho-primes.c

/*
 * Generate consecutive prime numbers using a modified sieve algorithm
 * that doesn't use much memory by using a windowing sieve buffer.
 * Works on up to 19 digit numbers.
 *
 * This is C or C++ program.
 * Copyright (C) 2004 George Gesslein II.
 *
 * Usage: matho-primes [start [stop]] ["pal" [base]]
 *
 * If "pal" is specified, display only palindromic primes of base "base",
 * default base is 10.
 *
 * Changes:
 *
 * 11/22/05 - converted everything to long doubles.
 */

#include <stdio.h>
#include <stdlib.h>
#if   UNIX
#include <libgen.h>
#endif
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include <errno.h>
#include <math.h>
#include <assert.h>

#define     true  1
#define     false 0

#define     MAX_K_INTEGER     1.0e19L           /* 19 digits for long doubles */

/* memory usage in bytes; can be set to any size; the larger, the faster */
#define MAX 500000
#if   MAX >= (INT_MAX / 2)
#error      MAX too big!
#endif

void  generate_primes(void);
int   test_pal(long double d, long double base);
void  usage(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 sq[] = {
      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
};

int         pal_flag;
long double pal_base = 10;

char  prime[MAX]; /* the boolean sieve array for the current window of numbers */

char  *prog_name = "matho-primes";

int
main(int argc, char *argv[])
{
      char        buf[100];
      int         start;

#if   UNIX
        prog_name = strdup(basename(argv[0]));
#endif
      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;
      start = 1;
      if (argc > start && isdigit(argv[start][0])) {
            if (get_long_double_int(argv[start], &start_value)) {
                  start++;
            } else {
                  usage();
            }
            if (argc > start && isdigit(argv[start][0])) {
                  if (get_long_double_int(argv[start], &end_value)) {
                        if (end_value < start_value) {
                              fprintf(stderr, "End value is less than start value.\n");
                              usage();
                        }
                        start++;
                        number = MAX_K_INTEGER;
                  } else {
                        usage();
                  }
            }
      }
      if (argc > start) {
            if (strncmp(argv[start], "pal", 3) == 0) {
                  pal_flag = true;
                  start++;
            } else {
                  usage();
            }
      }
      if (argc > start) {
            if (!get_long_double_int(argv[start], &pal_base)) {
                  usage();
            }
            if (pal_base < 2 || pal_base > 256) {
                  fprintf(stderr, "Palindromic base range is 2 to 256.\n");
                  usage();
            }
            start++;
      }
      if (argc > start) {
            usage();
      }
      if (start_value < 0.0) {
get1:
            fprintf(stderr, "Enter number (up to 19 digits) to start finding primes at: ");
            if (fgets(buf, sizeof(buf), stdin) == NULL)
                  exit(1);
            if (!get_long_double_int(buf, &start_value)) {
                  goto get1;
            }
      }
      if (number == 0) {
get2:
            fprintf(stderr, "Enter number of primes to output (default = 20): ");
            if (fgets(buf, sizeof(buf), stdin) == NULL)
                  exit(1);
            switch (buf[0]) {
            case '\0':
            case '\n':
            case '\r':
                  number = 20;
                  break;
            default:
                  if (!get_long_double_int(buf, &number)) {
                        goto get2;
                  }
            }
      }
      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 >= MAX)
            return;
      i = (int) d;
      assert(i >= 0);
      if (arg >= MAX) {
            prime[i] = false;
      } else {
            j = (int) arg;
            for (; i < MAX; i += j) {
                  prime[i] = false;
            }
      }
}

/*
 * Generate and display consecutive prime numbers.
 */
void
generate_primes()
{
      int         n;
      int         j;
      long double count, ii, vv;

      for (count = 0; count < number; start_value += MAX) {
            if (start_value > end_value) {
                  return;
            }
            vv = 1.0 + sqrtl(start_value + (long double) MAX);
            for (n = 0; n < MAX; n++)
                  prime[n] = true;
            elim_factor(2.0L);
            elim_factor(3.0L);
            elim_factor(5.0L);
            elim_factor(7.0L);
            ii = 1.0;
            for (;;) {
                  for (j = 0; j < 48; j++) {
                        ii += sq[j];
                        elim_factor(ii);
                  }
                  if (ii > vv)
                        break;
            }
            /* display the batch of generated prime numbers */
            for (n = 0; n < MAX && count < number; n++) {
                  if (prime[n]) {
                        ii = start_value + n;
                        if (ii > end_value) {
                              return;
                        }
                        if (ii != 0.0 && ii != 1.0) {
                              if (pal_flag && !test_pal(ii, pal_base))
                                    continue;
                              count++;
                              printf("%.0Lf\n", ii);
                        }
                  }
            }
      }
}

/*
 * Parse an 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.
 */
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:
            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  100

      int         i, j;
      long double d1;
      int         digits[MAX_DIGITS];

      d1 = d;
      for (i = 0; d1 >= 1.0; i++) {
            assert(i < MAX_DIGITS);
            digits[i] = (int) fmodl(d1, base);
            d1 /= base;
      }
      if (i <= 1)
            return false;
      for (j = 0, i--; j < i; j++, i--) {
            if (digits[i] != digits[j])
                  return false;
      }
      return true;
}

void
usage()
{
      fprintf(stderr, "Usage: %s [start [stop]] [\"pal\" [base]]\n\n", prog_name);
      fprintf(stderr, "Generate consecutive prime numbers up to 19 digits.\n");
      fprintf(stderr, "If \"pal\" is specified, display only palindromic primes.\n");
      fprintf(stderr, "The palindromic base may be specified, the default is base 10.\n");
      exit(1);
}

Generated by  Doxygen 1.6.0   Back to index