Sqrt (square root) function

Question

As an exercise in learning Scala, I implemented a square root function like this:

  def sqrt(x: Double): Double = {
    if (x < 0) throw new IllegalArgumentException("negative numbers not allowed")
    val threshold = if (x < 1) x / 1e15 else 1e-12
    def sqrt(x: Double, p: Double): Double = {
      if (p == x / p) p  // without this condition, non-termination with 1e50
      else if (Math.abs(p * p - x) < threshold) {
        def diff1 = Math.abs(x - p * p)
        def diff2 = Math.abs(x - x / p * x / p)
        if (diff1 < diff2) p else x / p
      }
      else sqrt(x, (p + x / p) / 2)
    }
    sqrt(x, x / 2)
  }

The implementation passes these unit tests:

  test("sqrt 2") {
    assert(sqrt(2) === Math.sqrt(2))
  }

  test("sqrt 1e-3") {
    assert(Math.abs(1e-3 - sqrt(1e-3) * sqrt(1e-3)) ===
      Math.abs(1e-3 - Math.sqrt(1e-3) * Math.sqrt(1e-3)))
  }

  test("sqrt 1e-20") {
    assert(sqrt(1e-20) === Math.sqrt(1e-20))
  }

  test("sqrt 1e-21") {
    assert(sqrt(1e-21) === Math.sqrt(1e-21))
  }

  test("sqrt 1e20") {
    assert(sqrt(1e20) === Math.sqrt(1e20))
  }

  test("sqrt 1e50") {
    assert(sqrt(1e50) === Math.sqrt(1e50))
  }

My questions:

1. How can I improve this?

2. Notice the unit test for the case of 1e-3. It's more complex than the others to compensate for the difference between sqrt and Math.sqrt. Although the both sqrt and Math.sqrt are equally incorrect (the square of both have the same discrepancy with x), I wonder if I can change the implementation to match the result of Math.sqrt.

3. I found the error thresholds for < 1 and >= 1 through trial and error: all unit tests pass with these values and some would fell if I stretched the limits further. I'm wondering if there is a better, proper way of setting suitable thresholds based on x and the numeric limits of the language.

Answer 1

I'm going to concentrate on the Scala rather than the maths side, because that's my area.

The internal helper function doesn't need two parameters, since the value of x never changes; pushing x repeatedly onto the stack is a waste of time. Also, better to use a case statement than a chain of if...else if...else - much less error prone.

def sqrt(x: Double): Double = {
  if (x < 0) throw new IllegalArgumentException("negative numbers not allowed")
  val threshold = if (x < 1) x / 1e15 else 1e-12
  def sqrtx(p: Double): Double = p match {
    case q if q == x / q => q // without this condition, non-termination with 1e50
    case q if Math.abs(q * q - x) < threshold => {
      def diff1 = Math.abs(x - p * p)
      def diff2 = Math.abs(x - x / p * x / p)
      if (diff1 < diff2) p else x / p
    }
    case _ => sqrtx((p + x / p) / 2)
  }
  sqrtx(x / 2)
}

I renamed the inner function sqrtx to avoid confusion about what it is doing.

Answer 2

It took a bit to remember where I'd seen this before. Implementing Newton's Method for taking a Square Root was one of the exercises in Odersky's Functional Programming in Scala course on Coursera.

Here is an implementation that's mostly based on the hints and suggestions in that assignment.

object NewtonsMethod {
  def sqrt(x: Double) = {
    require(x >= 0, "Square Root is undefined for negative numbers")

    def isGoodEnough(guess: Double) = Math.abs(guess * guess - x) / x < 0.001
    def improve(guess: Double) = (guess + x / guess) / 2

    @tailrec
    def iter(guess: Double): Double =
      if(isGoodEnough(guess)) guess
      else iter(improve(guess))

    iter(1.0)
  }
}

It has some very nice features, one of the most important is the helper function isGoodEnough, which scales the threshold with the size of x. This makes the non-termination check for 1e50 unnecessary.

It's tail-recursive, and flagged as such, honestly it shouldn't make a difference for this algorithm, but it's a nice touch.

Here's the test suite I used to verify it has the correct behavior.

import org.scalatest.WordSpec
import org.scalatest.Matchers
import org.scalacheck.Prop
import org.scalatest.prop.GeneratorDrivenPropertyChecks
import NewtonsMethodSpec._

class NewtonsMethodSpec extends WordSpec with Matchers with GeneratorDrivenPropertyChecks {
  "NewtonsMethod.sqrt" should {
    "work for very large numbers" in {
      NewtonsMethod.sqrt(1e50).round(2) shouldBe Math.sqrt(1e50).round(2)
    }
    "be equivalent to Math.sqrt" in {
      forAll {
        (d: Double) => whenever(d >= 0) {
          NewtonsMethod.sqrt(d).round(2) shouldBe Math.sqrt(d).round(2)
        }
      }
    }
  }
}

object NewtonsMethodSpec {
  implicit class Roundable(val d: Double) extends AnyVal {
    def round(scale: Int): BigDecimal =
      BigDecimal(d).round(new java.math.MathContext(scale, java.math.RoundingMode.FLOOR))
  }
}