solve_pell_equation.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 02 February 2019
# https://github.com/trizen

# Find the smallest solution in positive integers to Pell's equation: x^2 - d*y^2 = 1, where `d` is known.

# See also:
#   https://rosettacode.org/wiki/Pell%27s_equation
#   https://en.wikipedia.org/wiki/Pell%27s_equation

func solve_pell(n, w = 1) {

    return nil if n.is_square

    var x = n.isqrt
    var y = x
    var z = 1
    var r = 2*x

    var (e1, e2) = (1, 0)
    var (f1, f2) = (0, 1)

    n.times {

        y = (r*z - y)
        z = floor((n - y*y) / z)
        r = floor((x + y) / z)

        var A = (e2 + x*f2)
        var B = f2

        if ((z == w.abs) && (A**2 - n*B**2 == w)) {
            return (A, B)
        }

        (e1, e2) = (e2, r*e2 + e1)
        (f1, f2) = (f2, r*f2 + f1)
    }

    return nil
}

for d in [-3, -1, 1, 3] {
    for n in [61, 109, 181, 277] {
        var (x, y) = solve_pell(n, d)
        printf("x^2 - %3d*y^2 = %2s for x = %-21s and y = %s\n", n, x**2 - n*y**2, x, y)
    }
}

__END__
x^2 -  61*y^2 = -3 for x = 5639                  and y = 722
x^2 - 109*y^2 = -3 for x = 1399                  and y = 134
x^2 - 181*y^2 = -3 for x = 11262809              and y = 837158
x^2 - 277*y^2 = -3 for x = 233                   and y = 14
x^2 -  61*y^2 = -1 for x = 29718                 and y = 3805
x^2 - 109*y^2 = -1 for x = 8890182               and y = 851525
x^2 - 181*y^2 = -1 for x = 1111225770            and y = 82596761
x^2 - 277*y^2 = -1 for x = 8920484118            and y = 535979945
x^2 -  61*y^2 =  1 for x = 1766319049            and y = 226153980
x^2 - 109*y^2 =  1 for x = 158070671986249       and y = 15140424455100
x^2 - 181*y^2 =  1 for x = 2469645423824185801   and y = 183567298683461940
x^2 - 277*y^2 =  1 for x = 159150073798980475849 and y = 9562401173878027020
x^2 -  61*y^2 =  3 for x = 8                     and y = 1
x^2 - 109*y^2 =  3 for x = 9532                  and y = 913
x^2 - 181*y^2 =  3 for x = 148                   and y = 11
x^2 - 277*y^2 =  3 for x = 57427216              and y = 3450467