From: Takeshi NISHIMATSU on
Here is my "Pell's equation solver in Ruby".
Can you make it shorter or faster?

$ time ruby pell.rb
1
2 3 2
3 2 1
4
5 9 4
:
998 984076901 31150410
999 102688615 3248924
1000 39480499 1248483

real 0m5.415s
user 0m5.403s
sys 0m0.011s
$ ruby -v
ruby 1.8.6 (2008-08-11 patchlevel 287) [universal-darwin9.0]
(2.8 GHz Intel Xeon)

=============================================
#!/usr/bin/env ruby
# pell.rb - Pell's equation solver in Ruby
# x**2 - D*y**2 = 1 (D is a natural number but not a square number.) (1)
# c_i=[a_0; a_1, a_2, ..., a_i] is the continued fraction expression of sqrt(D).
# c=y/x will be the minimal integral solution of Eq.(1) for some i>=1.
# omega_0 = sqrt(D)
# a_i = [omega_i]
# omega_{i+1} = 1/(omega_i-a_i)
# omega_i = (omega_0-q_i)/p_i (q_i and p_i are integers.)
# q_{i+1} = p_i a_i - q_i
# p_{i+1} = (D-q_{i+1}^2)/p_i
# q_1 = a_0
# p_1 = D-a_0^2
# See http://en.wikipedia.org/wiki/Pell%27s_equation for more details.
# c_{i+1}=h_{i+1}/k_{i+1} may be writen with h_i and k_i...
# Time-stamp: <2009-12-01 16:34:54 takeshi>
# Author: Takeshi NISHIMATSU
##
require 'rational'
(1..1000).each do |d|
if d==((omega_0=Math::sqrt(d)).to_i)**2
puts d # skip if d is a square number
else
a=[(q=omega_0.floor)] # a_0 and q_1
p=d-q**2 # p_1
omega=1/(omega_0-q) # omega_1
i=0
while (c=Rational(1,a[i]); j=i; while (j-=1)>=0 do; c=1/(a[j]+c); end;
(c.denominator)**2-d*(c.numerator)**2!=1)
a[(i+=1)]=omega.floor
q=p*a[i]-q
p=(d-q**2)/p
omega=(omega_0+q)/p
end
printf("%d %d %d\n", d, c.denominator, c.numerator)
end
end
=============================================

Ciao, ciao,
--
Takeshi
love && peace && free_software
http://loto.sourceforge.net/feram/ Fast MD program for perovskite-type ferroelectrics