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
x2 - D*y2 = 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 Pell's equation - Wikipedia 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-q2 # 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-q2)/p
omega=(omega_0+q)/p
end
printf(“%d %d %d\n”, d, c.denominator, c.numerator)
end
end
Ciao, ciao,