# Pell's equation solver in Ruby **URL:** \ **Category:** Ruby\ **Created:** [December 1, 2009, 9:38am UTC](https://www.ruby-forum.com/t/pells-equation-solver-in-ruby/179654 "2009-12-01T09:38:19Z")\ **Posts on this page:** 1\ **Page:** 1
**Author:** ![Takeshi\_NISHIMATSU](https://www.ruby-forum.com/letter_avatar_proxy/v4/letter/t/ce73a5/32.png) [@Takeshi\_NISHIMATSU](https://www.ruby-forum.com/u/Takeshi_NISHIMATSU)\ **Post date:** [December 1, 2009, 9:38am UTC](https://www.ruby-forum.com/t/pells-equation-solver-in-ruby/179654/1 "2009-12-01T09:38:19Z")
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 [Pell's equation - Wikipedia](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,