new file: Math/pell_method_for_square_roots.sf

Author: trizen

Math/convergents_to_cube_root_of_2.sf

@@ -2,7 +2,10 @@
 
 # Compute the convergents to cube root of 2.
 
-# Algorithm due to Robert Israel, Oct 08 2017.
+# Formula due to Enrico Bombieri and Alfred J. van der Poorten:
+#   https://web.williams.edu/Mathematics/sjmiller/public_html/book/papers/vdp/BombieriPoorten_CFofAlgNumbs.pdf
+
+# Algorithm adapted after Robert Israel, Oct 08 2017.
 
 # See also:
 #   https://oeis.org/A002352 -- Numerators of convergents to cube root of 2.

Math/pell_method_for_square_roots.sf

@@ -0,0 +1,51 @@
+#!/usr/bin/ruby
+
+# Daniel "Trizen" Șuteu
+# Date: 04 February 2019
+# https://github.com/trizen
+
+# Approximate square roots, using Pell's method.
+
+# See also:
+#   https://rosettacode.org/wiki/Pell%27s_equation
+#   https://en.wikipedia.org/wiki/Pell%27s_equation
+#   https://en.wikipedia.org/wiki/Methods_of_computing_square_roots
+
+func pell_square_root(n, eps=1e-22) {
+
+    var x = n.isqrt
+    var y = x
+    var z = 1
+    var r = 2*x
+
+    var (e1, e2) = (1, 0)
+    var (f1, f2) = (0, 1)
+
+    loop {
+
+        y = (r*z - y)
+        z = ((n - y*y) / z)
+        r = round((x + y) / z)      # floor() also works
+
+        (e1, e2) = (e2, r*e2 + e1)
+        (f1, f2) = (f2, r*f2 + f1)
+
+        var A = (e2 + x*f2)
+        var B = f2
+
+        if (abs((A/B)**2 - n) <= eps) {
+            return A/B
+        }
+    }
+}
+
+for n in [61, 109, 181, 277] {
+    var s = pell_square_root(n)
+    say "sqrt(#{'%3d' % n}) =~ #{s.as_rat} =~ #{s}"
+}
+
+__END__
+sqrt( 61) =~ 5380205503727/688864726095 =~ 7.81024967590665439412972327539046655850673601909
+sqrt(109) =~ 5886776254306/563850903159 =~ 10.4403065089105501797577550769986386680615397574
+sqrt(181) =~ 2900300962727/215577672795 =~ 13.45362404707371031716308866094140080774035985437
+sqrt(277) =~ 4157003026204/249770104837 =~ 16.64331697709323806892821623002200061329832127014

Math/sqrt_convergents.sf

@@ -20,15 +20,13 @@ func sqrt_convergents (n, callback, count=10) {
     count.times {
 
         y = (r*z - y)
-        z = ((n - y*y) // z)                  #/
-        r = ((x + y)   // z)                  #/
+        z = ((n - y*y) / z)
+        r = ((x + y)  // z)                  #/
 
         callback(e2 + x*f2, f2)
 
         (f1, f2) = (f2, r*f2 + f1)
         (e1, e2) = (e2, r*e2 + e1)
-
-        y = x if (z == 1)
     }
 }
 

README.md

@@ -252,6 +252,7 @@ A simple collection of Sidef scripts.
     * [Partial sums of sigma function times k^m](./Math/partial_sums_of_sigma_function_times_k^m.sf)
     * [Partial sums sublinear formula](./Math/partial_sums_sublinear_formula.sf)
     * [Pell factorization](./Math/pell_factorization.sf)
+    * [Pell method for square roots](./Math/pell_method_for_square_roots.sf)
     * [Perfect squares in catalan's triangle](./Math/perfect_squares_in_catalan_s_triangle.sf)
     * [Permutations iter](./Math/permutations_iter.sf)
     * [Permutations rec](./Math/permutations_rec.sf)