 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 { # 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.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?