2008年11月30日星期日

Probabilistic Primitivity Testing

I have been learning Haskell and as an exercise once in a while I will pick up SICP and try to implement some of the programs in Haskell. 

Section 1.2.6 of SICP contains the definition of a function that probabilistically tests if a given natural number n is a prime based on Fermat's Little Theorem. Below is the Scheme code:

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))        

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else false)))

And here is the equivalent Haskell code that I wrote:

> import System.Random
> import System.IO.Unsafe
>
> fast_prime         :: Integer -> Int -> Bool
> fast_prime n times = and $ map fermat_test (mkRandomList n times)
>     where fermat_test a = expmod a n n == a
>           expmod base exp m 
>               | (exp == 0) = 1
>               | (even exp) = (square (expmod base (exp `div` 2) m)) `mod` m
>               | otherwise  = (base * (expmod base (exp -1) m)) `mod` m
>           square a = a * a
>           
> mkRandomList :: Integer -> Int -> [Integer]
> mkRandomList n len = 
>     take len $ 
>         unsafePerformIO $ 
>             do g <- getStdGen
>                return $ randomRs (1, n - 1) g

The major difference between the Haskell code and the Scheme code is the part dealing with the random number generator. Where in Scheme a simple call to random() can do, in Haskell, I have to resort to monad, which makes the Haskell code less cleaner than its Scheme counterpart (which is usually not the case in my experience).  But the part that makes me feel most uncomfortable is the call to unsafePerformIO whose very existence implies something "incorrect" here. 

But I have my reasons: I do not want mkRandomList itself to be a monad because that would force all the functions along the calling chain to be inside a monad which is way too clumsy from the point of view of the callers of that function (fast_prime in this case). After all, conceptually mkRandomList is a pure function whose very responsibility is to return a list of seemingly non-related numbers in a certain range. It is supposed to be used just once so it does not matter that subsequent call to it produce the same list.

Or maybe there are better approaches that did not occur to me?

没有评论: