Haskell code optimization for short Pi approximator

Question

Very new in Haskell, as a first easy program I went for an old algorithm of mine to approximate pi by counting points within a circle. The snippet below is what I could get working. I had quite an issue with types (you will notice very extensive use of fromIntegral) and would also like to know about performance, since I read that Haskell can be quite fast (the same algorithm I test in C++ and it is a lot faster).

So to sum up:

1. How could this code be written better/shorter in Haskell?

2. How would this code be written for better performance?

    module PiCalc (approximatePi) where
   
    type Point = (Int,Int)
   
    pointInCircle :: Point -> Int -> Bool
    pointInCircle (x,y) r =
        let xd = fromIntegral x
            yd = fromIntegral y in
        sqrt (xd*xd + yd*yd) <= fromIntegral r
   
    countPoinstInDistance :: Int -> Int
    countPoinstInDistance n = sum [1 | x <- [0..n], y <- [0..n], pointInCircle(x,y) n]
   
    approximatePi :: Int -> Double
    approximatePi n = 4 * fromIntegral (countPoinstInDistance n) / fromIntegral (n*n)
   

The algorithm works by passing in a large number like 10000 for example and checking the distance of all points (x,y) from the center.

Answer 1

Incorrect results

There's a bug in your code. It produces incorrect results. The line

[1 | x <- [0..n], y <- [0..n], pointInCircle(x,y) n]

should be

[1 | x <- [0..n], y <- [1..n], pointInCircle(x,y) n]

so that the x=0 column is included in the count but the y=0 row is not, because it's the same pixels as the x=0 column of the adjacent quarter circle. If you include both in the tally, you effectively count the border pixels twice. Your code returns 3.14199056 for approximatePi 10000 and with the bug fixed, it returns 3.14159052.

Style and Performance

Instead of clobbering the namespace with top level functions and types, define functions in where clauses. I removed the Point type in favor of taking separate x and y arguments and I put the pointInCircle function into a where clause in countPointsInDistance. Note how pointInCircle now captures n from the parent scope so that you only have to pass x and y as arguments. I also squared the right side of the comparison instead of taking the square root on the left side, which removes floating point calculations (and fromIntegral conversions) altogether.

countPointsInDistance :: Int -> Int
countPointsInDistance n = sum [1 | x <- [0..n], y <- [1..n], pointInCircle x y]
  where pointInCircle x y = x*x + y*y <= n*n

approximatePi :: Int -> Double
approximatePi n = 4 * fromIntegral (countPointsInDistance n) / fromIntegral (n*n)

You can take this one step further by removing countPointsInDistance as a top-level function as well, but you probably wanted to keep that one accessible for testing. If you want to remove it, your entire code simplifies down to:

approximatePi :: Int -> Double
approximatePi n = 4 * fromIntegral pointsInDistance / fromIntegral (n*n)
  where pointsInDistance = sum [1 | x <- [0..n], y <- [1..n], pointInCircle x y]
        pointInCircle x y = x*x + y*y <= n*n

Haskell can be rather fast, but you have to compile your program into an executable with optimizations turned on. Running code in GHCi is slow as hell. Also keep in mind that Int are 64 bit integers (fast, but can overflow) whereas Integer are unbounded, but slower. Using Integer by default is highly recommended, although it makes performance comparisons harder when the other language, often C++, uses overflowing integers.

Algorithm improvements

As usual, the biggest improvement comes from a better algorithm. Instead of running the inner loop over all possible values of y, take y from the previous iteration of the outer loop and adjust it slightly because x has been incremented. My code starts with x=0, y=n and increments x as part of the outer loop. After incrementing x, it decrements y until the condition x*x + y*y <= n*n is valid again. The remaining decrements of y down to zero are omitted. I just add that number to the counter instead of actually iterating down to zero and adding 1 on each iteration.

countPointsInDistance :: Int -> Int
countPointsInDistance n = go 0 n 0
  where go x y count
          | y == 0           = count
          | x*x + y*y <= n*n = go (x+1)  y  (count+y)
          | otherwise        = go   x  (y-1) count

This code can be improved further by not squaring the numbers, which prevents overflow issues for native integers, or performance issues for arbitrary precision integers. We can do this by introducing a variable z = x*x + y*y - n*n and updating z whenever we increment x or decrement y. We never recalculate z from scratch. We initialize with x=0, y=n, z=0

countPointsInDistance :: Int -> Int
countPointsInDistance n = go 0 n 0 0
  where go x y z count
          | y == 0    = count
          | z <= 0    = go (x+1)  y   (z+2*x+1) (count+y)
          | otherwise = go   x  (y-1) (z-2*y+1)  count

If you want to compare the performance of Haskell vs. C++, the equivalent C++ implementation of the last code sample would be:

int countPointsInDistance(int n) {
    int x = 0, y = n, z = 0, count = 0;
    while (y != 0) {
        if (z <= 0) {
            z += 2*x+1;
            x += 1;
            count += y;
        } else {
            z += -2*y+1;
            y -= 1;
        }
    }
    return count;
}