Monadic approach to estimating PI in scala

I'm trying to understand how to leverage monads in scala to solve simple problems as way of building up my familiarity. One simple problem is estimating PI using a functional random number generator. I'm including the code below for a simple stream based approach.

I'm looking for help in translating this to a monadic approach. For example, is there an idiomatic way convert this code to using the state (and other monads) in a stack safe way?

trait RNG {
    def nextInt: (Int, RNG)
    def nextDouble: (Double, RNG)
}

case class Point(x: Double, y: Double) {
    val isInCircle = (x * x + y * y) < 1.0
}

object RNG {
    def nonNegativeInt(rng: RNG): (Int, RNG) = {
      val (ni, rng2) = rng.nextInt
      if (ni > 0) (ni, rng2)
      else if (ni == Int.MinValue) (0, rng2)
      else (ni + Int.MaxValue, rng2)
    }

    def double(rng: RNG): (Double, RNG) = {
      val (ni, rng2) = nonNegativeInt(rng)
      (ni.toDouble / Int.MaxValue, rng2)
    }


    case class Simple(seed: Long) extends RNG {
      def nextInt: (Int, RNG) = {
      val newSeed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL
      val nextRNG = Simple(newSeed)
      val n = (newSeed >>> 16).toInt
      (n, nextRNG)
    }

    def nextDouble: (Double, RNG) = {
      val (n, nextRNG) = nextInt
      double(nextRNG)
    }
  }
}

object PI {
    import RNG._

    def doubleStream(rng: Simple):Stream[Double] = rng.nextDouble match {
        case (d:Double, next:Simple) => d #:: doubleStream(next)
    }

    def estimate(rng: Simple, iter: Int): Double = {
        val doubles = doubleStream(rng).take(iter)
        val inside = (doubles zip doubles.drop(3))
            .map { case (a, b) => Point(a, b) }
            .filter(p => p.isInCircle)
            .size * 1.0
        (inside / iter) * 4.0
    }
}

// > PI.estimate(RNG.Simple(10), 100000)
// res1: Double = 3.14944

I suspect I'm looking for something like replicateM from the Applicative monad in cats but I'm not sure how to line up the types or how to do it in a way that doesn't accumulate intermediate results in memory. Or, is there a way to do it with a for comprehension that can iteratively build up Points?

Solution

Id you want to iterate using monad in a stack safe way, then there is a tailRecM method implemented in Monad type class:

// assuming random generated [-1.0,1.0]
def calculatePi[F[_]](iterations: Int)
                     (random: => F[Double])
                     (implicit F: Monad[F]): F[Double] = {
  case class Iterations(total: Int, inCircle: Int)
  def step(data: Iterations): F[Either[Iterations, Double]] = for {
    x <- random
    y <- random
    isInCircle = (x * x + y * y) < 1.0
    newTotal = data.total + 1
    newInCircle = data.inCircle + (if (isInCircle) 1 else 0)
  } yield {
    if (newTotal >= iterations) Right(newInCircle.toDouble / newTotal.toDouble * 4.0)
    else Left(Iterations(newTotal, newInCircle))
  }
  // iterates until Right value is returned
  F.tailRecM(Iterations(0, 0))(step)
}
calculatePi(10000)(Future { Random.nextDouble }).onComplete(println)

It uses by-name param because you could try to pass there something like Future (even though the Future is not lawful), which are eager, so you would end up with evaluating the same thing time and time again. With by name param at least you have the chance of passing there a recipe for side-effecting random. Of course, if we use Option, List as a monad holding our "random" number, we should also expect funny results.

The correct solution would be using something that ensures that this F[A] is lazily evaluated, and any side effect inside is evaluated each time you need a value from inside. For that you basically have to use some of Effects type classes, like e.g. Sync from Cats Effects.

def calculatePi[F[_]](iterations: Int)
                     (random: F[Double])
                     (implicit F: Sync[F]): F[Double] = {
  ...
}
calculatePi(10000)(Coeval( Random.nextDouble )).value
calculatePi(10000)(Task( Random.nextDouble )).runAsync

Alternatively, if you don't care about purity that much, you could pass side effecting function or object instead of F[Int] for generating random numbers.

// simplified, hardcoded F=Coeval
def calculatePi(iterations: Int)
               (random: () => Double): Double = {
  case class Iterations(total: Int, inCircle: Int)
  def step(data: Iterations) = Coeval {
    val x = random()
    val y = random()
    val isInCircle = (x * x + y * y) < 1.0
    val newTotal = data.total + 1
    val newInCircle = data.inCircle + (if (isInCircle) 1 else 0)
    if (newTotal >= iterations) Right(newInCircle.toDouble / newTotal.toDouble * 4.0)
    else Left(Iterations(newTotal, newInCircle))
  }
  Monad[Coeval].tailRecM(Iterations(0, 0))(step).value
}