Problem 066: (source):

import scala.collection.mutable.{Map => MMap}
object Euler066 extends EulerApp {
  def leastx (d :Int) = {
    val amap = MMap[Int,BigInt]()
    val pmap = MMap[Int,BigInt]()
    val qmap = MMap[Int,BigInt]()
    val ppmap = MMap[Int,BigInt]()
    val qqmap = MMap[Int,BigInt]()

    def a (i :Int) :BigInt = amap.getOrElseUpdate(i, i match {
      case 0 => math.sqrt(d).toInt
      case _ => (a(0) + pp(i))/qq(i)
    })

    def p (i :Int) :BigInt = pmap.getOrElseUpdate(i, i match {
      case 0 => a(0)
      case 1 => a(0)*a(1) + 1
      case _ => a(i)*p(i-1) + p(i-2)
    })

    def q (i :Int) :BigInt = qmap.getOrElseUpdate(i, i match {
      case 0 => 1
      case 1 => a(1)
      case _ => a(i)*q(i-1) + q(i-2)
    })

    def pp (i :Int) :BigInt = ppmap.getOrElseUpdate(i, i match {
      case 0 => 0
      case 1 => a(0)
      case _ => a(i-1)*qq(i-1) - pp(i-1)
    })

    def qq (i :Int) :BigInt = qqmap.getOrElseUpdate(i, i match {
      case 0 => 1
      case 1 => BigInt(d) - a(0)*a(0)
      case _ => (BigInt(d) - pp(i)*pp(i))/qq(i-1)
    })

    val as = Stream.from(0) map(a)
    val r = (as.tail takeWhile(ai => ai != as.head*2) length)
    if (r % 2 == 1) p(r) else p(2*r+1)
  }

  def square (i :Int) = i*i
  def issquare (d :Int) = square(math.sqrt(d).toInt) == d
  def answer = (2 to 1000 filterNot(issquare) map(d => (d, leastx(d))) sortBy(_._2) last)._1
}

My efforts to post my Euler solutions with commentary remain woefully behind (I just finished problem 110), but when I took action to close this gap, I discovered that I had skipped over problem 66 so long ago. I recall being frustrated after beating my head against it for a while and making no headway.

Filled with renewed vigor, I tackled it again here under the invigorating Mexican sunshine. Alas, the sunshine proved to be of little help. I spent a long time wrangling with the equation in the form y² = (x+1)(x-1)/D. This led me to explore the family of solutions nD + 1 and nD – 1, which are very nice, but definitely not minimal.

I then played around with factoring D and trying to come up with some way to efficiently search for potential solutions using those factors, since D has to divide evenly into (x + 1)(x – 1). But it needn’t do so in any convenient way. For example, for D = 28, the minimal solution is (127 + 1)(127 – 1)/28 which factors into 2 · 2 · 2 · 2 · 2 · 2 · 2 × 2 · 3 · 3 · 7 / 7 · 2 · 2. Some of the factors go into the left numerator and some in the right, it’s mayhem.

Eventually I decided to do some Googling and discovered that this is closely related to the well-known Pell equation. For reasons not entirely clear to me, its solutions lie in the continued fraction expansion of the square root of D.

Implementing the “usual recurrence relations" was pretty straightforward, though a bit of memoization was needed (the unsightliness of that part of the code annoys me, but I am annoyed enough with this problem already that I don’t feel like twiddling it further). BigInt also turned out to be necessary, as the largest “minimal" x is 16421658242965910275055840472270471049, which doesn’t quite fit into 64 bits.