Randomized Prime (Perl)

From LiteratePrograms

Jump to: navigation, search

1 Overview
2 Operation
- 2.1 Jacobi Symbol
3 Big rand

[edit]

Overview

This script will determine if a number is composite or that it is prime within some known probability.

The output is the word 'composite' or a sentence describing the odds that n is prime.

[edit]

Operation

[edit]

Jacobi Symbol

$\left[\frac{a}{n}\right] = \prod_{i-1}^{t} \left[ \frac{a}{p_i}\right]^{k_i}$

$\left[\frac{ab}{n}\right] = \left[\frac{a}{n}\right] \left[\frac{b}{n}\right]$

property 1
if ($a % 2 == 0) {
  return jacobi(2, $n) * jacobi($a/2, $n);
}

For $a \equiv b \pmod{n}$ , $\left[\frac{a}{n}\right] = \left[\frac{b}{n}\right]$ .
For odd coprimes a and n, $\left[\frac{a}{n}\right] = (-1)^{\frac{a-1}{2}\frac{n-1}{2}} \left[\frac{n}{a}\right]$ .

property 3
if ($a % 4 == 3 && $n % 4 == 3) {
  -jacobi($n, $a);
} else {
  jacobi($n, $a);
}

$\left[\frac{1}{n}\right] = 1$ .

property 4
if ($a == 1) {
  return 1;
}

$\left[\frac{2}{n}\right] = \begin{cases} -1 & \mbox{for } n \equiv 3 \mbox{ or } 5\pmod{8} \\ 1 & \mbox{for } n \equiv 1 \mbox{ or } 7\pmod{8} \\ \end{cases}$ .

property 5
if ($a == 2) {
  my $x = $n % 8;
  if ($x == 1 || $x == 7) {
    return 1;
  } elsif ($x == 3 || $x == 5) {
    return -1;
  }
}

<<jacobi>>=

sub jacobi
{
  my ($a, $n) = @_;

  die "Bad input" unless defined $a and $n;

  if ($a == 0) {                # identity 1
    return ($n == 1) ? 1 : 0;
  }
  if ($a == 2) {                # identity 2
    my $x = $n % 8;
    if ($x == 1 || $x == 7) {
      return 1;
    } elsif ($x == 3 || $x == 5) {
      return -1;
    }
  }
  if ( $a >= $n ) {             # identity 3
    return jacobi($a % $n, $n);
  }
  die if $a % 2 == 0;
                                # identities 5 and 6
  return ($a % 4 == 3 && $n % 4 == 3) ? -jacobi($n, $a) : jacobi($n, $a);
}

[edit]

Big rand

Choose a uniformly at random from $\mathbb{Z}_n$ .

<<brand>>=
sub brand {
   my $number = shift;
   my $n = Math::BigInt->new($number);
   my $l = Math::BigFloat->new($n + 1)->blog(2)->bceil;
   my $r = 0;

   while ($r == 0 or $r >= $n) {
        my $x = Crypt::OpenSSL::Random::random_pseudo_bytes(10);
        $r = Math::BigInt->new('0b' . unpack("B$l", $x));
    }
    return $r;
}

<<check>>=

sub check
{
    my $number = shift;

    my $n = Math::BigInt->new($number);

    return 'Composite 0' if $n->copy->bmod(2) == 0;

    my $a = brand($number);

    return 'Composite 1' if (1 != $a->copy->bgcd($n));

    my $x = ($n - 1) / 2;

    my $q = $a->copy->bmodpow($x, $n);

    my $j = jacobi($a, $n);

    $j %= $n;

    if ($j == $q) {
        return 'Prime.';
    } else {
        return 'Composite 2';
    }
}

<<main>>=

my $number = shift || die "Need a number.";
my $loops = shift || 10;

die if ($loops < 0);

my $c = 0;

while ($c++ < $loops) {
    my $ret = check($number);
    if ($ret =~ /Composite/) {
        print $ret, "\n";
        exit;
    }
}

print "This is only a 1/2^", $loops, " chance that $number is not prime\n";

check
brand
jacobi

<<rand_primality.pl>>=
#!/usr/bin/perl

use strict;

use Math::BigInt::GMP;

# Precondition: a, n >= 0; n is odd
use strict;
use warnings;
#use diagnostics;

use Crypt::OpenSSL::Random;
use Data::Dumper;
use Math::BigInt;
use Math::BigFloat;

main