Velocity Reviews > Re: PRIME1 (SPOJ)

# Re: PRIME1 (SPOJ)

Dann Corbit
Guest
Posts: n/a

 10-20-2010
In article <4cbe825e\$0\$10754\$(E-Mail Removed)>,
http://www.velocityreviews.com/forums/(E-Mail Removed) says...
>
> superpollo ha scritto:
> > superpollo ha scritto:
> >> ok, the problem is:
> >>
> >> https://www.spoj.pl/problems/PRIME1/
> >>
> >> so i submitted this (C99) code:

> >
> > would it be any different (speed-wise) if i used pointer syntax instead
> > of array syntax for primes[] (since they are equivalent) ?

>

Your algorithm fails for numbers larger than 31607 * 31607 = 999002449

As long as your inputs are less than 2^64th, you could use Baillie-PSW
which has been proven correct up through 2^64th.

A similar concept to what you are doing, nicely implented might be worth
study. Here is a nice solution to this problem that was written by
James Dow Allen:

/*
* prime.c -- Find primes in a range.
* Copyright (c) 2006 by James Dow Allen
*
* This program can be distributed or modified subject to
* the terms of the GNU General Public License (version 2
* the Free Software Foundation.
*
* Supports only smallish integers (up to 4 billion or so).
* (I have routines, requiring no other libraries, that handle
* primes up to 10 billion billion. I'd probably send
* them to anyone who says "Pretty please.")
*/

#include <stdlib.h>
#include <stdio.h>

typedef unsigned long int Integer; /* assume this is 32-bits */
typedef unsigned long long int Largenum; /* assume this is 64-bits */

/*
* Following number must be a multiple of 30,
* at least 900 (for functionality) but even
* larger (for performance).
*/
#define MTHRESH 720000

/*
* Actually we have several versions
* of this routine, depending on whether
* intermediate results are 32-bit, 64-bit, or larger.
* Here, just one version, perhaps good enough for 32-bit n
*/
int fermtst2(int a, Integer n)
{
Integer i, b, resu;

b = n - 1;
for (i = 1 << 31; b && !(i & b); i >>= 1)
;
resu = a;
for (i >>= 1; i; i >>= 1) {
resu = (resu * (Largenum)resu) % n;
if (i & b) {
resu = (a * (Largenum)resu) % n;
resu %= n;
}
}
return resu == 1;
}

int ispcandid(Integer x)
{
return
/* Caller must do 2, 3, 5 herself */
(x % 7) &&
(x % 11) && (x % 13) && (x % 17) && (x % 19) &&
(x % 23) && (x % 29) && (x % 31) && (x % 37) &&
(x % 41) && (x % 43) && (x % 47) && (x % 53) &&
(x % 59) && (x % 61) && (x % 67) && (x % 71) &&
(x % 73) && (x % 79) && (x % 83) && (x % 89) &&
(x % 97) && (x % 101) && (x % 103) && (x % 107) &&
(x % 109) && (x % 113) && (x % 127) && (x % 131) &&
(x % 137) && (x % 139) && (x % 149) && (x % 151) &&
(x % 157) && (x % 163) && (x % 167) && (x % 173) &&
(x % 179) && (x % 181) && (x % 191) && (x % 193) &&
(x % 197) && (x % 199) && (x % 211) && (x % 223) &&
(x % 227) && (x % 229) && (x % 233) && (x % 239) &&
(x % 241) && (x % 251) && (x % 257) && (x % 263) &&
/*
* Prior to here, try all primes; after here,
* just those that kill specific pseudoprimes.
* Note that these aren't costly: trying to
* bypass fermtst2() is always a partial win.
*/
(x % 271) && (x % 277) && (x % 281) && (x % 293) &&
(x % 307) && (x % 311) && (x % 331) && (x % 337) &&
(x % 349) && (x % 367) && (x % 373) && (x % 379) &&
(x % 401) && (x % 421) && (x % 433) && (x % 461) &&
(x % 487) && (x % 541) && (x % 547) && (x % 601) &&
(x % 613) && (x % 661) &&
/* Preceding perfect for x < 2 billion;
* following improves that to 2.7 billion.
*/
(x % 353) && (x % 443) && (x % 499) && (x % 727) &&
(x % 739) && (x % 1093) &&
(x % 34501) /* nasty: 2380603501 = 34501 69001 */ &&
1;
}

