/* Sieve of Xuedong Luo (Algorithm3). Return a reference to an array
 * containing all prime numbers less than or equal to search_to.
 *
 *    A practical sieve algorithm for finding prime numbers.
 *    ACM Volume 32 Issue 3, March 1989, Pages 344-346
 *    http://dl.acm.org/citation.cfm?doid=62065.62072
 *
 * Avoid all composites that have 2 or 3 as one of their prime factors
 * where i is odd.
 *
 *    { 0, 5, 7, 11, 13, ... 3i + 2, 3(i + 1) + 1, ..., N }
 *      0, 1, 2,  3,  4, ... list indices (0 is not used)
 */

SV* primes( long search_to )
{
    AV* av = newAV();

    if( search_to < 2 )
        return newRV_noinc( (SV*) av ); // Return an empty list ref.

    sieve_size_t i, j, q = (sieve_size_t) sqrt((double) search_to) / 3;
    sieve_size_t M = (sieve_size_t) search_to / 3;
    sieve_size_t c = 0, k = 1, t = 2, ij;

    // Allocate space. Set bits to 1. Unset bit 0.
    sieve_type primes( M + 2, 1 );
    primes[0] = 0;

    // Unset bits greater than search_to.
    if ( 3 * M + 2 > search_to + ((sieve_size_t)search_to & 1) )
        primes[M] = 0;
    if ( 3 * (M + 1) + 1 > search_to + ((sieve_size_t)search_to & 1) )
        primes[M + 1] = 0;

    // Clear composites.
    for ( i = 1; i <= q; i++ ) {
        k  = 3 - k, c = 4 * k * i + c, j = c;
        ij = 2 * i * (3 - k) + 1, t = 4 * k + t;
        if ( primes[i] ) {
            while ( j <= M ) {
                primes[j] = 0;
                j += ij, ij = t - ij;
            }
        }
    }

    // Gather primes.
    if( search_to >= 2 )
        av_push( av, newSVuv( 2UL ) );
    if( search_to >= 3 )
        av_push( av, newSVuv( 3UL ) );

    for ( i = 1; i <= M; i += 2 ) {
        if ( primes[i] )
          av_push( av, newSVuv(static_cast<unsigned long>(3 * i + 2)) );
        if ( primes[i + 1] )
          av_push( av, newSVuv(static_cast<unsigned long>(3 * (i + 1) + 1)) );
    }

    return newRV_noinc( (SV*) av );
}