Beefy Boxes and Bandwidth Generously Provided by pair Networks
Just another Perl shrine
 
PerlMonks  

Re: Number functions I have lying around

by choroba (Bishop)
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 (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 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!
        Lady Aleena
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.

      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!
    Lady Aleena
      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.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1121951]
help
Chatterbox?
[Corion]: Discipulus: Ouch ... staff cuts are never good, but having to hunt for work isn't great either, especially with a family...

How do I use this? | Other CB clients
Other Users?
Others wandering the Monastery: (10)
As of 2018-05-23 21:18 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    Notices?