There's more than one way to do things PerlMonks

### Re: Number functions I have lying around

by choroba (Archbishop)
 on Mar 31, 2015 at 08:55 UTC ( #1121951=note: print w/replies, xml ) Need Help??

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 {
}

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.06e-02 12662%     --
لսႽ† ᥲᥒ⚪⟊Ⴙᘓᖇ Ꮅᘓᖇ⎱ Ⴙᥲ𝇋ƙᘓᖇ

Replies are listed 'Best First'.
Re^2: Number functions I have lying around
by Discipulus (Abbot) 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 Tk-Tartaglia. 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.

Eratosthenes was a clever chap!

```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.

Cheers,

JohnGG

"Eratosthenes was a clever chap"

I never met him. But another chap wrote some kind of shortcut:

```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»

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.

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!
Re^2: Number functions I have lying around
by Discipulus (Abbot) 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.

A modest improvement in speed may be obtained if one implements the wheel.

For example (the modulus is relatively cheap in perl):

```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
}

Actually, it's 0x208a2882, but you have to initialize primes to (2, 3, 5) then..

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!
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 non-prime \$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.
لսႽ† ᥲᥒ⚪⟊Ⴙᘓᖇ Ꮅᘓᖇ⎱ Ⴙᥲ𝇋ƙᘓᖇ
++ 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.

Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://1121951]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (10)
As of 2022-01-18 13:33 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
In 2022, my preferred method to securely store passwords is:

Results (53 votes). Check out past polls.

Notices?