6
\$\begingroup\$

Inspired by a tweet linked to me by a friend and a Haskell implementation by her for the same problem, I decided to try my hand at approximating the value of π using everything in the Haskell standard library I could find for the job. Here’s what I came up with:

module Pi where

import Data.List (genericLength)
import Control.Arrow (Arrow, (<<<), (***), arr)
import System.Random (newStdGen, randoms)

type Point a = (a, a)

chunk2 :: [a] -> [(a, a)]
chunk2 []      = []
chunk2 [_]     = error "list of uneven length"
chunk2 (x:y:r) = (x, y) : chunk2 r

both :: Arrow arr => arr a b -> arr (a, a) (b, b)
both f = f *** f

unsplit :: Arrow arr => (a -> b -> c) -> arr (a, b) c
unsplit = arr . uncurry

randomFloats :: IO [Float]
randomFloats = randoms <$> newStdGen

randomPoints :: IO [Point Float]
randomPoints = chunk2 <$> randomFloats

isInUnitCircle :: (Floating a, Ord a) => Point a -> Bool
isInUnitCircle (x, y) = x' + y' < 0.25
  where x' = (x - 0.5) ** 2
        y' = (y - 0.5) ** 2

lengthRatio :: (Fractional c) => [b] -> [b] -> c
lengthRatio = curry (unsplit (/) <<< both genericLength)

approximatePi :: [Point Float] -> Float
approximatePi points = circleRatio * 4.0
  where circlePoints = filter isInUnitCircle points
        circleRatio = circlePoints `lengthRatio` points

main :: IO ()
main = do
  putStrLn "How many points do you want to generate to approximate π?"
  numPoints <- read <$> getLine
  points <- take numPoints <$> randomPoints
  print $ approximatePi points

I’m interested in a general review, but I’m especially curious about my use of arrows: is there a better way to write lengthRatio? Are anything like both and unsplit provided anywhere in the standard library? If not, do any packages help?

\$\endgroup\$
5
\$\begingroup\$

About those arrows

I’m interested in a general review, but I’m especially curious about my use of arrows: is there a better way to write lengthRatio?

Compare the following two lines. Both do the same, but which one would you rather see if you need to change your code drunk in three months, with only 5% battery left?

lengthRatio = curry (unsplit (/) <<< both genericLength)

lengthRatio xs ys = genericLength xs / genericLength ys

Also, which one has which type, and which one is more general?

Arrows are great if you want to abstract functions. But throughout your small script, you're still just working with (->), not any other instance of Arrow. For a small script like this, Arrow is too much. For example the pointwise definition above is actually a character shorter than the pointfree one. Sure, the pointfree one is clever, but it's also very beginner-unfriendly.

About randomness

randomPoints introduces a dependency between your point coordinates \$x\$ and \$y\$, since both draw from the same sequence. This usually leads to points on hyperplanes (see disadvantages of LCG and spectral test). Your friends variant doesn't have this immediate problem:

randomTuples :: Int -> IO [(Float, Float)]
randomTuples n = do
  seed1 <- newStdGen
  seed2 <- newStdGen
  let xs = randoms seed1 :: [Float] -- two different
      ys = randoms seed2 :: [Float] -- generators being used
  return $ take n $ zipWith (,) xs ys

However, since newStdGen is merely a split, it's more or less hiding the dependency at another place. Still, it's something to keep in mind, if you don't want to end up with something like this.

But how would you check this? Well, you would run tests, over and over. Here's the second design critique on randomPoints, it doesn't take a RandomGen. Truth be told, if I say that Arrow is too much for a small script, then

randomPoints :: RandomGen g => g -> [Point Float]

is too much either.

Also, if you know you're going to generate Points, a newtype Point a together with

instance Random a => Random (Point a) where 

is feasible and doesn't introduce a potential error via chunk2. Keep possible problems with Random in mind, though.

About names

The function isInUnitCircle lies. It's not testing whether the point \$(x,y)\$ lies in the circle with radius \$r = 1\$ with center in the origin, e.g. $$ \sqrt{x^2 + y^2} \le 1^2 \Leftrightarrow x^2 + y^2 \le 1 $$ but in the circle with diameter \$d = 2r = 1\$ with center in \$(0.5, 0.5)\$. In the following picture, the green region is where you generate your random values. In the left one, you see the regular unit circle, in the right one, you see the circle size you're actually testing (after shifting your values from the green square into the red one):

enter image description here

Therefore, you're not calculating the "usual" fourth of a circle, but instead a circle with a fourth of the original size (\$\pi(\frac{1}{2})^2 = \frac{\pi}{4}\$)). Luckily, it doesn't matter for the convergence.

A real test that checks whether a point is in the unit circle is tremendously easier:

isInUnitCircle :: (Num a, Ord a) => Point a -> Bool
isInUnitCircle (x, y) = x ^ 2 + y ^ 2 <= 1

