Logo Search packages:      
Sourcecode: mathomatic version File versions

matho-primes.c

/*
 * Generate consecutive prime numbers using a modified sieve method
 * that doesn't use much memory.
 * Works on up to 14 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.
 */

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <assert.h>

#define     true  1
#define     false 0

/* memory usage in bytes; can be set to any size; the larger, the faster */
#define MAX 200000

void elim_factor ( double arg );
int test_pal ( double d , double base );
void usage ( void );
int get_double_int(char *cp, double *dp);

double nn, vv;
double huge_value = 1.0e14;
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
};
char  prime[MAX];
char  *prog_name;

main(int argc, char *argv[])
{
      int         n;
      double            count, number;
      int         j;
      double            ii;
      double            end_value;
      char        buf[50];
      int         pal_flag;
      double            pal_base;
      int         start;

      prog_name = argv[0];
      nn = -1.0;
      end_value = huge_value;
      number = 0;
      pal_flag = false;
      pal_base = 10;
      start = 1;
      if (argc > start && isdigit(argv[start][0])) {
            if (get_double_int(argv[start], &nn)) {
                  start++;
            } else {
                  usage();
            }
            if (argc > start && isdigit(argv[start][0])) {
                  if (get_double_int(argv[start], &end_value)) {
                        start++;
                        number = huge_value;
                  } else {
                        usage();
                  }
            }
      }
      if (argc > start) {
            if (strncmp(argv[start], "pal", 3) == 0) {
                  pal_flag = true;
                  start++;
            } else {
                  usage();
            }
      }
      if (argc > start) {
            if (sscanf(argv[start], "%lf", &pal_base) != 1) {
                  usage();
            }
            if (pal_base < 2 || pal_base > 256) {
                  fprintf(stderr, "Palindromic base range is 2 to 256.\n");
                  exit(1);
            }
            start++;
      }
      if (argc > start) {
            usage();
      }
      if (nn < 0.0) {
get1:
            fprintf(stderr, "Enter number (up to 14 digits) to start finding primes at: ");
            if (fgets(buf, sizeof(buf), stdin) == NULL)
                  exit(1);
            if (!get_double_int(buf, &nn)) {
                  goto get1;
            }
      }
      if (number == 0) {
get2:
            fprintf(stderr, "Enter number of primes to find (default = 20): ");
            if (fgets(buf, sizeof(buf), stdin) == NULL)
                  exit(1);
            if (buf[0] == '\0' || buf[0] == '\n') {
                  number = 20;
            } else {
                  if (!get_double_int(buf, &number)) {
                        goto get2;
                  }
            }
      }
      if (pal_flag)
            fprintf(stderr, "Finding base %g palindromic primes starting at %.0f:\n", pal_base, nn);
      count = 0;
      for (; count < number; nn += MAX) {
            if (nn > end_value) {
                  fprintf(stderr, "Maximum reached.  %.0f primes found.\n", count);
                  exit(0);
            }
            vv = 1.0 + sqrt(nn + (double) MAX);
            for (n = 0; n < MAX; n++)
                  prime[n] = true;
            elim_factor(2.0);
            elim_factor(3.0);
            elim_factor(5.0);
            elim_factor(7.0);
            ii = 1.0;
            for (;;) {
                  for (j = 0; j < 48; j++) {
                        ii += sq[j];
                        elim_factor(ii);
                  }
                  if (ii > vv)
                        break;
            }
            for (n = 0; n < MAX && count < number; n++) {
                  if (prime[n]) {
                        ii = nn + n;
                        if (ii > end_value) {
                              fprintf(stderr, "Maximum reached.  %.0f primes found.\n", count);
                              exit(0);
                        }
                        if (ii != 0.0 && ii != 1.0) {
                              if (pal_flag && !test_pal(ii, pal_base))
                                    continue;
                              count++;
                              printf("%.0f\n", ii);
                        }
                  }
            }
      }
      exit(0);
}

int
get_double_int(cp, dp)
char  *cp;
double      *dp;
{
      if (sscanf(cp, "%lf", dp) != 1) {
            fprintf(stderr, "That is not a number!\n");
            return false;
      }
      if (*dp > huge_value) {
            fprintf(stderr, "Number is too large!\n");
            return false;
      }
      if (*dp < 0.0 || fmod(*dp, 1.0) != 0.0) {
            fprintf(stderr, "Number must be a positive integer!\n");
            return false;
      }
      return true;
}

/*
 * Eliminate all multiples of "arg" from the sieve array.
 */
void
elim_factor(double arg)
{
      double            d;
      unsigned int      i, j;

      d = ceil(nn / arg);
      if (d < 2.0)
            d = 2.0;
      d *= arg;
      d -= nn;
      if (d >= MAX)
            return;
      assert(d >= 0.0);
      i = d;
      if (arg >= MAX) {
            prime[i] = false;
      } else {
            j = arg;
            for (; i < MAX; i += j) {
                  prime[i] = false;
            }
      }
}

/*
 * Return true if "d" is a palindrome base "base".
 */
int
test_pal(double d, double base)
{
#define     MAX_DIGITS  100

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

      d1 = d;
      for (i = 0; d1 >= 1.0; i++) {
            assert(i < MAX_DIGITS);
            digits[i] = fmod(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 14 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(2);
}

Generated by  Doxygen 1.6.0   Back to index