in reply to Number functions I have lying around
Note that the primes subroutine is quite inefficient and returns 1 as well, which is usually not considered prime.
Here's a faster one:
#! /usr/bin/perl
use warnings;
use strict;
use feature qw{ say };
sub primes {
my $n = shift;
return if $n < 2;
my @primes = (2);
for my $i (3 .. $n) {
my $sqrt = sqrt $i;
my $notprime;
for my $p (@primes) {
last if $p > $sqrt;
$notprime = 1, last if 0 == $i % $p;
}
push @primes, $i unless $notprime;
}
return @primes
}
use List::Util qw{ sum };
sub primes_la {
# Copy your code here.
}
use Test::More tests => 1;
is_deeply([1, primes(10000)], [primes_la(10000)], 'same');
use Benchmark qw{ cmpthese };
cmpthese(10,
{ ch => 'primes(10000)',
la => 'primes_la(10000)',
});
__END__
1..1
ok 1  same
s/iter la ch
la 1.35  99%
ch 1.06e02 12662% 
Re^2: Number functions I have lying around
by Discipulus (Monsignor) on Mar 31, 2015 at 09:23 UTC

ack! ouch! choroba now i need to replace the code for primality check taken from this node used in my TkTartaglia. I think your code is worth to put in the previously mentioned thread.
Rate la di ch
la 0.325/s  99% 99%
di 25.0/s 7572%  49%
ch 48.6/s 14822% 95% 
L*
There are no rules, there are no thumbs..
Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.
 [reply] [d/l] 

use strict;
use warnings;
use Benchmark qw{ cmpthese };
use Test::More tests => 1;
is_deeply(
[ primes( 10000 ) ],
[ primes_jg( 10000 ) ],
'same'
);
cmpthese(
10,
{
ch => 'primes( 10000 )',
jg => 'primes_jg( 10000 )',
}
);
sub primes_jg
{
my $limit = shift;
my $sqrtLimit = sqrt $limit;
my $sieve = q{};
vec( $sieve, 0, 1 ) = 1;
vec( $sieve, 1, 1 ) = 1;
vec( $sieve, $limit, 1 ) = 0;
my @primes;
my $reached = 1;
while( $reached < $sqrtLimit )
{
my $prime = $reached + 1;
++ $prime while vec( $sieve, $prime, 1 );
push @primes, $prime;
my $fill = 2 * $prime;
while( $fill <= $limit )
{
vec( $sieve, $fill, 1 ) = 1;
$fill += $prime;
}
$reached = $prime;
}
foreach my $value ( $reached + 1 .. $limit )
{
push @primes, $value unless vec( $sieve, $value, 1 );
}
return @primes;
}
sub primes {
my $n = shift;
return if $n < 2;
my @primes = (2);
for my $i (3 .. $n) {
my $sqrt = sqrt $i;
my $notprime;
for my $p (@primes) {
last if $p > $sqrt;
$notprime = 1, last if 0 == $i % $p;
}
push @primes, $i unless $notprime;
}
return @primes
}
1..1
ok 1  same
Rate ch jg
ch 71.0/s  25%
jg 94.7/s 33% 
I hope this is of interest.
 [reply] [d/l] [select] 

use Math::Prime::XS qw( sieve_primes );
my @primes = sieve_primes(10_000);
I'm to tired to write the benchmarks...
Best regards, Karl
«The Crux of the Biscuit is the Apostrophe»
 [reply] [d/l] 

choroba optimezed is still little faster..
Rate la di ch jg ch_opt
la 0.327/s  99% 99% 99% 99%
di 25.3/s 7628%  48% 52% 54%
ch 48.6/s 14733% 92%  7% 13%
jg 52.3/s 15877% 107% 8%  6%
ch_opt 55.6/s 16875% 120% 14% 6% 
There are no rules, there are no thumbs..
Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.
 [reply] [d/l] 

johngg, would you please explain what you are doing with vec. I read the doc on it and got lost in the first sentence. I don't know how you are using bit vectors (this is the first time I've ever heard of them). Everything between my $sqrtLimit = sqrt $limit; to return @primes; is a big mystery to me.
You are also not eliminating numbers which end with 2, 4, 5, 6, 8, and 0 right off the bat, why?
Thanks for stopping by.
No matter how hysterical I get, my problems are not time sensitive. So, relax, have a cookie, and a very nice day!
Lady Aleena
 [reply] [d/l] [select] 

Re^2: Number functions I have lying around
by Discipulus (Monsignor) on Mar 31, 2015 at 11:04 UTC

And you can profit of an enhencemt if you too add if($i%2==0){next} before eleborating the square root, as you can see in the ch_opt row.
Rate la di ch ch_opt
la 0.329/s  99% 99% 99%
di 25.3/s 7600%  50% 56%
ch 50.4/s 15230% 99%  12%
ch_opt 57.6/s 17417% 127% 14% 
L*
There are no rules, there are no thumbs..
Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.
 [reply] [d/l] [select] 

sub primes {
my $n = shift;
return if $n < 2;
my @primes = (2);
I: for my $i (3 .. $n) {
next unless 0x208a28aa & (1 << $i % 30);
my $sqrt = int sqrt $i;
for my $p (@primes) {
next I unless $i % $p;
last if $p > $sqrt;
}
push @primes, $i;
}
return @primes
}
 [reply] [d/l] 

Actually, it's 0x208a2882, but you have to initialize primes to (2, 3, 5) then..
 [reply] 