int fermtrial[] = {
2, 3, 5, 7, 11, 13, 17,
10111,
};

int isprime(Integer x)
{
int fix, f, fbnd;

/*
* Save some cycles when x is small.
* It relies on divisibility by 601, 661
* (and others) having been tested elsewhere.
* Otherwise we'd need
* fbnd = x < 721801 ? 5 : x < 42702661 ? 11 : 17;
*/
fbnd = x < 9006401 ? 5 : x < 42702661 ? 11 : 17;
for (fix = 0; (f = fermtrial[fix]) <= fbnd; fix++)
if (! fermtst2(f, x))
return 0;
return 1;
}

void emit(Integer x)
{
printf("%lu\n", x);
}

void cemit(Integer x)
{
if (ispcandid(x) && isprime(x))
emit(x);
}

int main(int argc, char **argv)
{
Integer x, last;
int resid, y, z;

if (argc == 2) {
x = last = (Integer) atof(argv[1]);
} else if (argc == 3) {
x = (Integer) atof(argv[1]);
last = (Integer) atof(argv[2]);
} else {
printf("Usage: primelist #a #b -- list primes from a to b\n");
printf(" Or: primelist #a -- show a if prime\n");
exit(1);
}
if (x < MTHRESH) {
if (x < 3) {
if (last >= 2)
emit((Integer)2);
x = 3;
} else {
x += !(x % 2);
}
/*
* We "should" do the following with a variation
* of (Duff's) sieve as below, but
* o quite fast anyway, for MTHRESH < 1 million.
* o for real speed we might use table for
* x < 100 million or so.
* o because range needs segmentation based on
* last and 30-mult, 3 loops would be needed.
* o most of the small odds are prime anyway,
* so avoiding 9,15,21,... isn't quite as
* big a win as one might think.
*/
for (; x < MTHRESH; x += 2) {
if (x > last)
exit(0);
for (y = 3; z = x / y, y <= z; y += 2)
if (z * y == x)
goto bigbreak;
emit(x);
bigbreak:
;
}
}
resid = x % 30;
x -= resid;
if (x + 30 <= last) {
/* Erastothenes' sieve is done a la Duff's Device! */
switch (resid) {
for (; x + 30 <= last; x += 30) {
case 0:
case 1:
cemit(x + 1);
case 2:
case 3:
case 4:
case 5:
case 6:
case 7:
cemit(x + 7);
case 8:
case 9:
case 10:
case 11:
cemit(x + 11);
case 12:
case 13:
cemit(x + 13);
case 14:
case 15:
case 16:
case 17:
cemit(x + 17);
case 18:
case 19:
cemit(x + 19);
case 20:
case 21:
case 22:
case 23:
cemit(x + 23);
case 24:
case 25:
case 26:
case 27:
case 28:
case 29:
cemit(x + 29);
}
}
} else
x += resid;
for (x += !(x % 2); x <= last; x += 2)
if (x % 3 && x % 5)
cemit(x);
exit(0);
}

Ben Bacarisse
Guest
Posts: n/a

 10-20-2010
Dann Corbit <(E-Mail Removed)> writes:
<snip>
> A similar concept to what you are doing, nicely implented might be worth
> study. Here is a nice solution to this problem that was written by
> James Dow Allen:

It's not quite portable and it might be hard for the OP to discover the
problem. The simplest fix is to change

1 << 31

to

1LU << 31

in fermtst2.

<snip>
--
Ben.

Dann Corbit
Guest
Posts: n/a

 10-20-2010
In article
<0.8a8325db8e725b8b496f.20101020154634BST.878w1sex (E-Mail Removed)>,
(E-Mail Removed) says...
>
> Dann Corbit <(E-Mail Removed)> writes:
> <snip>
> > A similar concept to what you are doing, nicely implented might be worth
> > study. Here is a nice solution to this problem that was written by
> > James Dow Allen:

>
> It's not quite portable and it might be hard for the OP to discover the
> problem. The simplest fix is to change
>
> 1 << 31
>
> to
>
> 1LU << 31
>
> in fermtst2.
>
> <snip>

Thanks for that correction.