All the Perl that's Practical to Extract and Report

is_prime(number) in Perl

Ovid on 2008-02-23T14:54:05

As a general rule, if you want to check if a number is prime, it turns out that the deterministic solutions are rather slow. In fact, they tend to be too slow for general use. So I figured, while writing Number::Properties that is_prime() should use the Miller-Rabin primality tests. My code had some issues, but I suspect this had something to do with Perl's lack of precision. As a result, and since Perl is a bit slow for work like this, I figured I would write the code in XS. I used the C from Top Coder's primality testing article, but quickly got stuck with the following error:

Error: 'long long' not in typemap in Properties.xs, line 64
Please specify prototyping behavior for Properties.xs (see perlxs manual)

Much googling suggests to me that there's not a simple way around this that's also robust. Damn. I may have to switch back to Perl.

That implies...

AndyArmstrong on 2008-02-23T15:09:30

...that you're trying to move a long long between C and Perl. Do you need to do that? Can't the long longs just live inside the C?

It's probably not too hard to write the code to thunk between long long and Math::BigInt - would that help?

Re:That implies...

Ovid on 2008-02-23T16:02:55

I don't know since my C and XS skills are so poor. The relevant bit of XS code is this:

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include <stdbool.h>

#include "ppport.h"

/* This function calculates (ab)%c */
int modulo(int a,int b,int c){
    long long x=1,y=a; // long long is taken to avoid overflow
                       // of intermediate results
    // skipping implementation
}

/* this function calculates (a*b)%c taking into account
   that a*b might overflow */
long long mulmod(long long a,long long b,long long c){
    long long x = 0,y=a%c;
    // skipping implementation
}

/* Miller-Rabin primality test, iteration signifies the accuracy of the test */
bool Miller(long long p,int iteration){
    if(p<2){
        return false;
    }
    if(p!=2 && p%2==0){
        return false;
    }
    long long s=p-1;
    while(s%2==0){
        s/=2;
    }
    int i;
    for(i=0;i<iteration;i++){
        long long a=rand()%(p-1)+1,temp=s;
        long long mod=modulo(a,temp,p);
        while(temp!=p-1 && mod!=1 && mod!=p-1){
            mod=mulmod(mod,mod,p);
            temp *= 2;
        }
        if(mod!=p-1 && temp%2==0){
            return false;
        }
    }
    return true;
}

MODULE = Number::Properties     PACKAGE = Number::Properties

void Miller(p, iteration)
    long long p
    int iteration

In other words, I just need to figure out the interface to the Miller() function and this should work, but I'm unsure.

Re:That implies...

AndyArmstrong on 2008-02-23T16:19:08

If you're prepared to limit the input to 32 bit integers you can probably just change the long long to long in the prototype of Miller(). It should still use long long internally.

To be more strictly correct you could use the output of

perl -MConfig -le 'print $Config{ivtype}'

as the type of that argument. That will be long for 32 bit Perl and long long for 64 bit Perl.

And if you want to be really sneaky write a version that bypasses XS type mapping and correctly handles Math::BigInt objects :)

Re:That implies...

bart on 2008-02-23T18:10:16

If you're prepared to limit the input to 32 bit integers
A plain float in Perl has 53 bits of precision, which means you can represent 53 bit integers exactly, which is better than merely 32 bits.

Or, am I missing a point?

Re:That implies...

AndyArmstrong on 2008-02-23T18:21:06

No, I think you're quite right :)

This little program:

#!/usr/bin/perl

my $num  = 1;
my $bits = 0;
while ( $num + 1 != $num ) {
    $num *= 2;
    $bits++;
}

print "$bits\n";

prints 53.

So Ovid - make that argument a double and then cast it to a long long inside the function.

Re:That implies...

renodino on 2008-02-23T18:19:48

If (and its a big if) you need the full 64 bit input, might I suggest passing the number as a string/PV, then doing a strtoll() or its equivalent in your non-XS functions ?

Then your XS looks like:

void Miller(char *p, int iteration)

and the Miller function becomes

bool Miller(char *p_as_str,int iteration){
    long long p = stroll(p_as_str);
    if(p<2){

Caveat: obviously, if you put something other than a number in p_as_str, things will go awry. Also, different platforms (e.g., Windows) have variants of the stroll().

Re:That implies...

renodino on 2008-02-23T21:52:33

DOH! That should be strtoll()...

The Updated AKS primality test is slow?

Limbic Region on 2008-02-23T17:04:49

I assume you were referring to

http://en.wikipedia.org/wiki/AKS_primality_test

The key significance of AKS is that it was the first published primality-proving algorithm to be simultaneously general, polynomial, deterministic, and unconditional. Previous algorithms have achieved any three of these properties, but not all four.

I haven't seen an implementation (in perl or otherwise), but I thought it put to bed the problem of determining if a number was prime or not. The real problem is factorizing a composite.

Fast primality testing with Pari

pjf on 2008-02-23T23:00:46

If you're after a fast way of checking for prime numbers from Perl, then it's hard to go past PARI/GP, which has a Math::Pari interface. With the generation of large primes, strong random numbers, support for integers up to (2^29)^29 out of the box, and a host of other math functions, Pari is essentially crack for number theorists. It provides both isprime() and ispseudoprime() functions.

There's also a brief mention of using Math::Pari for primality testing on perlmonks.

Cheerio,

Paul