Re^2: Number functions I have lying around
by Lady_Aleena (Curate) on Mar 31, 2015 at 18:28 UTC

A few questions for you choroba...
 Why aren't you eliminating numbers which end with 2, 4, 5, 6, 8, and 0 right off the bat?
 If sqrt($number) == int(sqrt($number), then you wouldn't have to go through the @primes loop, right?
 What is the @primes loop doing?
 I think $n stands for "number", but what do $i and $p represent??
Thank you for stopping by.
No matter how hysterical I get, my problems are not time sensitive. So, relax, have a cookie, and a very nice day!
Lady Aleena
 [reply] [d/l] [select] 

Here's the algorithm in plain words: Let's create the list of primes up to $n. We start with just 2 as the known prime. Then, for each number $i between 3 and $n, we do the following: we try to divide the number $i by all the known primes up to sqrt $i. If any of them divides the number, then it can't be prime. If none of them divides it, it is a prime, though: because a) if a nonprime $d divides $i, then $d = $p1 * $d1, where $p1 is prime, and $p1 divides $i; b) if a number $p2 greater than sqrt $i divides $i, then $i / $p2 must be less than sqrt $i, and it must divide $i. If we find a new prime, we push it to the list.
 I don't eliminate numbers ending with 2, 4, 5, 6, 8, and 0, because they get eliminated in the 0 == $i % $p test.
 testing every number for sqrt $i == int sqrt $i wouldn't help us much, as it happens rarely.
 the @primes loop, as described above, tries to divide the candidate $i by all the known primes up to sqrt $i, to check its primality.
 $n represents the highest number, we are interested in primes less or equal $n. $i is the candidate, i.e. the number we might include in the @primes list if it passes the primality test. $p is a known prime less or equal sqrt $i.
 [reply] [d/l] [select] 

++ to everyone who contributed to this thread, choroba, johngg, atcroft among others for their examples and efforts to explain.
As already said I found the choroba's code faster than the johngg's marvellous one, given the first one was optimized with a premature return in the foreach loop adding if($i%2==0){next}
.
It seems reasonably; there are a lot of numbers (in fact an half of the whole..) that divides by 2. Also a lot that divides by 3. In fact adding this premature return check for number 3 was even faster.
So my (erroneous) thought was: maybe a check against all yet obtained primes can be the best optimization. I added a succint map { next if $i % $_ == 0 } @primes; at the top of the for loop and it revealed to be slow like a dead snail:
Rate primes_all ch jg ch_opt2
primes_all 1.75/s  96% 97% 97%
ch 49.8/s 2751%  4% 11%
jg 51.9/s 2869% 4%  8%
ch_opt2 56.2/s 3117% 13% 8% 
So the fastest code need to check for a few primes not all. But how many? Adding a check for 5, 7 .. seemed to made the code faster. But i'm to lazy to cut and paste a lot of code so i wrote a code generator to have all optimized subs for primes from 2 to 149.
UPDATE:The section below is based on a wrong assumption! as spotted by choroba
And i had a simpatic 1446 lines program ready to overheat my CPUs. I run it few times and results vary but the winner seems to be a prime number between 61 and 89.
In fact preformance increse constantly for primes between 2 and 53 and decrease constantly for primes from 101 to 149.
So choosen 79 as the winner, the fastest code to check primality for numbers from 1 to 10000 can be as ugly as:
sub primes_opt_till_79 {
my $n = shift;
return if $n < 2;
my @primes = (2);
for my $i (3 .. $n) {
if($i%2==0){next}
if($i%3==0){next}
if($i%5==0){next}
if($i%7==0){next}
if($i%11==0){next}
if($i%13==0){next}
if($i%17==0){next}
if($i%19==0){next}
if($i%23==0){next}
if($i%29==0){next}
if($i%31==0){next}
if($i%37==0){next}
if($i%41==0){next}
if($i%43==0){next}
if($i%47==0){next}
if($i%53==0){next}
if($i%59==0){next}
if($i%61==0){next}
if($i%67==0){next}
if($i%71==0){next}
if($i%73==0){next}
if($i%79==0){next}
my $sqrt = sqrt $i;
my $notprime;
for my $p (@primes) {
last if $p > $sqrt;
$notprime = 1, last if 0 == $i % $p;
}
push @primes, $i unless $notprime;
}
return @primes
}
And for completeness an example of benchmark results (cutted to fit here)
Rate
primes_original 74.2/s
opt_till_2 91.8/s
opt_till_3 108/s
opt_till_5 119/s
opt_till_7 128/s
opt_till_11 134/s
opt_till_13 140/s
opt_till_17 142/s
opt_till_19 151/s
opt_till_23 156/s
opt_till_29 158/s
opt_till_31 164/s
opt_till_37 165/s
opt_till_149 169/s
opt_till_139 169/s
opt_till_41 172/s
opt_till_131 174/s
opt_till_137 174/s
opt_till_127 175/s
opt_till_47 176/s
opt_till_113 178/s
opt_till_43 180/s
opt_till_109 181/s
opt_till_53 181/s
opt_till_107 184/s
opt_till_59 184/s
opt_till_103 186/s
opt_till_101 188/s
opt_till_61 188/s
opt_till_97 189/s
opt_till_67 190/s
opt_till_71 190/s
opt_till_73 191/s
opt_till_89 191/s
opt_till_83 192/s
opt_till_79 193/s
L*
There are no rules, there are no thumbs..
Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.
 [reply] [d/l] [select] 








