/* Sieve of Eratosthenes.  Return a reference to an array containing all
 * prime numbers less than or equal to search_to.  Uses an optimized sieve
 * that requires one bit per odd from 0 .. n.  Evens aren't represented in the
 * sieve.  2 is just handled as a special case.
 */

SV* primes( long search_to )
{
    AV* av = newAV();
    if( search_to < 2 )
        return newRV_noinc( (SV*) av ); // Return an empty list ref.
    av_push( av, newSVuv( 2UL ) );
    // Allocate space for odd numbers (15 bits per 30 values)
    sieve_type primes( search_to/2 + 1, 0 );
    // Sieve over the odd numbers
    for( sieve_size_t i = 3; i * i <= search_to; i+=2 )
        if( ! primes[i/2] )
            for( sieve_size_t k = i*i; k <= search_to; k += 2*i)
                primes[k/2] = 1;
    // Add each prime to the list ref
    for( sieve_size_t i = 3; i <= search_to; i += 2 )
        if( ! primes[i/2] )
          av_push( av, newSVuv( static_cast<unsigned long>( i ) ) );
    return newRV_noinc( (SV*) av );
}