Matthew van Eerde's web log

  • Matthew van Eerde's web log

    Teaching someone to fish and the AKS primality test

    • 0 Comments

    This morning my wife (whom I love and adore) woke me up at 3:00 AM with an urgent question.

    "Hey!" she said, shaking me awake.

    "Is 19 prime?"

    ...

    Like a fool, I answered her.  "Yes.  Yes it is."  Off she went.

    This is a true response, but not a correct response.  I realized shortly afterwards that a correct response would look more like:

    I'm glad you asked me that. dear.  Eratosthenes, the Greek mathematician, discovered a very efficient way to list primes in about 200 BC that is still in use today.  You start by writing out all the numbers from 1 to 19: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19.  1 is a very special number (it's the multiplicative identity, or what algebraists would call a unit) so we put a square around it.  The first number we didn't consider was 2, so we circle it - that means it's prime - and then cross out every multiple of 2 after that.  Going back, the first number we didn't consider was 3... and so on until we get 1 2 3 X 5 X 7 X X X 11 X 13 X X X 17 X 19.  A common optimization is to realize that after circling a prime p, the first number you cross out (that wasn't crossed out before) is always p2, which means that after circling a number you can immediately jump to its square, and also means you can stop crossing out altogether once you hit p> N...

    This would allow her to fish rather than waking me up when she wanted a fish.

    An even better response would have been:

    It's funny you've asked me that.  Number theorists and cryptanalysts have considered this question for thousands of years.  Eratosthenes' method (see above) is a very simple way to find all the primes below a given number, but an efficient way to determine whether a given number is prime was found only very recently.

    In practice, the test that is usually used is the randomized version of the Miller-Rabin test.  Although this is nondeterministic, it is very fast indeed, and will tell you to a very high degree of certainty whether the given number is prime.  This usually suffices.

    There is a deterministic version of the Miller-Rabin test too, which is guaranteed to tell you with perfect certainty whether the given number is prime.  But it only works if you believe in the generalized Riemann hypothesis.  Most mathematicians nowadays believe the hypothesis, but no-one has (yet) been able to prove it.

    Amazingly, in 2002 three mathematicians named Manindra Agrawal, Neeraj Kayal, and Nitin Saxena came up with a deterministic, proven, polynomial-time (specifically, polynomial in the number of digits in the input) method for telling whether a given number is prime.   This is known as the AKS primality test.  The most striking thing about this test is its simplicity - if something this straightforward can be found after thousands of years of looking, what else out there remains to be found?

    Such a response would probably prevent her from waking me again for any mathematical problem at all.  Boo-ya.

    Here's Agrawal, Kayal, and Saxena's "PRIMES is in P" paper.

    Here's Yves Gallot's C++ implementation of AKS.

    My own Perl implementation follows:

    use strict;
    # use bignum;

    print <<USAGE and exit 0 unless @ARGV;
    $0 [-v] n
        Use the AKS primality test to check whether n is prime
        -v adds verbose log spew
    USAGE

    sub is_power($$);
    sub ceil_log2($);
    sub first_r($$);
    sub check_gcds($$);
    sub check_polynomials($$$);
    sub gcd($$);
    sub totient($);
    sub polypow($$\@);
    sub polymult($$\@\@);
    sub polyeq(\@\@);

    my $verbose = $ARGV[0] eq "-v";
    shift @ARGV if $verbose;

    die "Expected only one argument" unless 1 == @ARGV;
    my $n = shift;

    # step 0: restrict to integers >= 2
    print "$n is not an integer and so is NEITHER PRIME NOR COMPOSITE\n" and exit 0 unless int($n) == $n;
    print "$n < 2 and so is NEITHER PRIME OR COMPOSITE\n" and exit 0 unless $n >= 2;

    # step 1: check if the number is a power of some lower number.
    # this can be done quickly by iterating over the exponent (2, 3, ...)
    # and doing a binary search on the base.
    # we start at the top and work down for performance reasons;
    # several subroutines need to know ceil(log2(n)) so we calculate it once and pass it around.
    my $log2_n = ceil_log2($n);
    is_power($n, $log2_n) and exit 0;
    print "Not a power.\n";

    # step 2: find the smallest r such that o_r(n) > (log2 n)^2
    # where o_r(n) is the multiplicative order of n mod r
    # that is, the smallest k such that n^k == 1 mod r
    my $r = first_r($n, $log2_n);
    print "r = $r\n";

    # step 3: for all a between 2 and r inclusive, check whether gcd(a, n) > 1
    check_gcds($n, $r) or exit 0;

    # step 4: if r >= n, we're done
    if ($r >= $n) {
        print "$r >= $n so $n is PRIME\n";
        exit 0;
    }

    # step 5: for all a between 1 and floor( sqrt(phi(r)) log2(n) )
    # check whether (x + a)^n = x^n + a mod x^r - 1, n
    check_polynomials($n, $r, $log2_n) or exit 0;

    # step 6: if we got this far, n is prime
    print "$n is PRIME\n";

    sub is_power($$) {
        my $n = shift;
        my $log2_n = shift; # actually ceil(log2(n))

        print "Checking for power-ness...\n";

        # we consider numbers of the form b^i
        # we iterate over the exponent i
        # starting at i = ceil(log2(n)) and working down to i = 2
        #
        # for each exponent we do a binary search on the base
        # the lowest the base can be is 2
        # and the highest the base can be (initially) is 2
        #
        # we set up bounds on the base that are guaranteed to
        # surround the actual base
        my $b_low = 1; # 1 ^ ceil(log2(n)) = 1 < n
        my $b_high = 3; # 3 ^ ceil(log2(n)) > 2 ^ log2(n) = n

        for (my $i = $log2_n; $i >= 2; $i--) {
            print "\tb^$i\n" if $verbose;

            # let's check that the bounds are really correct
            die "$b_low ^ $i is not < $n" unless $b_low ** $i < $n;
            die "$b_high ^ $i is not > $n" unless $b_high ** $i > $n;

            # do a binary search to find b such that b ^ i = n
            while ($b_high - $b_low > 1) {
                print "\t\tb^$i: b is between $b_low and $b_high\n" if $verbose;
                my $b = int(($b_low + $b_high)/2);

                my $t = $b ** $i;
                if ($t == $n) {
                    print "$n = $b^$i; $n is COMPOSITE\n";
                    return 1;
                }

                ($t > $n ? $b_high : $b_low) = $b;
            }

            # as we pass from the exponent (say, 5)
            # to the exponent below (say, 4)
            # we need to reconsider our bounds
            #
            # b_low can remain the same because b ^ (i - 1) is even less than b ^ i
            # OPEN ISSUE: can we even raise b_low?
            #
            # but we need to raise b_high since b ^ i > n does NOT imply b ^ (i - 1) > n
            #
            # we'll square b_high; b ^ i > n => (b ^ 2) ^ (i - 1) = b ^ (2 i - 2) > n
            # since i >= 2
            #
            # OPEN ISSUE: is there a better way to raise this higher bound? Does this help much?
            $b_high *= $b_high;
        }

        # nope, not a power
        return 0;
    }

    sub ceil_log2($) {
        my $n = shift;

        my $i = 0;
        my $t = 1;

        until ($t >= $n) {
            $i++;
            $t *= 2;
        }

        return $i;
    }

    sub first_r($$) {
        my $n = shift;
        my $log2_n = shift; # actually ceil(log2(n))

        my $s = $log2_n ** 2;
        print "Looking for the first r where o_r($n) > $s...\n";

        # for each r we want to find the smallest k such that
        # n^k == 1 mod r

        my $r;
        for ($r = 2; ; $r++) {
            # print "\tTrying $r...\n";

            # find the multiplicative order of n mod r
            my $k = 1;
            my $t = $n % $r;

            until (1 == $t or $k > $s) {
                $t = ($t * $n) % $r;
                $k++;
            }

            if ($k > $s) {
                # print "\to_$r($n) is at least $k\n";
                last;
            } else {
                # print "\to_$r($n) = $k\n";
            }
        }

        return $r;
    }

    sub check_gcds($$) {
        my ($n, $r) = @_;

        print "Checking GCD($n, a) for a = 2 to $r...\n";

        for (my $a = 2; $a <= $r; $a++) {
            my $g = gcd($n, $a);

            next if ($g == $n); # this is OK

            if (1 != $g) {
                print "gcd($n, $a) = $g; $n is COMPOSITE\n";
                return 0;
            }
        }

        print "All GCDs are 1 or $n\n";

        return 1;
    }

    sub gcd($$) {
        my ($x, $y) = @_;

        ($x, $y) = ($y, $x) unless $x > $y;

        while ($y) {
            ($x, $y) = ($y, $x % $y);
        }

        return $x;
    }

    sub check_polynomials($$$) {
        my $n = shift;
        my $r = shift;
        my $log2_n = shift; # actually ceil(log2(n))

        # iterate over a from 1 to floor( sqrt(phi(r)) log2(n) )
        # for each a, check whether the polynomial equality holds:
        # (x + a)^n = x^n + a mod (x^r - 1, n)
        # if it fails to hold, the number is composite
        #
        # first we need to evaluate phi(r) so we can determine the upper bound
        # OPEN ISSUE: this seems to be a potential weakness in the algorithm
        # because the usual way to evaluate phi(r) is to find the prime factorization of r
        # and then form the product r*PI(1 - 1/p) where the product ranges over all primes
        # which divide r

        my $phi = totient($r);

        # a < sqrt(phi(r)) * log2(n) => a^2 < phi(r) * (log2(n))^2
        my $a2_max = $phi * $log2_n * $log2_n;
        print "Checking polynomials up to roughly ", int sqrt($a2_max), "...\n";

        for (my $a = 1; $a * $a <= $a2_max; $a++) {
            print "\ta = $a...\n" if $verbose;

            # polynomials are of the form (c0, c1, c2, ..., ci, ...)
            # which corresponds to c0 + c1 x + c2 x^2 + ... + ci x^i + ...)
            my @x = (0, 1);
            my @x_plus_a = ($a % $n, 1);

            my @lhs = polypow($n, $r, @x_plus_a);

            # POTENTIAL OPTIMIZATION:
            # x^n + a mod (x^r - 1) is just x^(n % r) + a
            # and we know n % r != 0
            my @rhs = polypow($n, $r, @x); # x^n
            $rhs[0] = ($rhs[0] + $a) % $n; # + a

            next if polyeq(@lhs, @rhs);

            print "(x + $a)^$n is not equal to x^$n + $a mod(x^$r - 1, $n)\n";
            print "So $n is COMPOSITE\n";
            return 0;
        }

        return 1;
    }

    sub totient($) {
        my $r = shift;

        print "Finding the Euler totient of $r\n";

        # we'll do a trial division to find the totient
        # there are faster ways that use a sieve
        # but we don't know how big r is
        my $t = $r;

        # by construction p will always be prime when it is used
        # OPEN ISSUE: this might be slow
        for (my $p = 2; $r > 1; $p++) {
            next if $r % $p;

            print "\t$p is a factor\n" if $verbose;
            # decrease the totient
            $t /= $p;
            $t *= $p - 1;

            # decrease r
            $r /= $p; # we know there's at least one factor of p
            $r /= $p until $r % $p; # there might be more
        }

        print "Totient is $t\n";

        return $t;
    }

    sub polypow($$\@) {
        my $n = shift; # this is both the mod and the exponent
        my $r = shift;
        my @base = @{ +shift };

        my $exp = $n;
        my @result = (1); # 1

        # print "\t(", join(" ", @base), ")^$exp mod (x^$r - 1, $n)\n" if $verbose;

        # basic modpow routine, but with polynomials
        while ($exp) {
            if ($exp % 2) {
                @result = polymult($n, $r, @result, @base);
            }

            $exp = int ($exp / 2);
            @base = polymult($n, $r, @base, @base);
        }

        # print "\t= (", join(" ", @result), ")\n" if $verbose;
        return @result;
    }

    sub polymult($$\@\@) {
        my $n = shift;
        my $r = shift;
        my @first = @{ +shift };
        my @second = @{ +shift };

        # print "\t\t(", join(" ", @first), ") * (", join(" ", @second), ") mod (x^$r - 1, $n)\n" if $verbose;

        my @result = ();

        # first do a straight multiplication first * second
        my $s = @second - 1;
        for (my $i = @first - 1; $i >= 0; $i--) {
            for (my $j = $s; $j >= 0; $j--) {
                my $k = $i + $j;
                $result[$k] += $first[$i] * $second[$j];
                $result[$k] %= $n;
            }
        }

        # then do a straight mod x^r - 1
        # consider a polynomial
        # c0 + ... + ck x^k
        # with k >= r
        # we can subtract ck (x^r - 1)
        # without changing the mod value
        # the net effect is to eliminate the x^k term
        # and add ck to the x^(k - r) term

        for (my $i = @result - 1; $i >= $r; $i--) {
            my $j = $i - $r;
            $result[$j] += $result[$i];
            $result[$j] %= $n;

            pop @result;
        }

        # eliminate any leading zero terms
        for (my $i = @result - 1; 0 == $result[$i]; $i--) {
            pop @result;
        }

        # print "\t\t= (", join(" ", @result), ")\n" if $verbose;
        return @result;
    }

    sub polyeq(\@\@) {
        my @lhs = @{ +shift };
        my @rhs = @{ +shift };

        # print "(", join(" ", @lhs), ") = (", join(" ", @rhs), ")?\n" if $verbose;

        return 0 unless @lhs == @rhs;

        for (my $i = @lhs - 1; $i >= 0; $i--) {
            return 0 unless $lhs[$i] == $rhs[$i];
        }

        return 1;
    }

    Here's the output when I run it on 19:

    >perl -w aks.pl 19
    Checking for power-ness...
    Not a power.
    Looking for the first r where o_r(19) > 25...
    r = 19
    Checking GCD(19, a) for a = 2 to 19...
    All GCDs are 1 or 19
    19 >= 19 so 19 is PRIME

    And here's the output with a bigger input:

    >perl -w aks.pl 99
    Checking for power-ness...
    Not a power.
    Looking for the first r where o_r(997) > 100...
    r = 103
    Checking GCD(997, a) for a = 2 to 103...
    All GCDs are 1 or 997
    Finding the Euler totient of 103
    Totient is 102
    Checking polynomials up to roughly 100...
    997 is PRIME
  • Matthew van Eerde's web log

    Programmatically setting a local user account to never expire its password

    • 1 Comments

    As a Windows tester, I install Windows on my own machines a lot (this is known internally as "selfhosting", or "dogfooding", or "ice cream-ing".)

    One of my little idiosyncracies is I like to run as a non-administrative user.  That is, I don't add my domain account to the local Administrators group.

    Instead, I create a local "Admin" account with a known (to me) password; every time I need to elevate, I get a prompt that asks for credentials rather than just "Yes/No".  To this prompt I pass the credentials of the local "Admin" account.

    Although I usually install fresh builds regularly (on my multiple machines), sometimes one machine gets a little stale.  In fact, it happened once that my local .\Admin account got so stale that I had to change the password!  This was annoying enough that I devoted some energy into figuring out how to check the "Password never expires" box on the local account properties programmatically.

    The result was the following script: call as cscript.exe never-expire-admin-password.wsf  This version hardcodes the username "Admin"; a production version would probably allow passing a username in via the command line.

    If the Admin password already has the box checked, this script does nothing.

    ' LDAP doesn't work for controlling local users
    ' (unless you're a domain controller, of course)
    '
    ' have to use WinNT provider instead

    Const ADS_UF_DONT_EXPIRE_PASSWD = &H10000

    ' hardcoding "Admin" username
    Dim admin: Set admin = GetObject("WinNT://localhost/Admin,user")

    WScript.Echo "Admin's userFlags are 0x" & Hex(admin.userFlags)

    If Not admin.userFlags And ADS_UF_DONT_EXPIRE_PASSWD Then
        WScript.Echo "Setting local admin account to never expire password"
        admin.userFlags = (admin.userFlags Or ADS_UF_DONT_EXPIRE_PASSWD)

        ' Save
        admin.SetInfo
    End If

Page 1 of 1 (2 items)

February, 2012