perl-scripts/Math/square_root_convergents.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 17 April 2018
# https://github.com/trizen

# Find the convergents of a square root for a non-square positive integer.

# See also:
#    https://en.wikipedia.org/wiki/Pell%27s_equation#Solutions
#    https://en.wikipedia.org/wiki/Continued_fraction#Infinite_continued_fractions
#    https://www.wolframalpha.com/input/?i=Convergents%5BSqrt%5B61%5D%5D

use 5.020;
use strict;
use warnings;

use experimental qw(signatures);
use Math::AnyNum qw(:overload isqrt idiv);

sub sqrt_convergents ($n, $callback, $count = 10) {
    my $x = isqrt($n);
    my $y = $x;
    my $z = 1;

    my $r = $x + $x;

    my ($e1, $e2) = (1, 0);
    my ($f1, $f2) = (0, 1);

    for (1 .. $count) {
        $y = $r * $z - $y;
        $z = idiv($n - $y * $y, $z);
        $r = idiv($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);
    }
}

sqrt_convergents(61, sub ($n, $d) {
        printf("%20s / %-20s =~ %s\n", $n, $d, ($n / $d)->as_dec);
}, 20)

__END__
           7 / 1                    =~ 7
           8 / 1                    =~ 8
          39 / 5                    =~ 7.8
         125 / 16                   =~ 7.8125
         164 / 21                   =~ 7.80952380952380952380952380952380952380952380952
         453 / 58                   =~ 7.81034482758620689655172413793103448275862068966
        1070 / 137                  =~ 7.8102189781021897810218978102189781021897810219
        1523 / 195                  =~ 7.81025641025641025641025641025641025641025641026
        5639 / 722                  =~ 7.81024930747922437673130193905817174515235457064
       24079 / 3083                 =~ 7.81024975673045734674018812844631852092118066818
       29718 / 3805                 =~ 7.81024967148488830486202365308804204993429697766
      440131 / 56353                =~ 7.81024967614856351924476070484268805564921122212
      469849 / 60158                =~ 7.81024967585358555803051963163669004953622128395
     2319527 / 296985               =~ 7.81024967590955772177045978753135680253211441655
     7428430 / 951113               =~ 7.81024967590601747636716142035699228167420695543
     9747957 / 1248098              =~ 7.81024967590685987799035011673762797472634360443
    26924344 / 3447309              =~ 7.81024967590662745927330564216900776808809422074
    63596645 / 8142716              =~ 7.81024967590666308391450714970287555159728031777
    90520989 / 11590025             =~ 7.81024967590665248780740334900054141384509524354
   335159612 / 42912791             =~ 7.81024967590665449842216042298437312082544339752