Haskell Prime-testing confusion

Question

So I'm totally new to Haskell, hope it doesn't show too much. By any means, I was trying to create a function to determine if a number is prime. The basic idea is this: Given a number, see if it is divisible by any other number less than it. If so, return false. If not, it's prime, return true. The code so far (that is known to be working) is this:

divisible :: Integer -> Integer -> Bool
divisible x 2 = even x
divisible x y = ((x `mod` y) /= 0) && not (divisible x (y-1))
isPrime :: Integer -> Bool
isPrime x = not (even x) && not (divisible x (x-1))

Produces:

ghci> isPrime 9
False
ghci> isPrime 13
True

What I'd like to do is optimize this a bit, since I only need to check values less than or equal to sqrt(x). The problem is, when I try and implement this, stuff gets crazy:

isPrime x = not (even x) && not (divisible x (ceiling(sqrt(fromIntegral(x-1)))))

Besides the fact that it looks terrible (I said I was new), it doesn't give the correct result:

ghci> isPrime 9
False
ghci> isPrime 13
False

I'm trying to figure out what changed, because:

ghci> ceiling(sqrt(13))
4

Seems to be giving me the correct number. I know this is a small problem, but I'm seriously confused...

Answer 1

You got your conditions mixed up:

divisible x y = ((x `mod` y) /= 0) && not (divisible x (y-1))

should be

divisible x y = (x `mod` y) == 0 || divisible x (y-1)

for the test to work.

As is, your divisible function expands e.g.

divisible 21 5 = (21 `mod` 5 /= 0) && not (divisible 21 4)
               = (21 `mod` 5 /= 0) && not ((21 `mod` 4 /= 0) && not (divisible 21 3))
               = not ((21 `mod` 4 /= 0) && not ((21 `mod` 3 /= 0) && not (divisible 21 2)))
               = not (True && not (False && not (divisible 21 3)))
               = not (not False)
               = False

since 21 `mod` 3 == 0, and isPrime 21 evaluates to True with the square root bound.

However, I get

*StrangePrime> isPrime 9
True
*StrangePrime> isPrime 13
True

with your code using the square root.

Without the square root, it happened to work for odd numbers, because the difference between an odd composite and any of its divisors is always even. Unfolding divisible a few steps for n = p*m, where p is the smallest prime factor of the odd composite n, we see

divisible n (n-1) = n `mod` (n-1) /= 0 && not (divisible n (n-2))
                  = not (divisible n (n-2))
                  = not (n `mod` (n-2) /= 0 && not (divisible n (n-3)))
                  = not (not (divisible n (n-3)))
                  = not^2 (divisible n (n-3))

and inductively

divisible n (n-1) = not^(k-1) (divisible n (n-k))

if there are no divisors of n larger than n-k. Now, in the above situation, the largest divisor of n is m = n - (p-1)*m, so we obtain

divisible n (n-1) = not^((p-1)*m-1) (divisible n m)
                  = not^((p-1)*m-1) (n `mod` m /= 0 && not (...))

But n `mod` m == 0, so we have

divisible n (n-1) = not^((p-1)*m-1) False

Since p is odd, p-1 is even, and hence so is (p-1)*m, so altogether we have an odd number of nots, which is the same as one not, giving

divisible n (n-1) = True
isPrime n = not (even n) && not (divisible n (n-1)) = True && not True = False

If p is an odd prime, the unfolding reaches divisible p (p-1) = not^(p-3) (divisible p (p - (p-2))). p-3 is even, divisible p 2 is even p, which is False.

Generally, consider divisible n s for an odd n, and let d be the largest divisor of n not exceeding s, if n is composite, or d = 2 if n is prime. The unfolding of divisible n s still proceeds the same way

divisible n s = not^k (divisible n (s-k))

while no divisor has been found and s-k > 2. So in the case of a composite n, we find

divisible n s = not^(s-d) (divisible n d)
              = not^(s-d) (n `mod` d /= 0 && not (...))
              = not^(s-d) False
              = odd (s-d)
              = even s     -- since d is odd, as a divisor of an odd number

and in the case of an odd prime n,

divisible n s = not^(s-2) (divisible n 2)
              = not^(s-2) (even n)
              = not^(s-2) False
              = odd s

So divisible n s measures the parity of the distance of s to the next smaller divisor of n or to 2, whichever is larger. When s was n-1, the starting point was always even, so it worked out correctly, but ceiling (sqrt (fromIntegral (n-1))) can be odd, in which case the results are flipped and composites are declared prime and vice versa.

You can make your divisible function work for the primality test of odd numbers with a square root bound if you make sure that the first call gets an even second argument (so if ceiling (sqrt (fromIntegral (n-1))) is odd, start at ceiling (sqrt (fromIntegral (n-1))) + 1), but the logic of that function is confusing, and its name doesn't correctly describe its results.

A more idiomatic way to write it would be

isPrime n = and [n `mod` k /= 0 | k <- [2 .. ceiling (sqrt $ fromIntegral n)]]

The test becomes more efficient when one skips candidate divisors that are already known to be nondivisors from prior tests, easy is skipping all even numbers except 2,

isPrime 2 = True
isPrime n = all ((/= 0) . (n `mod`)) (2 : [3, 5 .. ceiling (sqrt (fromIntegral n))])

slightly more involved, but still more efficient is also skipping multiples of 3

isPrime n = all ((/= 0) . (n `mod`)) (takeWhile (<= bound) (2:3:scanl (+) 5 (cycle [2,4])))
  where
    bound = ceiling (sqrt (fromIntegral (n-1)))

In the same vein one can eliminate the multiples of more small primes from the trial divisors, each gaining a bit of efficiency, but at the cost of a more complicated wheel, e.g. also eliminating multiples of 5 leads to

isPrime n = all ((/= 0) . (n `mod`)) (takeWhile (<= bound) (2:3:5: scanl (+) 7 (cycle [4,2,4,2,4,6,2,6])))
  where
    bound = ceiling (sqrt (fromIntegral (n-1)))

Answer 2

Here's how I'd do it:

divisibleBy :: (Integral a) => a -> a -> Bool
divisibleBy x y = mod x y == 0

isPrime :: (Integral a) => a -> Bool
isPrime x = or $ map (divisibleBy x) [2..(x-1)]

divisibleBy is a simple test of divisibility. isPrime performs this test on all integers between 1 and x, returning true if x is divisible by any of those integers. You might change the upper bound to root x, as you've done in your code, but otherwise this works.