Pell's equation solver in Ruby
in [Ruby]
|
From: Takeshi NISHIMATSU on 1 Dec 2009 03:34 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
|
Pages: 1 Prev: one more regular expression, please Next: How to use "eval" to create method parameters? |