(chibi math prime) Miller-Rabin primality test is not correctly implemented

[wrog]

(in /lib/chibi/math/prime.scm)
modular-root-of-one? is not doing the right thing and this will result in further false positives (i.e., numbers deemed prime that actually aren't) beyond what the Miller-Rabin test is normally expected to produce.

First of all, the docstring:

;; Returns true iff any (modular-expt a odd*2^i n) for i=0..twos-1
;; returns 1 modulo n.

does not match up with what the routine is actually doing, which is that it returns true if any of the a^(odd*2^i) powers encountered are 1 or -1.

But what we actually need this to do is return true if either a^odd = 1 (and thence all squares of it are 1 and we're getting no information on square roots of 1) or that we encounter -1 first (and again all squares will be 1). If we encounter 1 first and it's the square of something that's not -1 or 1, then 1 has a third square root, we're not in a field, and so we have a witness, a, to n not being prime.

An example where modular-root-of-one? is going wrong would be
n=15 (neg1=14, twos=1, odd=7) and a=4,
which is an actual witness to 15 not being prime because a^7=4 and a^14=1 (mod 15) and thus we have 4 being a third square root of 1 (along with 1 and 14)

Since this routine isn't exported and is only used once, it would probably best for the sake of readability to reverse the sense of it, rename it to miller-rabin-composite-witness?, i.e., so that it's returning true when it finds a witness. Meaning we replace

(not (modular-root-of-one? ...))

with

(miller-rabin-composite-witness? ...)

in the body of miller-rabin-composite? And then

(define (miller-rabin-composite-witness? twos odd a n neg1)
  ;; Returns true iff (1) a^(n-1)!=1 mod n OR
  ;; (2) the powers of a include a square root of 1 other than {1, -1}.
  ;; (which then means Z/(n) cannot be a field and n cannot be prime)
  (let ((b (modular-expt a odd n)))
    (let lp ((i 0) (b b))
      (cond ((= b neg1)  #f)               ; -1 is an expected sqrt(1)
            ((= b 1)     (not (zero? i)))  ; !! previous b is bad sqrt(1)
            ((>= i twos))                  ; !! a^(n-1)!=1 mod n 
            (else (lp (+ i 1) (remainder (* b b) n))))))))

[wrog]

In other news, factor is quite a bit slower than it needs to be. (once you've found a factor you don't need to continue searching all the way up to sqrt(original n), sqrt(new n = n/p^k)` suffices).

(factor 0) should raise an error (every prime divides 0, there's no reasonable result you could return here; nobody's going to be doing this except by mistake.

(factor 1) should be an empty list, not (1), the elements of the returned list should be primes and 1 isn't (also, being able to include 1s in a factor list would break unique factorization). product of the list should be the original number; product of an empty list is 1, so there's no need to put 1s in there.

And there are well-known formulas for 'totient' and 'aliquot' based on the factorizations that give multiple orders of magnitude speedup (searching for factors is way faster than doing a gcd on every number in the range (gcd takes multiple divisions)).

(totient 1) = 1, not 0. Totient is the number of integers k in the range 1 ≤ k ≤ n for which the greatest common divisor gcd(n, k) is equal to 1 (see wikipedia). Things people care about include having this function be multiplicative for numbers that are relatively prime (i.e., φ(ab) = φ(a)φ(b) if gcd(a,b)=1), and also having the sum of φ(d) over all d that divide n (including d=1 and d=n) is supposed to give you n and these only work if φ(1) = 1.

(totient 0) has no agreed definition and should probably raise an error (and using the Euler product formula based on the factorization will make this automatic).

(aliquot 1) = 0, not 1. It's a sum of proper divisors, which are the ones that are strictly less, so 1 doesn't have any. (see wikipedia).

(aliquot 0) is arguably nonsensical (if it's a sum of divisors, then it's infinite), and so should almost certainly raise an error (and using the formula based on the factorization will make this automatic).

I think I'll be doing a pull request.

[ashinn]

Hi Roger, thanks very much for the detailed comments and sorry for the slow reply. I had actually been meaning to make the factor optimization you suggest, though had not really thought about optimizing totient or aliquot.

I still need to go through the miller rabin code and come up with a suitable test case.

[wrog]

I still need to go through the miller rabin code and come up with a suitable test case.

That may be hard. Unless there's a feature of the test framework to allow sabotaging the innards of the module with fluid-let (specifically to make the random number generator really unlucky and reduce the prime-table list). That's about the only route I see to getting something deterministic/reproducible.

[ashinn]

Non-portable but explicitly setting the state of default-random-source should be deterministic on a given platform.

But as you say changing the innards is an option - we can make all of the random-integer calls refer to a random-source from an exported parameter instead of default-random-source, and parameterize that in the tests.

[wrog]

I think I'll be doing a pull request.

and done.

[wrog]

(While I have this swapped in, I thought I'd look at the rest of the module.)
So there's this issue with random-coprime:

> (let* ((p 12)   (hist (make-vector (+ p 1) 0)))
>   (let loop ((n 10000))
>     (let ((i (random-coprime p)))
>       (vector-set! hist i (+ (vector-ref hist i) 1)))
>       (unless (zero? n) (loop (- n 1))))
>    hist))

#(0 0 0 0 0 3617 0 1829 0 0 0 4555 0)

Two aspects of this that are not ideal:

• The possible outcomes (in this case, 1, 5, 7, 11) are not equally probable
• We're not seeing 1 at all, even though 1 is coprime to everything.

I'll grant that uniformity isn't actually promised but I'm having trouble seeing how a routine like will actually be useful without either a guarantee of uniformity or some kind of statement concerning what the distribution of outcomes is going to be like.

And, as it happens, there's a fairly easy way to make this uniform:

(define (random-coprime n)
  (let ((n-1 (- n 1)))
    (let loop ()
      (let ((p (+ 1 (random-integer n-1))))
        (if (coprime? p n) p (loop))))))

it's slower but not that much slower than the current implementation (worst case for numbers < 2^60 is 614889782588491410, (square-free with lots of distinct prime factors, where we seem to be about 3/4 as fast and the expected number of RNG calls is around 8)., and it seems like there ought to be something we can do with the Chinese Remainder Theorem to speed this up..

I'm also pretty sure a similar issue exists with random-prime.

[ashinn]

The distribution will be skewed by the natural distribution of primes and coprimes - (co)primes coming after longer gaps are more likely.

Your approach is probably better, even if slower, but please make it a separate pull request.

random-prime is used by the RSA implementation for very large primes, where the slightly skewed distribution is unlikely to be an effective attack vector.