Well, about the onlything I can add is that I looked into Knuth's AoP and found algoithm P (1.3.2) which is for calculating the first N primes. I found it interesting as it was bit hard to convert the pseudocode* he had to perl. I switched around your code to do the same, and here they are:

sub make_primes_ambrus {
    my $count=shift;
    my @primes=(2);
    my( $n, $p, $sqrp);
    NLOOP:
    for ($n=3; @primes<$count; $n+=2) {
        $sqrp = sqrt($n) + 0.5;
        for $p (@primes) {
                $sqrp < $p and last;
                0 == $n % $p and next NLOOP;
        }
        push @primes, $n
    }
    return \@primes;
}

sub make_primes_knuth  {
    my $count=shift;
    my @primes=(2);
    my $n=3;
    PRIME:
    while (@primes<$count) {
        push @primes,$n;
        while ( $n+=2 ) {
            for my $p (@primes) {
                my $r=$n % $p;
                last if !$r;
                # According to Knuth the proof of the vaildity 
                # of the following test "is interesting and a little 
                # unusual".
                next PRIME if int($n/$p) <= $p;
            }
        }
    }
    return \@primes;
}
[download]

Incidentally, part of the issue for my original problem is that im generating a random prime relatively rarely, and as a fresh process every time. So the time (and hassle coding) to produce a table of primes to use for the div technique is not really necessary. Doing a brute force check on each random number in turn is just fine for me here.

* (Incidentally is there any easy way to get $Q and $R from $x/$y where Q is the quotient and R is the remainder? I used $Q=int($x/$y); and $R=$X % $Y; but i wonder if there is a better way, obviously repeated subtraction would leave an integer remainer somewhere, but given todays chip design is that actually more efficient than doing div and mod as seperate tasks? )

I'm replying to your question about division an modulo together.

In perl, there's no reason to do this, as the opcodes to do the separating of the two results would take much more time than what you'd gain from calculating the two together. However, this might not be the case with big integers, and indeed, the Math::BigInt module defined the bdiv method in such a way that it can store the remainder at once. (Btw the iquo function of Maple works the same way: it stores the other result. This is not surprising as Maple handles big integers.)

In C, you can hope that the compiler notices that you use division and modulus on the same number and will optimize them to use just one operation if that helps. If you don't trust the compiler, you can also use the div (or ldiv, lldiv, imaxdiv) functions of the C library, which return both the quotient and the remainder. I think this doesn't have much merit anyway, as if the compiler cannot do the optimization of a division and a modulus together, then either it wouldn't give any speed gain in the particular case (eg because there's no free register to store the other result) or your compiler is crap and you'd gain more speed by installing a new compiler or fiddling with the compiler flags than doing micro-optimizations on division.

My favourite way to test primality is:

sub is_prime {
    ('1' x shift) !~ /^(11+)\1+$/
}
[download]

:-)

See also prime factorization using base 1.

Maybe I'm missing something, but if you want a "small" (whatever that means) random prime number then generating a random number, testing whether it's a prime and then generating tens and hundreds more until you accidentaly tramp over a prime sounds horribly inefficient to me. And actually not gauranteed to finish at all. Unless your "small" means something like "less than 20" primes are fairly sparse so on average you have to generate quite a few random numbers to find a single prime. For a reasonable definition of "small" you'd be much better off if you just pregenerate a list of all "small" primes and then keep selecting a random item from the list.

primes are fairly sparse so on average you have to generate quite a few random numbers to find a single prime

The density of primes is high enough that this really isn't a big problem. The prime number theorem says that the density of primes around N is roughtly 1/ln(N). So if you want to find a k-digit prime, you should expect to find one after k tries on average. (You can of course be clever and only consider odd numbers, and then you'll find one after k/2 tries on average). So Rabin-Miller takes O(k) time on k-bit numbers, and actually finding a prime takes O(k) tries on average. That's quite a good tradeoff actually.

I'm pretty certain this is how most modern crypto software generates RSA primes (guess and then check with Rabin-Miller). No, it's not super fast, but it's quite reasonable, and you usually only have to do it once.

blokhead

From the Prime number theorem it will take k tries on average for a k digit number only if you are working in base e and looking at numbers right around e**k. For mortals who work in base 10, it takes an average of about 2.3*k+1 tries or so.