About optimization

Last, but not least, there's an issue with approximatePi, or rather the use of lengthRatio on the same list twice. Actually, taking the length of the list again is a litle bit strange, since you know how large the sample is:

numPoints <- read <$> getLine                -- sample size
points    <- take numPoints <$> randomPoints
print $ approximatePi points                 -- sample size still known (?)

But let's say that you don't actually know how many points you have. Let's assume that someone wants to check a many points. Suddenly, the memory usage of your program explodes:

$ echo 10000000 | ./CalcPi +RTS -s
How many points do you want to generate to approximate π?
3.141744
  33,724,505,920 bytes allocated in the heap
   5,288,250,096 bytes copied during GC
   1,319,621,976 bytes maximum residency (17 sample(s))
       5,554,344 bytes maximum slop
            2587 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0     63626 colls,     0 par    1.606s   3.155s     0.0000s    0.0006s
  Gen  1        17 colls,     0 par    2.732s   3.373s     0.1984s    1.2720s

  INIT    time    0.000s  (  0.000s elapsed)
  MUT     time   17.158s  ( 15.542s elapsed)
  GC      time    4.337s  (  6.528s elapsed)
  EXIT    time    0.019s  (  0.166s elapsed)
  Total   time   21.514s  ( 22.236s elapsed)

  %GC     time      20.2%  (29.4% elapsed)

  Alloc rate    1,965,530,983 bytes per MUT second

  Productivity  79.8% of total user, 77.2% of total elapsed

Even though randoms generates a lazy list, approximatePi needs to hold onto it completely due to lengthRatio. A classic space leak. The altnerative version of lengthRatio won't save you from that. Instead, provide a function to check the ratio of filtered elements:

-- Rational from Data.Ratio
filterRatio :: (a -> Bool) -> [a] -> Rational
filterRatio p xs = -- exercise

That way, you can define a version of approximatePi that works for large lists:

approximatePi :: [Points Float] -> Double
approximatePi points = circleRatio * 4
  where 
    circleRatio = fromRational $ filterRatio isInUnitCircle points
$ echo 10000000 | ./GenPIRatio +RTS -s
How many points do you want to generate to approximate π?
3.1421592
  24,445,866,792 bytes allocated in the heap
      15,555,552 bytes copied during GC
          77,896 bytes maximum residency (2 sample(s))
          21,224 bytes maximum slop
               1 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0     46874 colls,     0 par    0.012s   0.113s     0.0000s    0.0001s
  Gen  1         2 colls,     0 par    0.000s   0.000s     0.0001s    0.0001s

  INIT    time    0.000s  (  0.000s elapsed)
  MUT     time   10.809s  ( 10.746s elapsed)
  GC      time    0.012s  (  0.113s elapsed)
  EXIT    time    0.000s  (  0.000s elapsed)
  Total   time   10.821s  ( 10.859s elapsed)

  %GC     time       0.1%  (1.0% elapsed)

  Alloc rate    2,261,657,279 bytes per MUT second

  Productivity  99.9% of total user, 99.5% of total elapsed

Summary

Food for thought:

  • Use the right level of abstraction for your problem. Arrow is an overkill for such a small script, but alright for learning.
  • Try to decrease the amount of IO wherever possible, but again, that might be too abstract for a small script.

Bad:

  • Don't lie, give things the right name.
  • Don't overcomplicate, keep pointfree to a sane minimum.
  • Major space leak in approximatePi. Read the linked section of RWH and try to define filterRatio or a similar function.

Good:

  • Type signatures! Yay!
  • Explicit imports!
  • Type synonym instead of (a, a) everywhere!

So beside the slight arrow-overkill, well done.

\$\endgroup\$
  • \$\begingroup\$ Thanks a lot for the thorough answer! I agree that arrows are pretty heavily overkill here, though I find it a little disappointing there isn’t a simpler way to represent that sort of computation in a point-free way (one that’s actually readable). It would seem arrows aren’t terribly useful for things besides things like FRP, from what I’ve seen. Also, thanks for the tip for filterRatio, though I’m not 100% sure what you have in mind... do you just recommend folding and summing? In that case, wouldn’t that just be length . filter, hoping stream fusion takes care of things? \$\endgroup\$ – Alexis King Jan 21 '16 at 22:57
  • \$\begingroup\$ length . filter only tells you how many things passed the predicate, not how many things were originally in there. If you want to get a ratio, you need the original length. Either calculate the length of the list along the number of elements that passed the predicate, or add Int -> to approximatePi. But stream fusion cannot get rid of the space leak if you use the same list twice. \$\endgroup\$ – Zeta Jan 22 '16 at 5:06
  • \$\begingroup\$ Ahh, I see what you’re getting at, that makes sense. Thanks for pointing that out—my ability to inuit about laziness is not as robust as I’d like. \$\endgroup\$ – Alexis King Jan 23 '16 at 3:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.