# Haskell code optimization for short Pi approximator

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.

• There seems to be some review context missing here. Maybe you defined r as some large integer, ten raised to roughly the number of decimal places to which you'd like to find pi ? From your "counting points within a circle" description I immediately went to rolling a pair of random numbers in the unit interval and asking whether they land within (quadrant I of) the unit circle. But across integers?!? I'm not yet seeing it. Anyway, there's never a need for a time consuming sqrt() to produce bool return value, since you can stick to comparing squared quantities. Just square both sides.
– J_H
Oct 30 at 3:54
• @J_H You understood right. Unless arithmetic gets slower with big numbers, working with integers should yield the same result as with a unit circle. The observation about sqrt is spot on and I wonder how I didn't think about it earlier. Oct 30 at 8:13
• Examining benchmark timings would be entertaining. My gut says to stick with floats on the unit interval: simple to understand, and {fast, compact} to compute. Python's bignums (int) are wonderful and they'll go forever, for as long as you've got RAM to malloc. But they are slightly weird under the hood. Count to a million? No problem, very simple representation. Use numbers in the vicinity of 2 ** 53, comparable to FP 53 bits of significand? Now we have spilled into two objects with max of 2 ** 30. Back in 2008 it would have been four, with max of 2 ** 15. Carry and such propagates.
– J_H
Oct 30 at 16:07

### 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;
}