Things like being clever about odds works from this starting point.

Firstly, let me explain why random generating a prime isn't that bad.

The first reason is the same as blokhead has mentioned, that primes aren't so rare. As he has explained, if you pick a random 31-bit number, it's a prime with 1/21 probability.

The second reason is that even when we don't pick a prime, in most of the cases, it is divisable with a small prime, and in this case, the primality tester algorithm should return false very quickly. Note that even before the R-M test, I divide the number with the first eight primes (this should really be the first 50 primes in real-world applications, but then I wouldn't get much of the R-M algorithm here) first thing I do. Even without those divisions, the R-M test is correct, so they wouldn't strictly be needed, but it really helps finding a random prime number. As a special case, blokhead notes that you can generate a random even number. I do this, but it doesn't help much as the primality tests test the parity of the number anyway. (Update: it was stupid to say that I do this as I haven't posted the code to generate a random prime. Anyway, here it is:

Read more... (1212 Bytes)
end update)

Now let's look at whether building a table of large primes is feasable, esp as spiritway has asked this too. There are 5 million 8-digit primes, so if we really needed five-digit primes, so we could build a table. However, I'd only do this if you needed five-digit primes very often. (demerphq's original application was a random hash seed, so I don't think he needs it that often.) Also, if you need 9-digit primes or even primes less then 2^31, the table gets larger and larger.

My original requirements were for primes with 7-8 digits, or more generally, those that can be represented by a long.

testing whether it's a prime and then generating tens and hundreds more until you accidentaly tramp over a prime sounds horribly inefficient to me.

Surprisingly it doesnt take very long at all for my particular task. But you are right it isn't efficient. But efficiency in this particular case is not required at all. If the program takes a few seconds then fine, but as it is it finishes in a heartbeat.

---
$world=~s/war/peace/g

I expected slightly smaller "small primes" :-) And I stand corrected regarding the sparseness of primes. Thanks, you never stop learning :-)

(For a minute I thought whether it wouldn't be better to generate just one random number and then keep subtracting/incrementing until you find a prime, so that the algorithm is garanteed to finish, but that's not a good idea. Primes are not evenly distributed so some would come up more often than others.)

Jenda


XML sucks. Badly. SOAP on the other hand is the most powerfull vacuum pump ever invented.

sub is_prime (Int $n) {
    $n % all(2 .. sqrt($n)+1);
}
[download]

Credit: Autrijus in the Apocalypse Now talk.

I don't have perl6 handy, but wouldn't that fail for 2? In that case, sqrt(1) + 1 is greater then 2 so 2 is included in the range.

You are correct. Nice catch!

The all() method originally comes from the Quantum::Superpositions module.

This one and other primality tests were also discussed in Prime Number Finder.

Gu

You might consider generating a list of primes below a given number, and keeping it on disk. Then randomly access one of those numbers. I don't know exactly how many primes are less than 10**8 or so, but it's not altogether overwhelming.

See the other answer I gave in this thread: Re^2: Simple primality testing.

#!/usr/bin/perl -w

# a not too bad approximation to the number of primes 
# less than N is approximately N/log(N).  Better methods 
# exist, but none are as simple...

my $n1 = 10**8;
my $n2 = 2**32;

$\="\n";

# ~ 5428681 
print "Number of primes less than $n1 is ",int($n1/log($n1));

# ~ 193635250
print "Number of primes less than $n2 is ",int($n2/log($n2));
[download]

Actually since these are signed longs it only goes up to 2^31-1, which is a prime, a Mersenne Prime actually. (For those that missed this bit, a Mersenne prime is any prime expressable in the form 2^N-1.)

I dont know if I'm the only one that finds it interesting that the highest number you can represent with a long is a prime, but I do. :-)

Now if only they would use questions like this on pub-quizzes. :-)

---
$world=~s/war/peace/g

I used Inline::C for a 64-bit integer calculation that couldn't easily be done on a perl compiled with 32-bit integers (the precision of doubles wouldn't be enough). Obviously, Math::BigInt or its fast backend Math::BigInt::GMP would work here. However, salva just mentioned his module Math::Int64 which I think would fit here exactly.

See also Rabin Miller.

my $random_prime = [17, 23]->[int rand 2];
[download]

I thought the canonical example prime number was 91.

No, no, no, that's the "hard to factor" product of two random primes.