/* Gcc 4.1.2 + GMP
List,count primes, sum primes.
Origional Written by tege while on holiday in Rodupp, August 2001.
Between 10 and 500 times faster than previous program.
Modified code by Cino hilliard Oct 2007.
Original program was primes.c in the GMP 4.2.2 distribution
C:\gmp_build\gmp-4.2.2\demos\primes.c after building GMP on your system.
See the instructions at the end of the code.
Cino Hilliard add:
*Date and timing and timer output.
*Change total_primes to unsigned long long
*chr64() macro to printf unsigned long long as a string
*Instructions for installing GMP 4.2.2 Mingw on a P4 Windows XP.
Cino Hilliard delete:
* Debug stuff
* command line -c and -g tests.
Cino Comment:
*This program is twice as fast as the sieve program I've been using.
*Summing primes < 10^13 on p4 2.53 ghz 2 gigs ram takes 40 hours vs 76
*hours on Achim Flamkamp's sieve that I modified to sum primes.
Copyright 2001, 2002, 2006 Free Software Foundation, Inc.
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any later
version.
This program 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 General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see http://www.gnu.org/licenses/.
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <windows.h>
#define timer GetTickCount()/1000.0
#define ld long double;
#define chr64(a) _ui64toa(a,str,10)
#include "gmp.h"
/* IDEAS:
* Do not fill primes[] with real primes when the range [fr,to] is small,
when fr,to are relatively large. Fill primes[] with odd numbers instead.
[Probably a bad idea, since the primes[] array would become very large.]
* Separate small primes and large primes when sieving. Either the Montgomery
way (i.e., having a large array a multiple of L1 cache size), or just
separate loops for primes <= S and primes > S. The latter primes do not
require an inner loop, since they will touch the sieving array at most once.
* Pre-fill sieving array with an appropriately aligned ...00100100... pattern,
then omit 3 from primes array. (May require similar special handling of 3
as we now have for 2.)
* A large SIEVE_LIMIT currently implies very large memory usage, mainly due
to the sieving array in make_primelist, but also because of the primes[]
array. We might want to stage the program, using sieve_region/find_primes
to build primes[]. Make report() a function pointer, as part of achieving
this.
* Store primes[] as two arrays, one array with primes represented as delta
values using just 8 bits (if gaps are too big, store bogus primes!)
and one array with "rem" values. The latter needs 32-bit values.
* A new entry point, mpz_probab_prime_likely_p, would be useful.
* Improve command line syntax and versatility. "primes -f FROM -t TO",
allow either to be omitted for open interval. (But disallow
"primes -c -f FROM" since that would be infinity.) Allow printing a
limited *number* of primes using syntax like "primes -f FROM -n NUMBER".
* When looking for maxgaps, we should not perform any primality testing until
we find possible record gaps. Should speed up the searches tremendously.
*/
/* mpf_out_raw and mpf_inp_raw function by jim white */
int mpf_out_raw (FILE *f,mpf_t X) {
int expt;
mpz_t Z;
size_t nz;
expt = X->_mp_exp;
fwrite(&expt, sizeof(int), 1, f);
nz = X->_mp_size;
Z->_mp_alloc = nz;
Z->_mp_size = nz;
Z->_mp_d = X->_mp_d;
return (mpz_out_raw(f, Z) + sizeof(int));
}
void mpf_inp_raw (FILE *f, mpf_t X) {
int expt;
mpz_t Z;
size_t nz;
mpz_init (Z);
fread(&expt, sizeof(int), 1, f);
mpz_inp_raw (Z, f);
mpf_set_z (X, Z);
X->_mp_exp = expt;
mpz_clear (Z);
}
static float t;
struct primes
{
unsigned int prime;
int rem;
};
struct primes *primes;
unsigned long n_primes;
void find_primes __GMP_PROTO ((unsigned char *, mpz_t, unsigned long, mpz_t));
void sieve_region __GMP_PROTO ((unsigned char *, mpz_t, unsigned long));
void make_primelist __GMP_PROTO ((unsigned long));
int flag_print = 1;
int flag_count = 0;
int flag_maxgap = 0;
unsigned long maxgap = 0;
unsigned long long total_primes = 0;
mpz_t sum,r;
void
report (mpz_t prime)
{
total_primes += 1;
mpz_add(sum,sum,prime); //Sum the primes
}
char myfrom[24],myto[24];
int
main (int argc, char *argv[])
{
char *progname = argv[0],str[24];
strcpy(myfrom,argv[1]);
strcpy(myto,argv[2]);
mpz_t fr, to;
mpz_t fr2, to2;
unsigned long sieve_lim;
unsigned long est_n_primes;
unsigned char *s;
char *pEnd;
mpz_t tmp;
mpz_t siev_sqr_lim;
t=timer;
system("date/t");
system("time/t");
if (flag_count || flag_maxgap)
flag_print--; /* clear unless an explicit -p */
mpz_init (fr);
mpz_init (to);
mpz_init (fr2);
mpz_init (to2);
mpz_init (sum);
if (argc == 3)
{
mpz_set_str (fr, argv[1], 0);
if (argv[2][0] == '+')
{
mpz_set_str (to, argv[2] + 1, 0);
mpz_add (to, to, fr);
}
else
mpz_set_str (to, argv[2], 0);
}
else if (argc == 2)
{
mpz_set_ui (fr, 0);
mpz_set_str (to, argv[1], 0);
}
else
{
printf("Usage: sumprimesGmp start end");
exit (1);
}
mpz_set (fr2, fr);
if (mpz_cmp_ui (fr2, 3) < 0)
{
mpz_set_ui (fr2, 2);
report (fr2);
mpz_set_ui (fr2, 3);
}
mpz_setbit (fr2, 0); /* make odd */
mpz_sub_ui (to2, to, 1);
mpz_setbit (to2, 0); /* make odd */
mpz_init (tmp);
mpz_init (siev_sqr_lim);
mpz_sqrt (tmp, to2);
#define SIEVE_LIMIT 1000000000
if (mpz_cmp_ui (tmp, SIEVE_LIMIT) < 0)
{
sieve_lim = mpz_get_ui (tmp);
}
else
{
sieve_lim = SIEVE_LIMIT;
mpz_sub (tmp, to2, fr2);
if (mpz_cmp_ui (tmp, sieve_lim) < 0)
sieve_lim = mpz_get_ui (tmp); /* limit sieving for small ranges */
}
mpz_set_ui (siev_sqr_lim, sieve_lim + 1);
mpz_mul_ui (siev_sqr_lim, siev_sqr_lim, sieve_lim + 1);
est_n_primes = (size_t) (sieve_lim / log((double) sieve_lim) * 1.13) + 10;
primes = malloc (est_n_primes * sizeof primes[0]);
make_primelist (sieve_lim);
assert (est_n_primes >= n_primes);
#define S (1 << 17) /* FIXME: Figure out L1 cache size */
s = malloc (S/2);
while (mpz_cmp (fr2, to2) <= 0)
{
unsigned long rsize;
rsize = S;
mpz_add_ui (tmp, fr2, rsize);
if (mpz_cmp (tmp, to2) > 0)
{
mpz_sub (tmp, to2, fr2);
rsize = mpz_get_ui (tmp) + 2;
}
sieve_region (s, fr2, rsize);
find_primes (s, fr2, rsize / 2, siev_sqr_lim);
mpz_add_ui (fr2, fr2, S);
}
free (s);
t=timer-t;
printf("\n"); //Newline after primes are listed
printf ("Pi(interval) = %s\n", chr64(total_primes));
printf("Sum of primes = ");
mpz_out_str (stdout, 10, sum);
printf("\n");
printf ("Time in sec = %.06f\n",t);
return 0;
}
/* Find primes in region [fr,fr+rsize). Requires that fr is odd and that
rsize is even. The sieving array s should be aligned for "long int" and
have rsize/2 entries, rounded up to the nearest multiple of "long int". */
void
sieve_region (unsigned char *s, mpz_t fr, unsigned long rsize)
{
unsigned long ssize = rsize / 2;
unsigned long start, start2, prime;
unsigned long i;
mpz_t tmp;
mpz_init (tmp);
#if 0
/* initialize sieving array */
for (ii = 0; ii < (ssize + sizeof (long) - 1) / sizeof (long); ii++)
((long *) s) [ii] = ~0L;
#else
{
long k;
long *se = (long *) (s + ((ssize + sizeof (long) - 1) & -sizeof (long)));
for (k = -((ssize + sizeof (long) - 1) / sizeof (long)); k < 0; k++)
se[k] = ~0L;
}
#endif
for (i = 0; i < n_primes; i++)
{
prime = primes[i].prime;
if (primes[i].rem >= 0)
{
start2 = primes[i].rem;
}
else
{
mpz_set_ui (tmp, prime);
mpz_mul_ui (tmp, tmp, prime);
if (mpz_cmp (fr, tmp) <= 0)
{
mpz_sub (tmp, tmp, fr);
if (mpz_cmp_ui (tmp, 2 * ssize) > 0)
break; /* avoid overflow at next line, also speedup */
start = mpz_get_ui (tmp);
}
else
{
start = (prime - mpz_tdiv_ui (fr, prime)) % prime;
if (start % 2 != 0)
start += prime; /* adjust if even divisable */
}
start2 = start / 2;
}
#if 0
for (ii = start2; ii < ssize; ii += prime)
s[ii] = 0;
primes[i].rem = ii - ssize;
#else
{
long k;
unsigned char *se = s + ssize; /* point just beyond sieving range */
for (k = start2 - ssize; k < 0; k += prime)
se[k] = 0;
primes[i].rem = k;
}
#endif
}
mpz_clear (tmp);
}
/* Find primes in region [fr,fr+rsize), using the previously sieved s[]. */
void
find_primes (unsigned char *s, mpz_t fr, unsigned long ssize,
mpz_t siev_sqr_lim)
{
unsigned long j,ij;
mpz_t tmp;
mpz_init (tmp);
for (j = 0; j < (ssize + sizeof (long) - 1) / sizeof (long); j++)
{
if (((long *) s) [j] != 0)
{
for (ij = 0; ij < sizeof (long); ij++)
{
if (s[j * sizeof (long) + ij] != 0)
{
if (j * sizeof (long) + ij >= ssize)
goto out;
mpz_add_ui (tmp, fr, (j * sizeof (long) + ij) * 2);
if(mpz_cmp(tmp, siev_sqr_lim) < 0 ||
mpz_probab_prime_p (tmp,10))
{
report (tmp);
}
}
}
}
}
out:
mpz_clear (tmp);
}
/* Generate a list of primes and store in the global array primes[]. */
void
make_primelist (unsigned long maxprime)
{
#if 1
unsigned char *s;
unsigned long ssize = maxprime / 2;
unsigned long i, ii, j;
s = malloc (ssize);
memset (s, ~0, ssize);
for (i = 3; ; i += 2)
{
unsigned long isqr = i * i;
if (isqr >= maxprime)
break;
if (s[i * i / 2 - 1] == 0)
continue; /* only sieve with primes */
for (ii = i * i / 2 - 1; ii < ssize; ii += i)
s[ii] = 0;
}
n_primes = 0;
for (j = 0; j < ssize; j++)
{
if (s[j] != 0)
{
primes[n_primes].prime = j * 2 + 3;
primes[n_primes].rem = -1;
n_primes++;
}
}
/* FIXME: This should not be needed if fencepost errors were fixed... */
if (primes[n_primes - 1].prime > maxprime)
n_primes--;
free (s);
#else
unsigned long i;
n_primes = 0;
for (i = 3; i <= maxprime; i += 2)
{
if (i < 7 || (i % 3 != 0 && i % 5 != 0 && i % 7 != 0))
{
primes[n_primes].prime = i;
primes[n_primes].rem = -1;
n_primes++;
}
}
#endif
}
/***************************************************************************
Gcc MinGW (Building dynamic GMP library in Windows XP)
****************************************************************************
1. Install latest MinGW first, and then MSYS into the c drive.
http://www.mingw.org/download.shtml#hdr1
You can use another drive such as I did with f.
2. Create the folders c:\gmp_download, c:\gmp_build, c:\gmp_install
3. Download latest GMP from GNU MP to c:\gmp_download.
4. Run or click on "MSYS", unzip and untar GMP into c:\gmp_build using
following command:
cd /c/gmp_build
tar xzvf /c/gmp_download/gmp-x.x.x.tar.gz
5. Configure GMP for compilation: Dynamic Library
cd gmp-x.x.x
./configure --prefix=/c/gmp_install --disable-static --enable-shared
6. Build GMP:
make
This takes a while so be patient.
7. Install GMP header files and lib files:
make install
The libgmp-3.dll file winds up in c:\gmp_install\lib\bin. Copy this dll to
your working dir or to c:\mingw\bin. The gmp.h file winds up in
c:\gmp_install\include. Copy this into c:\mingw\include
The GMP Demos programs are found in C:\gmp_build\gmp-4.2.2\demos
Your batch file should be, say, gcgmp.bat
gcc %1.c -o %1.exe -O3 -s -mtune=pentium4 libgmp-3.dll
gcc -dumpversion
So c:\home\gcc sumprimes will compile and link your program to the GMP library.
Note: mingw/bin must be in the path. Porting to other computers is tricky if
they are not the same cpu and OS.
Cino hilliard
*/