One of my favorite ways to sharpen my skills is to implement algorithms. There's something about transferring an idea from words and equations into code that cements the knowledge of how it works in my brain. It's also a way to avoid cheating yourself into believing that you understand something just because you read it — if you can translate an algorithm into code and execute it, then you know you really understand it.

In this article I'm going to implement two algorithms in Ruby — the Euclidean algorithm for finding the greatest common divisor of two integers, and the Chakravala method for solving Pell's equation.

The Euclidean algorithm

The Euclidean algorithm is an extremely straightforward way to calculate the greatest common divisor of two natural numbers. It works because of the observation by Euclid around 300 BCE that a divisor of two numbers will also divide evenly into their difference.

def euclid(x, y)

  # make sure that x is larger than y
  if x < y 
    x, y = y, x
  end

  loop do 
    r = x % y
    if r.zero?
      return y
    else
      x, y = y, r
    end
  end

end

The Euclidean algorithm has two steps:

1. Divide x by y and examine the remainder r.
2. If r == 0, then y is the GCD. If not, let x = y and y = r and repeat.

Because all natural numbers are divisible by one, the algorithm will always result in an answer.

Not only is the Euclidean algorithm simple and elegant, it also wins a spot as one of the oldest algorithms still in common use — there are faster methods to calculating the GCD of two extremely large numbers, but for most cases the Euclidean algorithm is sufficient, and thus is still being used today.

The Chakravala method

The Chakravala method is a gorgeous, easily computable algorithm for finding the minimal integer solution for indeterminate quadratic equations, including Pell's equation x^2 – Dy^2 = 1. It was published in the 12th century, about 600 years before similar problems were solved in Europe, and was called “the finest thing achieved in the theory of numbers before Lagrange” by Hermann Hankel.

It also happens to be remarkably simple to implement, though I found the explanation in its Wikipedia entry pretty confusing. Some cursory googling lead me this explanation, though, which was much clearer, and is what the code below is based upon (though in this version we don't bother with the trivial solution as a starting point).

# returns the array [x, y] that solves the equation x^2 = Ny^2 + 1 for a given N
def chakravala(n)

  raise ArgumentError, "n cannot be a square!" if Math.sqrt(n) % 1 == 0

  current_p = current_x = optimal_p = Math.sqrt(n).round
  current_k = (current_p ** 2) - n

  current_y = 1

  # puts "#{current_p}, #{current_k}, #{current_x}, #{current_y}"

  return [current_x, current_y] if current_k == 1 # simplest answer was right!

  loop do 

    current_k_abs   = current_k.abs
    diff            = (current_p + optimal_p) % current_k_abs
    next_p_low      = optimal_p - diff
    next_p_high     = next_p_low + current_k_abs

    if next_p_low < 1
      next_p = next_p_high
    else
      next_p = ((next_p_low ** 2) - n).abs < ((next_p_high ** 2) - n).abs ? next_p_low : next_p_high
    end

    next_k = ((next_p ** 2) - n) / current_k
    next_x = ((next_p * current_x) + (n * current_y)) / current_k_abs
    next_y = ((next_p * current_y) + current_x) / current_k_abs

    # puts "#{next_p}, #{next_k}, #{next_x}, #{next_y}"
    
    return [next_x, next_y] if next_k == 1

    current_p = next_p
    current_k = next_k
    current_x = next_x
    current_y = next_y

  end
end

The algorithm has just three steps, and the last two of which are repeated until a solution is found (and it is proven that a solution will be found for all non-square values of N).

1. Choose a value for current_p such that the absolute value of current_p ** 2 - n is minimal (as close to zero as possible). This is most easily done by rounding the square root of n and using that for current_p.
   • From that, it follows that current_k = (current_p ** 2) - n. current_x = current_p at this point, and current_y = 1.
   • If current_k == 1 then you're done and you can return immediately (this is the case for some values of n, such as n = 8).
2. Choose a value for next_p such that (current_p + next_p) % current_k.abs == 0 and the absolute value of next_p ** 2 - n is minimal. Again rounding the square root of n gives the best choice to satisfy the second requirement, but as we have to ensure that current_p + next_p can be divided evenly by the absolute value of current_k, we likely can't use it directly. Instead, we need to test two values closest to the square root of n (which in the code we call optimal_p) that meet the first criteria, and then decide which one to use.
   • First, find out how far away optimal_p is from a valid value for next_p. Do this by getting the remainder of (current_p + optimal_p) / current_k_abs (that is, (current_p + optimal_p) % current_k_abs) where current_k_abs is just the absolute value of current_k, and storing it as diff for use in the next calculation.
   • The first potential value for next_p is diff less than optimal_p. We'll call that next_p_low. This makes sense because if optimal_p was diff lower, then there would have been no remainder on (current_p + optimal_p) / current_k_abs.
   • The second potential value for next_p is just next_p_low + current_k_abs, which is greater than optimal_p, and which we'll call next_p_high.
   • If next_p_low is less than one, then automatically choose next_p_high. Otherwise, choose the one that results in the smallest absolute value when used as p in the equation (p ** 2) - n.
3. Finally, calculate new values for k, x, and y. If the new value of k equals 1, return. Otherwise, loop back to Step 2 with updated values and do it all over again.
   • next_k = ((next_p ** 2) - n) / current_k
   • next_x = ((next_p * current_x) + (n * current_y)) / current_k_abs
   • next_y = ((next_p * current_y) + current_x) / current_k_abs

And that's it. For any non-square n, the algorithm will eventually return the minimal integer solution. Perhaps more importantly, the algorithm will, if you let it go, return every integer solution (of which there are infinitely many). Even more interestingly, the values that it returns are a good rational approximation of the square root of n, with larger results being increasingly accurate. This is helpful if you need square roots but want to avoid floating point math for whatever reason (and there are plenty of reasons!).

Here are a few examples with intermediate values to help you verify your work if you're trying to implement this algorithm yourself.

Now what?

What's the point of implementing algorithms like this? I think there are a number of benefits:

1. You can learn / relearn bits of math that you may have forgotten (or never learned in the first place).
2. You can sharpen your programming skills (many mathematical algorithms rely on dynamic programming or recursion, both skills that you need to have in your toolbox).
3. You can take insights away from them and use those insights to solve future problems.

In the case of the Euclidean algorithm it really is reinventing the wheel — Ruby's Integer class has a GCD method that is probably implemented in C and is quite a bit faster than our pure Ruby approach — but the Chakravala method solves a problem that most programming languages can't solve right out of the box, and in an elegant way.

Plus, it's fun, and that's what's really important.