Re: PRIME1 (SPOJ)

Discussion in 'C Programming' started by Dann Corbit, Oct 20, 2010.

  1. Dann Corbit

    Dann Corbit Guest

    In article <4cbe825e$0$10754$>,
    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) ?

    >
    > i am referring to the latest version: http://pastebin.org/320104


    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
    * or, at your option, any later version) as published by
    * 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);
    }
     
    Dann Corbit, Oct 20, 2010
    #1
    1. Advertising

  2. Dann Corbit <> 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.
     
    Ben Bacarisse, Oct 20, 2010
    #2
    1. Advertising

  3. Dann Corbit

    Dann Corbit Guest

    In article
    <>,
    says...
    >
    > Dann Corbit <> 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.
     
    Dann Corbit, Oct 20, 2010
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Shriphani
    Replies:
    2
    Views:
    572
    Ian Kelly
    Jun 1, 2008
  2. bert

    Re: PRIME1 (SPOJ)

    bert, Oct 18, 2010, in forum: C Programming
    Replies:
    4
    Views:
    456
    Michael Foukarakis
    Oct 20, 2010
  3. Ben Bacarisse

    Re: PRIME1 (SPOJ)

    Ben Bacarisse, Oct 19, 2010, in forum: C Programming
    Replies:
    2
    Views:
    451
    Ben Bacarisse
    Oct 19, 2010
  4. BartC

    Re: PRIME1 (SPOJ)

    BartC, Oct 19, 2010, in forum: C Programming
    Replies:
    5
    Views:
    532
    BartC
    Oct 19, 2010
  5. Ben Bacarisse

    Re: PRIME1 (SPOJ)

    Ben Bacarisse, Oct 19, 2010, in forum: C Programming
    Replies:
    8
    Views:
    562
    Dann Corbit
    Oct 20, 2010
Loading...

Share This Page