Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2)

Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2)
You are encouraged to solve this task according to the task description, using any language you may know.

This task performs the basic mathematical functions on 2 continued fractions. This requires the full version of matrix NG:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$

I may perform perform the following operations:

Input the next term of continued fraction N1
Input the next term of continued fraction N2
Output a term of the continued fraction resulting from the operation.

I output a term if the integer parts of ${\displaystyle {\frac {a}{b}}}$ and ${\displaystyle {\frac {a_{1}}{b_{1}}}}$ and ${\displaystyle {\frac {a_{2}}{b_{2}}}}$ and ${\displaystyle {\frac {a_{12}}{b_{12}}}}$ are equal. Otherwise I input a term from continued fraction N1 or continued fraction N2. If I need a term from N but N has no more terms I inject ${\displaystyle \infty }$.

When I input a term t from continued fraction N1 I change my internal state:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$ is transposed thus ${\displaystyle {\begin{bmatrix}a_{2}+a_{12}*t&a+a_{1}*t&a_{12}&a_{1}\\b_{2}+b_{12}*t&b+b_{1}*t&b_{12}&b_{1}\end{bmatrix}}}$

When I need a term from exhausted continued fraction N1 I change my internal state:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$ is transposed thus ${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{12}&a_{1}\\b_{12}&b_{1}&b_{12}&b_{1}\end{bmatrix}}}$

When I input a term t from continued fraction N2 I change my internal state:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$ is transposed thus ${\displaystyle {\begin{bmatrix}a_{1}+a_{12}*t&a_{12}&a+a_{2}*t&a_{2}\\b_{1}+b_{12}*t&b_{12}&b+b_{2}*t&b_{2}\end{bmatrix}}}$

When I need a term from exhausted continued fraction N2 I change my internal state:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$ is transposed thus ${\displaystyle {\begin{bmatrix}a_{12}&a_{12}&a_{2}&a_{2}\\b_{12}&b_{12}&b_{2}&b_{2}\end{bmatrix}}}$

When I output a term t I change my internal state:

${\displaystyle {\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}}$ is transposed thus ${\displaystyle {\begin{bmatrix}b_{12}&b_{1}&b_{2}&b\\a_{12}-b_{12}*t&a_{1}-b_{1}*t&a_{2}-b_{2}*t&a-b*t\end{bmatrix}}}$

When I need to choose to input from N1 or N2 I act:

if b and b2 are zero I choose N1
if b is zero I choose N2
if b2 is zero I choose N2
if abs(${\displaystyle {\frac {a_{1}}{b_{1}}}-{\frac {a}{b}})}$ is greater than abs(${\displaystyle {\frac {a_{2}}{b_{2}}}-{\frac {a}{b}})}$ I choose N1
otherwise I choose N2

When performing arithmetic operation on two potentially infinite continued fractions it is possible to generate a rational number. eg ${\displaystyle {\sqrt {2}}}$ * ${\displaystyle {\sqrt {2}}}$ should produce 2. This will require either that I determine that my internal state is approaching infinity, or limiting the number of terms I am willing to input without producing any output.

Translation of: Python
Works with: GCC version 12.2.1
pragma ada_2022;                -- When big_integers were introduced.

package CONTINUED_FRACTIONS is

type memoization_storage is array (natural range <>) of big_integer;
type memoization_access is access memoization_storage;

type continued_fraction_record is abstract tagged
record
terminated : boolean := false;   -- Are there no more terms?
memo_count : natural := 0;       -- How many terms are memoized?
memo       : memoization_access  -- Memoized terms.
:= new memoization_storage (0 .. 31);
end record;

procedure generate_term (cf : in out continued_fraction_record;
output_exists : out boolean;
output        : out big_integer) is abstract;

type continued_fraction is access all
continued_fraction_record'class; -- The 'class notation is important.

function term_exists (cf : in continued_fraction;
i  : in natural)
return boolean;

function get_term (cf : in continued_fraction;
i  : in natural)
return big_integer
with pre => i < cf.memo_count;

function cf2string (cf        : in continued_fraction;
max_terms : in positive := 20)
return unbounded_string;

end CONTINUED_FRACTIONS;

package body CONTINUED_FRACTIONS is

function term_exists (cf : in continued_fraction;
i  : in natural)
return boolean is
procedure resize_if_necessary is
memo1 : memoization_access;
begin
if cf.memo'length <= i then
memo1 := new memoization_storage(0 .. 2 * (i + 1));
for i in 0 .. cf.memo_count - 1 loop
memo1(i) := cf.memo(i);
end loop;
cf.memo := memo1;
end if;
end;
exists : boolean;
term   : big_integer;
begin
if i < cf.memo_count then
exists := true;
elsif cf.terminated then
exists := false;
else
resize_if_necessary;
while cf.memo_count <= i and not cf.terminated loop
generate_term (cf.all, exists, term);
if exists then
cf.memo(cf.memo_count) := term;
cf.memo_count := cf.memo_count + 1;
else
cf.terminated := true;
end if;
end loop;
exists := term_exists (cf, i);
end if;
return exists;
end;

function get_term (cf : in continued_fraction;
i  : in natural)
return big_integer is
begin
return cf.memo(i);
end;

function cf2string (cf        : in continued_fraction;
max_terms : in positive := 20)
return unbounded_string is
s    : unbounded_string := null_unbounded_string;
done : boolean;
i    : natural;
term : big_integer;
begin
s := s & "[";
i := 0;
done := false;
while not done loop
if not term_exists (cf, i) then
s := s & "]";
done := true;
elsif i = max_terms then
s := s & ",...]";
done := true;
else
if i = 1 then
s := s & ";";
elsif i /= 0 then
s := s & ",";
end if;
term := get_term (cf, i);
s := s & trim (term'image, left);
i := i + 1;
end if;
end loop;
return s;
end;

end CONTINUED_FRACTIONS;

package CONSTANT_TERM_CONTINUED_FRACTIONS is

use CONTINUED_FRACTIONS;

type constant_term_continued_fraction_record is
new continued_fraction_record with
record
term : big_integer;
end record;

type constant_term_continued_fraction is access all
constant_term_continued_fraction_record;

function constant_term_cf (term : in big_integer)
return continued_fraction;

overriding
procedure generate_term (cf : in out constant_term_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer);

end CONSTANT_TERM_CONTINUED_FRACTIONS;

package body CONSTANT_TERM_CONTINUED_FRACTIONS is

function constant_term_cf (term : in big_integer)
return continued_fraction is
cf : constant_term_continued_fraction;
begin
cf := new constant_term_continued_fraction_record;
cf.term := term;
return continued_fraction (cf);
end;

overriding
procedure generate_term (cf : in out constant_term_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer) is
begin
output_exists := true;
output := cf.term;
end;

end CONSTANT_TERM_CONTINUED_FRACTIONS;

package INTEGER_CONTINUED_FRACTIONS is

use CONTINUED_FRACTIONS;

type integer_continued_fraction_record is
new continued_fraction_record with
record
term : big_integer;
done : boolean := false;
end record;

type integer_continued_fraction is access all
integer_continued_fraction_record;

function i2cf (term : in big_integer)
return continued_fraction;

overriding
procedure generate_term (cf : in out integer_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer);

end INTEGER_CONTINUED_FRACTIONS;

package body INTEGER_CONTINUED_FRACTIONS is

function i2cf (term : in big_integer)
return continued_fraction is
cf : integer_continued_fraction;
begin
cf := new integer_continued_fraction_record;
cf.term := term;
return continued_fraction (cf);
end;

overriding
procedure generate_term (cf : in out integer_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer) is
begin
output_exists := not (cf.done);
cf.done := true;
if output_exists then
output := cf.term;
end if;
end;

end INTEGER_CONTINUED_FRACTIONS;

package NG8_CONTINUED_FRACTIONS is

use CONTINUED_FRACTIONS;

stopping_processing_threshold : big_integer := 2 ** 512;
infinitization_threshold      : big_integer := 2 ** 64;

type ng8_continued_fraction_record is
new continued_fraction_record with
record
a12, a1, a2, a : big_integer;
b12, b1, b2, b : big_integer;
x, y           : continued_fraction;
ix, iy         : natural;
xoverflow      : boolean;
yoverflow      : boolean;
end record;

type ng8_continued_fraction is access all
ng8_continued_fraction_record;

function apply_ng8 (a12, a1, a2, a : in big_integer;
b12, b1, b2, b : in big_integer;
x, y           : in continued_fraction)
return continued_fraction;

function "+" (x, y : in continued_fraction)
return continued_fraction;
function "+" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "+" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;

-- Keeping the same sign. (Effectively clones x as an
-- ng8_continued_fraction.)
function "+" (x : in continued_fraction)
return continued_fraction;

-- Subtraction.
function "-" (x, y : in continued_fraction)
return continued_fraction;
function "-" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "-" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;

-- Negation.
function "-" (x : in continued_fraction)
return continued_fraction;

-- Multiplication.
function "*" (x, y : in continued_fraction)
return continued_fraction;
function "*" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "*" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;

-- Division.
function "/" (x, y : in continued_fraction)
return continued_fraction;
function "/" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "/" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;

-- A rational number as a continued fraction. The terms are
-- memoized, so this implementation will not be as inefficient as
-- one might suppose.
function r2cf (n, d : in big_integer)
return continued_fraction;

overriding
procedure generate_term (cf : in out ng8_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer);

end NG8_CONTINUED_FRACTIONS;

package body NG8_CONTINUED_FRACTIONS is

use CONTINUED_FRACTIONS;
use CONSTANT_TERM_CONTINUED_FRACTIONS;

-- An arbitrary infinite source of non-zero finite terms.
forever_cf : continued_fraction := constant_term_cf (1234);

function apply_ng8 (a12, a1, a2, a : in big_integer;
b12, b1, b2, b : in big_integer;
x, y           : in continued_fraction)
return continued_fraction is
cf : ng8_continued_fraction;
begin
cf := new ng8_continued_fraction_record;
cf.a12 := a12;
cf.a1  := a1;
cf.a2  := a2;
cf.a   := a;
cf.b12 := b12;
cf.b1  := b1;
cf.b2  := b2;
cf.b   := b;
cf.x   := x;
cf.y   := y;
cf.ix  := 0;
cf.iy  := 0;
cf.xoverflow := false;
cf.yoverflow := false;
return continued_fraction (cf);
end;

function "+" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1, x, y);
end;

function "+" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, y, 0, 0, 0, 1, x, forever_cf);
end;

function "+" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 1, x, 0, 0, 0, 1, forever_cf, y);
end;

function "+" (x : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 1, 0, 0, 0, 0, 1, forever_cf, x);
end;

function "-" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1, x, y);
end;

function "-" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, -y, 0, 0, 0, 1, x, forever_cf);
end;

function "-" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, -1, x, 0, 0, 0, 1, forever_cf, y);
end;

function "-" (x : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, -1, 0, 0, 0, 0, 1, forever_cf, x);
end;

function "*" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1, x, y);
end;

function "*" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, y, 0, 0, 0, 0, 0, 1, x, forever_cf);
end;

function "*" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, x, 0, 0, 0, 0, 1, forever_cf, y);
end;

function "/" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0, x, y);
end;

function "/" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, 0, 0, 0, 0, y, x, forever_cf);
end;

function "/" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 0, x, 0, 0, 1, 0, forever_cf, y);
end;

function r2cf (n, d : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 0, 0, n, 0, 0, 0, d, forever_cf, forever_cf);
end;

procedure possibly_infinitize_output (q             : in big_integer;
output_exists : out boolean;
output        : out big_integer) is
begin
output_exists := abs (q) < abs (infinitization_threshold);
if output_exists then
output := q;
end if;
end;

procedure divide (a, b : in big_integer;
q, r : out big_integer) is
begin
if b /= 0 then
q := a / b;
r := a rem b;
end if;
end;

function too_big (num : big_integer)
return boolean is
begin
return (abs (stopping_processing_threshold) <= abs (num));
end;

function any_too_big (a, b, c, d, e, f, g, h : in big_integer)
return boolean is
begin
return (too_big (a) or else
too_big (b) or else
too_big (c) or else
too_big (d) or else
too_big (e) or else
too_big (f) or else
too_big (g) or else
too_big (h));
end;

overriding
procedure generate_term (cf : in out ng8_continued_fraction_record;
output_exists : out boolean;
output        : out big_integer) is

a12, a1, a2, a        : big_integer;
b12, b1, b2, b        : big_integer;
q12, q1, q2, q        : big_integer;
r12, r1, r2, r        : big_integer;
done                  : boolean;

function all_b_are_zero
return boolean is
begin
return (b12 = 0 and b1 = 0 and b2 = 0 and b = 0);
end;

function all_q_are_equal
return boolean is
begin
return (q = q1 and q = q2 and q = q12);
end;

procedure compare_fractions is
n1, n2, n : big_integer;
begin
-- Rather than compare fractions, we will put the numerators over
-- a common denominator of b*b1*b2, and then compare the new
-- numerators.
n1 := a1 * b2 * b;
n2 := a2 * b1 * b;
n  := a  * b1 * b2;
absorb_y_instead_of_x := (abs (n1 - n) <= abs (n2 - n));
end;

procedure absorb_x_term is
term                           : big_integer;
new_a12, new_a1, new_a2, new_a : big_integer;
new_b12, new_b1, new_b2, new_b : big_integer;
begin
new_a2 := a12;
new_a  := a1;
new_b2 := b12;
new_b  := b1;
if not cf.xoverflow and then term_exists (cf.x, cf.ix) then
term := get_term (cf.x, cf.ix);
new_a12 := a2 + (a12 * term);
new_a1  := a  + (a1  * term);
new_b12 := b2 + (b12 * term);
new_b1  := b  + (b1  * term);
if any_too_big (new_a12, new_a1, new_a2, new_a,
new_b12, new_b1, new_b2, new_b) then
cf.xoverflow := true;
new_a12 := a12;
new_a1  := a1;
new_b12 := b12;
new_b1  := b1;
else
cf.ix := cf.ix + 1;
end if;
else
new_a12 := a12;
new_a1  := a1;
new_b12 := b12;
new_b1  := b1;
end if;
a12 := new_a12;
a1  := new_a1;
a2  := new_a2;
a   := new_a;
b12 := new_b12;
b1  := new_b1;
b2  := new_b2;
b   := new_b;
end;

procedure absorb_y_term is
term                           : big_integer;
new_a12, new_a1, new_a2, new_a : big_integer;
new_b12, new_b1, new_b2, new_b : big_integer;
begin
new_a1 := a12;
new_a  := a2;
new_b1 := b12;
new_b  := b2;
if not cf.yoverflow and then term_exists (cf.y, cf.iy) then
term := get_term (cf.y, cf.iy);
new_a12 := a1 + (a12 * term);
new_a2  := a  + (a2  * term);
new_b12 := b1 + (b12 * term);
new_b2  := b  + (b2  * term);
if any_too_big (new_a12, new_a1, new_a2, new_a,
new_b12, new_b1, new_b2, new_b) then
cf.yoverflow := true;
new_a12 := a12;
new_a2  := a2;
new_b12 := b12;
new_b2  := b2;
else
cf.iy := cf.iy + 1;
end if;
else
new_a12 := a12;
new_a2  := a2;
new_b12 := b12;
new_b2  := b2;
end if;
a12 := new_a12;
a1  := new_a1;
a2  := new_a2;
a   := new_a;
b12 := new_b12;
b1  := new_b1;
b2  := new_b2;
b   := new_b;
end;

procedure absorb_term is
begin
absorb_y_term;
else
absorb_x_term;
end if;
end;

begin
a12 := cf.a12;
a1  := cf.a1;
a2  := cf.a2;
a   := cf.a;
b12 := cf.b12;
b1  := cf.b1;
b2  := cf.b2;
b   := cf.b;

done := false;
while not done loop
if all_b_are_zero then
-- There are no more terms.
output_exists := false;
done := true;
elsif b2 = 0 and b = 0 then
null;
elsif b2 = 0 or b = 0 then
elsif b1 = 0 then
null;
else
divide (a12, b12, q12, r12);
divide (a1, b1, q1, r1);
divide (a2, b2, q2, r2);
divide (a, b, q, r);
if b12 /= 0 and then all_q_are_equal then
-- Output a term.
cf.a12 := b12;
cf.a1  := b1;
cf.a2  := b2;
cf.a   := b;
cf.b12 := r12;
cf.b1  := r1;
cf.b2  := r2;
cf.b   := r;
possibly_infinitize_output (q, output_exists, output);
done := true;
else
compare_fractions;
end if;
end if;
absorb_term;
end loop;
end;

end NG8_CONTINUED_FRACTIONS;

use CONTINUED_FRACTIONS;
use CONSTANT_TERM_CONTINUED_FRACTIONS;
use INTEGER_CONTINUED_FRACTIONS;
use NG8_CONTINUED_FRACTIONS;

procedure show (expression : in string;
cf         : in continued_fraction;
note       : in string := "") is
expr     : string := 19 * ' ';
contfrac : string := 48 * ' ';
begin
move (source => expression,
target => expr,
justify => right);
put (expr);
put (" =>  ");
if note = "" then
put_line (to_string (cf2string (cf)));
else
move (source => to_string (cf2string (cf)),
target => contfrac,
justify => left);
put (contfrac);
put_line (note);
end if;
end;

golden_ratio : continued_fraction := constant_term_cf (1);
silver_ratio : continued_fraction := constant_term_cf (2);
one          : continued_fraction := i2cf (1);
two          : continued_fraction := i2cf (2);
three        : continued_fraction := i2cf (3);
four         : continued_fraction := i2cf (4);
sqrt2        : continued_fraction := silver_ratio - 1;

begin

show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "(1 + sqrt(2))");
show ("sqrt(2)", sqrt2, "silver ratio minus 1");
show ("13/11", r2cf (13, 11));
show ("22/7", r2cf (22, 7), "approximately pi");
show ("1", one);
show ("2", two);
show ("3", three);
show ("4", four);
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (0, 1, 0, 0, 0, 0, 2, 0, silver_ratio, sqrt2),
"method 1");
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (1, 0, 0, 1, 0, 0, 0, 8, silver_ratio, silver_ratio),
"method 2");
show ("(1 + 1/sqrt(2))/2",
(one / 2) * (one + (1 / sqrt2)),
"method 3");
show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2);
show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2);
show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2);
show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2);

-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end:

Output:
$gnatmake -q -g bivariate_continued_fraction_task.adb && ./bivariate_continued_fraction_task golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2 silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] (1 + sqrt(2)) sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] silver ratio minus 1 13/11 => [1;5,2] 22/7 => [3;7] approximately pi 1 => [1] 2 => [2] 3 => [3] 4 => [4] (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3 sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] sqrt(2) - sqrt(2) => [0] sqrt(2) * sqrt(2) => [2] sqrt(2) / sqrt(2) => [1]  ATS Translation of: Python Using 128-bit integers (Margin note: This program is a bug-fix of an ATS program on which I based Python code. It does not really matter, however, which program came first. So I am calling this a translation from Python.) The following program uses GNU C extensions: 128-bit integers, and integer operations that detect overflow. I add the 128-bit integers as a new g0int(int128knd) type. The overflow-detecting operations I call +!, -!, and *!. There are no multiple-precision integers or rationals. Rather than detect numbers "getting too large", I catch integer overflows. I use the most negative 128-bit value to represent "no finite term". I have not included any code to weed out large terms that show up where "infinities" belong. (Margin note: Sometimes infinite sequences get truncated due to overflow, though I have not seen this happen very near the beginning of a continued fraction. You can give a "maximum number of terms to print" argument to this program, to see the phenomenon for yourself. It happens relatively quickly with 128-bit integers.) (*------------------------------------------------------------------*) #include "share/atspre_staload.hats" staload UN = "prelude/SATS/unsafe.sats" %{# #include <stdint.h> #include <limits.h> #include <float.h> #include <math.h> %} #define NIL list_nil () #define :: list_cons exception gint_overflow of () (*------------------------------------------------------------------*) extern fn {tk : tkind} g0int_neginf : () -<> g0int tk extern fn {tk : tkind} g0int_add_overflow_exn : (g0int tk, g0int tk) -< !exn > g0int tk extern fn {tk : tkind} g0int_sub_overflow_exn : (g0int tk, g0int tk) -< !exn > g0int tk extern fn {tk : tkind} g0int_mul_overflow_exn : (g0int tk, g0int tk) -< !exn > g0int tk infixl ( + ) +! infixl ( - ) -! infixl ( * ) *! overload neginf with g0int_neginf overload +! with g0int_add_overflow_exn overload -! with g0int_sub_overflow_exn overload *! with g0int_mul_overflow_exn (*------------------------------------------------------------------*) (* 128-bit integers. *) %{# /* The most negative int128 will be treated as "neginf" or "negative infinity". For our purposes the sign will not matter, though. */ #define neginf_int128() (((__int128) 1) << 127) #define neg_c(x) (-(x)) #define add_c(x, y) ((x) + (y)) #define sub_c(x, y) ((x) - (y)) #define mul_c(x, y) ((x) * (y)) #define div_c(x, y) ((x) / (y)) #define mod_c(x, y) ((x) % (y)) #define eq_c(x, y) (((x) == (y)) ? atsbool_true : atsbool_false) #define neq_c(x, y) (((x) != (y)) ? atsbool_true : atsbool_false) #define lt_c(x, y) (((x) < (y)) ? atsbool_true : atsbool_false) #define lte_c(x, y) (((x) <= (y)) ? atsbool_true : atsbool_false) #define gt_c(x, y) (((x) > (y)) ? atsbool_true : atsbool_false) #define gte_c(x, y) (((x) >= (y)) ? atsbool_true : atsbool_false) /* GNU extensions for detection of integer overflow. */ #define add_overflow(x, y, pz) \ (__builtin_add_overflow ((x), (y), (pz)) ? \ atsbool_true : atsbool_false) #define sub_overflow(x, y, pz) \ (__builtin_sub_overflow ((x), (y), (pz)) ? \ atsbool_true : atsbool_false) #define mul_overflow(x, y, pz) \ (__builtin_mul_overflow ((x), (y), (pz)) ? \ atsbool_true : atsbool_false) %} tkindef int128_kind = "__int128" (* A GNU extension. *) stadef int128knd = int128_kind typedef int128_0 = g0int int128knd typedef int128_1 (i : int) = g1int (int128knd, i) stadef int128 = int128_1 // 2nd-select stadef int128 = int128_0 // 1st-select stadef Int128 = [i : int] int128_1 i extern fn g0int_neginf_int128 : () -<> int128 = "mac#neginf_int128" extern fn g0int_neg_int128 : int128 -<> int128 = "mac#neg_c" extern fn g0int_add_int128 : (int128, int128) -<> int128 = "mac#add_c" extern fn g0int_sub_int128 : (int128, int128) -<> int128 = "mac#sub_c" extern fn g0int_mul_int128 : (int128, int128) -<> int128 = "mac#mul_c" extern fn g0int_div_int128 : (int128, int128) -<> int128 = "mac#div_c" extern fn g0int_mod_int128 : (int128, int128) -<> int128 = "mac#mod_c" extern fn g0int_eq_int128 : (int128, int128) -<> bool = "mac#eq_c" extern fn g0int_neq_int128 : (int128, int128) -<> bool = "mac#neq_c" extern fn g0int_lt_int128 : (int128, int128) -<> bool = "mac#lt_c" extern fn g0int_lte_int128 : (int128, int128) -<> bool = "mac#lte_c" extern fn g0int_gt_int128 : (int128, int128) -<> bool = "mac#gt_c" extern fn g0int_gte_int128 : (int128, int128) -<> bool = "mac#gte_c" implement g0int_neginf<int128knd> () = g0int_neginf_int128 () implement g0int_neg<int128knd> x = g0int_neg_int128 x implement g0int_add<int128knd> (x, y) = g0int_add_int128 (x, y) implement g0int_sub<int128knd> (x, y) = g0int_sub_int128 (x, y) implement g0int_mul<int128knd> (x, y) = g0int_mul_int128 (x, y) implement g0int_div<int128knd> (x, y) = g0int_div_int128 (x, y) implement g0int_mod<int128knd> (x, y) = g0int_mod_int128 (x, y) implement g0int_eq<int128knd> (x, y) = g0int_eq_int128 (x, y) implement g0int_neq<int128knd> (x, y) = g0int_neq_int128 (x, y) implement g0int_lt<int128knd> (x, y) = g0int_lt_int128 (x, y) implement g0int_lte<int128knd> (x, y) = g0int_lte_int128 (x, y) implement g0int_gt<int128knd> (x, y) = g0int_gt_int128 (x, y) implement g0int_gte<int128knd> (x, y) = g0int_gte_int128 (x, y) implement g0int2int<intknd,int128knd> i =$UN.cast i
implement g0int2int<int128knd,intknd> i = $UN.cast i implement g0int2float<int128knd,ldblknd> i =$UN.cast i

implement g0int_iseqz<int128knd> x = (x = g0i2i 0)
implement g0int_isneqz<int128knd> x = (x <> g0i2i 0)
implement g0int_isltz<int128knd> x = (x < g0i2i 0)
implement g0int_isltez<int128knd> x = (x <= g0i2i 0)
implement g0int_isgtz<int128knd> x = (x > g0i2i 0)
implement g0int_isgtez<int128knd> x = (x >= g0i2i 0)

implement g0int_abs<int128knd> x = (if isltz x then ~x else x)

local

extern fn
(int128, int128, &int128? >> int128) -< !wrt > bool

extern fn
sub_overflow_int128 :
(int128, int128, &int128? >> int128) -< !wrt > bool
= "mac#sub_overflow"

extern fn
mul_overflow_int128 :
(int128, int128, &int128? >> int128) -< !wrt > bool
= "mac#mul_overflow"

in (* local *)

fn
g0int_add_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt add_overflow_int128 (x, y, z) in if ~overflow then z else$raise gint_overflow ()
end

fn
g0int_sub_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt sub_overflow_int128 (x, y, z) in if ~overflow then z else$raise gint_overflow ()
end

fn
g0int_mul_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt mul_overflow_int128 (x, y, z) in if ~overflow then z else$raise gint_overflow ()
end

end (* local *)

implement

implement
g0int_sub_overflow_exn<int128knd> (x, y) =
g0int_sub_overflow_exn_int128 (x, y)

implement
g0int_mul_overflow_exn<int128knd> (x, y) =
g0int_mul_overflow_exn_int128 (x, y)

(*------------------------------------------------------------------*)

(* We will truncate quotients towards zero. *)
infixl ( / ) div divrem
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod
macdef divrem (a, b) =
let
val x = ,(a)
and y = ,(b)
in
(* Optimizing C compilers will compute both quotient and remainder
at the same time. *)
@(x \g0int_div y, x \g0int_mod y)
end

(*------------------------------------------------------------------*)
(* Continued fractions.

cf_generator tk -- A closure that produces terms of type g0int tk,
sequentially.

cf tk           -- A structure from which one can get the ith
term of a continued fraction. It gets terms
from a cf_generator tk.                       *)

typedef integer = int128

typedef cf_generator = () -<cloref1> integer

local

typedef _cf (terminated : bool,
m          : int,
n          : int) =
[m <= n]
@{
terminated = bool terminated, (* No more terms? *)
m = size_t m,         (* The number of terms computed so far. *)
n = size_t n,         (* The size of memo storage.*)
memo = arrayref (integer, n), (* Memoized terms. *)
gen = cf_generator            (* A thunk to generate terms. *)
}
typedef _cf (m : int) =
[terminated : bool]
[n : int | m <= n]
_cf (terminated, m, n)
typedef _cf =
[m : int]
_cf m

fn
_cf_get_more_terms
{terminated : bool}
{m          : int}
{n          : int}
{needed     : int | m <= needed; needed <= n}
(cf         : _cf (terminated, m, n),
needed     : size_t needed)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)

macdef memo = cf.memo

fun
loop {i : int | m <= i; i <= needed}
.<needed - i>.
(i : size_t i)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
if i = needed then
@{
terminated = false,
m = needed,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
else
let
val term = (cf.gen) ()
in
if term <> neginf<integerknd> () then
begin
memo[i] := term;
loop (succ i)
end
else
@{
terminated = true,
m = i,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
end
in
loop (cf.m)
end

fn
_cf_update {terminated : bool}
{m          : int}
{n          : int | m <= n}
{needed     : int}
(cf         : _cf (terminated, m, n),
needed     : size_t needed)
: [terminated1 : bool]
[m1 : int | m <= m1]
[n1 : int | m1 <= n1]
_cf (terminated1, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)
macdef memo = cf.memo
macdef gen = cf.gen
in
if (cf.terminated) then
cf
else if needed <= (cf.m) then
cf
else if needed <= (cf.n) then
_cf_get_more_terms (cf, needed)
else
let                     (* Provides twice the room needed. *)
val n1 = needed + needed
val memo1 = arrayref_make_elt (n1, g0i2i 0)
val () =
let
var i : [i : nat] size_t i
in
for (i := i2sz 0; i < (cf.m); i := succ i)
memo1[i] := memo[i]
end
val cf1 =
@{
terminated = false,
m = (cf.m),
n = n1,
memo = memo1,
gen = (cf.gen)
}
in
_cf_get_more_terms (cf1, needed)
end
end

in (* local *)

typedef cf = ref _cf

extern fn cf_make : cf_generator -> cf

extern fn cf_get_at_size : {i : int} (cf, size_t i) -> integer
extern fn cf_get_at_int : {i : nat} (cf, int i) -> integer

extern fn cf2generator : cf -> cf_generator

(* The precedence of the overloads has to exceed that of ref[] *)
overload cf_get_at with cf_get_at_size of 1
overload cf_get_at with cf_get_at_int of 1
overload [] with cf_get_at of 1

implement
cf_make gen =
let
#ifndef CF_START_SIZE #then
#define CF_START_SIZE 32
#endif
in
ref
@{
terminated = false,
m = i2sz 0,
n = i2sz CF_START_SIZE,
memo = arrayref_make_elt (i2sz CF_START_SIZE, g0i2i 0),
gen = gen
}
end

implement
cf_get_at_size (cf, i) =
let
prval () = lemma_g1uint_param i
val cf1 = _cf_update (!cf, succ i)
in
!cf := cf1;
if i < (cf1.m) then
arrayref_get_at<integer> (cf1.memo, i)
else
neginf<integerknd> ()
end

implement cf_get_at_int (cf, i) = cf_get_at_size (cf, g1i2u i)

implement
cf2generator cf =
let
val i : ref Size_t = ref (i2sz 0)
in
lam () =>
let
val term = cf[!i]
in
!i := succ !i;
term
end
end

end (* local *)

(*------------------------------------------------------------------*)
(* Make a string from an int128. *)

fn
int128_2string (i : int128) : string =
let
fun
loop (i : int128, accum : List0 charNZ) : List0 charNZ =
if iseqz i then
accum
else
let
val @(i, digit) = i divrem (g0i2i 10)
val digit = (g0i2i digit) : int
val digit = g1ofg0 (int2char0 (digit + (char2int0 '0')))
val () = assertloc (digit <> '\0')
in
loop (i, digit :: accum)
end

val minus_sign : Char = '-'
val () = assertloc (minus_sign <> '\0')
in
if iseqz i then
"0"
else if i = neginf<integerknd> () then
"neginf"
else if isltz i then
strnptr2string (string_implode (minus_sign :: (loop (~i, NIL))))
else
strnptr2string (string_implode (loop (i, NIL)))
end

implement tostring_val<integer> i = $effmask_all int128_2string i implement fprint_val<integer> (f, i) = fprint! (f, tostring_val<integer> i) fn fprint_integer (f : FILEref, i : integer) = fprint_val<integer> (f, i) fn print_integer (i : integer) = fprint_integer (stdout_ref, i) fn prerr_integer (i : integer) = fprint_integer (stderr_ref, i) overload fprint with fprint_integer overload print with print_integer overload prerr with prerr_integer (*------------------------------------------------------------------*) (* Converting a continued fraction to a string. *) extern fn cf2string_with_default_max_terms : cf -> string extern fn cf2string_given_max_terms : (cf, Size_t) -> string overload cf2string with cf2string_with_default_max_terms overload cf2string with cf2string_given_max_terms val cf2string_default_max_terms : ref Size_t = ref (i2sz 20) implement cf2string_with_default_max_terms cf = cf2string_given_max_terms (cf, !cf2string_default_max_terms) implement cf2string_given_max_terms (cf, max_terms) = let fun loop (i : Size_t, accum : List0 string) : List0 string = let val term = cf[i] in if i = max_terms then begin if term = neginf<integerknd> () then "]" :: accum else ",...]" :: accum end else if term = neginf<integerknd> () then "]" :: accum else let val separator = (if i = i2sz 0 then "" else if i = i2sz 1 then ";" else ",") and term_str = tostring_val<integer> term in loop (succ i, term_str :: separator :: accum) end end val string_lst = loop (i2sz 0, "[" :: NIL) in strptr2string (stringlst_concat (list_vt2t (reverse string_lst))) end (*------------------------------------------------------------------*) (* The continued fraction for a rational number or integer. *) typedef ratnum = @(integer, integer) fn r2cf (ratnum : ratnum) : cf = cf_make let val ratnum_ref : ref ratnum = ref ratnum in lam () =<cloref1> let val @(n, d) = !ratnum_ref in if iseqz d then neginf<integerknd> () else let val @(q, r) = n divrem d in !ratnum_ref := @(d, r); q end end end fn i2cf (intnum : integer) : cf = r2cf @(intnum, g0i2i 1) (*------------------------------------------------------------------*) (* Application of a homographic function to a continued fraction. *) typedef ng4 = @(integer, integer, integer, integer) fn apply_ng4 (ng4 : ng4, x : cf) : cf = cf_make let val state : ref @(ng4, Size_t) = ref @(ng4, i2sz 0) in lam () =<cloref1> let fnx loop (a1 : integer, a : integer, b1 : integer, b : integer, i : Size_t) : integer = if (iseqz b1) * (iseqz b) then neginf<integerknd> () else if (iseqz b1) + (iseqz b) then absorb_term (a1, a, b1, b, i) else let val @(q, r) = a divrem b and @(q1, r1) = a1 divrem b1 in if q1 <> q then absorb_term (a1, a, b1, b, i) else begin !state := @(@(b1, b, r1, r), i); q end end and absorb_term (a1 : integer, a : integer, b1 : integer, b : integer, i : Size_t) : integer = let val term = x[i] in if term <> neginf<integerknd> () then loop (a + (a1 * term), a1, b + (b1 * term), b1, succ i) else loop (a1, a1, b1, b1, succ i) end val @(@(a1, a, b1, b), i) = !state in loop (a1, a, b1, b, i) end end (*------------------------------------------------------------------*) (* Some basic operations involving only one continued fraction. *) fn cf_neg (x : cf) : cf = apply_ng4 (@(g0i2i ~1, g0i2i 0, g0i2i 0, g0i2i 1), x) fn cf_add_ratnum (x : cf, y : ratnum) : cf = let val @(n, d) = y in apply_ng4 (@(d, n, g0i2i 0, d), x) end fn ratnum_add_cf (x : ratnum, y : cf) : cf = cf_add_ratnum (y, x) fn cf_add_integer (x : cf, y : integer) : cf = cf_add_ratnum (x, @(y, g0i2i 1)) fn integer_add_cf (x : integer, y : cf) : cf = cf_add_ratnum (y, @(x, g0i2i 1)) fn cf_sub_ratnum (x : cf, y : ratnum) : cf = cf_add_ratnum (x, @(~(y.0), (y.1))) fn ratnum_sub_cf (x : ratnum, y : cf) : cf = let val @(n, d) = x in apply_ng4 (@(~d, n, g0i2i 0, d), y) end fn cf_sub_integer (x : cf, y : integer) : cf = cf_add_ratnum (x, @(~y, g0i2i 1)) fn integer_sub_cf (x : integer, y : cf) : cf = ratnum_sub_cf (@(x, g0i2i 1), y) fn cf_mul_ratnum (x : cf, y : ratnum) : cf = let val @(n, d) = y in apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), x) end fn ratnum_mul_cf (x : ratnum, y : cf) : cf = cf_mul_ratnum (y, x) fn cf_mul_integer (x : cf, y : integer) : cf = cf_mul_ratnum (x, @(y, g0i2i 1)) fn integer_mul_cf (x : integer, y : cf) : cf = cf_mul_ratnum (y, @(x, g0i2i 1)) fn cf_div_ratnum (x : cf, y : ratnum) : cf = cf_mul_ratnum (x, @(y.1, y.0)) fn ratnum_div_cf (x : ratnum, y : cf) : cf = let val @(n, d) = x in apply_ng4 (@(g0i2i 0, n, d, g0i2i 0), y) end fn cf_div_integer (x : cf, y : integer) : cf = cf_mul_ratnum (x, @(g0i2i 1, y)) fn integer_div_cf (x : integer, y : cf) : cf = ratnum_div_cf (@(x, g0i2i 1), y) overload ~ with cf_neg overload + with cf_add_ratnum overload + with ratnum_add_cf overload + with cf_add_integer overload + with integer_add_cf overload - with cf_sub_ratnum overload - with ratnum_sub_cf overload - with cf_sub_integer overload - with integer_sub_cf overload * with cf_mul_ratnum overload * with ratnum_mul_cf overload * with cf_mul_integer overload * with integer_mul_cf overload / with cf_div_ratnum overload / with ratnum_div_cf overload / with cf_div_integer overload / with integer_div_cf (*------------------------------------------------------------------*) (* Application of a bihomographic function to a continued fraction. *) typedef ng8 = @(integer, integer, integer, integer, integer, integer, integer, integer) fn apply_ng8 (ng8 : ng8, x : cf, y : cf) : cf = cf_make let val ng : ref ng8 = ref ng8 and xsource : ref cf_generator = ref (cf2generator x) and ysource : ref cf_generator = ref (cf2generator y) fn neginf_source () :<cloref1> integer = neginf<integerknd> () in lam () =<cloref1> let fnx recurs () : integer = let val @(a12, a1, a2, a, b12, b1, b2, b) = !ng val bz = iseqz b and b1z = iseqz b1 and b2z = iseqz b2 and b12z = iseqz b12 in if b12z * b1z * b2z * bz then neginf<integerknd> () else if bz * b2z then absorb_x_term () else if bz + b2z then absorb_y_term () else if b1z then absorb_x_term () else let val @(q, r) = a divrem b and @(q1, r1) = a1 divrem b1 and @(q2, r2) = a2 divrem b2 and @(q12, r12) = (if ~b12z then a12 divrem b12 else @(neginf (), neginf ())) : @(integer, integer) in if (~b12z) * (q = q1) * (q = q2) * (q = q12) then begin (* Output a term. *) !ng := @(b12, b1, b2, b, r12, r1, r2, r); q end else let (* Put numerators over a common denominator, then compare the numerators. *) val n = a *! b1 *! b2 and n1 = a1 *! b *! b2 and n2 = a2 *! b *! b1 in if abs (n1 -! n) > abs (n2 -! n) then absorb_x_term () else absorb_y_term () end end end and absorb_x_term () : integer = let val @(a12, a1, a2, a, b12, b1, b2, b) = !ng val term = (!xsource) () in if term <> neginf<integerknd> () then begin try let val a12_ = a2 +! (a12 *! term) and a1_ = a +! (a1 *! term) and a2_ = a12 and a_ = a1 and b12_ = b2 +! (b12 *! term) and b1_ = b +! (b1 *! term) and b2_ = b12 and b_ = b1 in !ng := @(a12_, a1_, a2_, a_, b12_, b1_, b2_, b_); (* Be aware: this is not a tail recursion! *) recurs () end with | ~ gint_overflow () => begin (* Replace the sources with ones that return no terms. (You have to replace BOTH sources.) *) !ng := @(a12, a1, a12, a1, b12, b1, b12, b1); !xsource := neginf_source; !ysource := neginf_source; recurs () end end else begin !ng := @(a12, a1, a12, a1, b12, b1, b12, b1); recurs () end end and absorb_y_term () : integer = let val @(a12, a1, a2, a, b12, b1, b2, b) = !ng val term = (!ysource) () in if term <> neginf<integerknd> () then begin try let val a12_ = a1 +! (a12 *! term) and a1_ = a12 and a2_ = a +! (a2 *! term) and a_ = a2 and b12_ = b1 +! (b12 *! term) and b1_ = b12 and b2_ = b +! (b2 *! term) and b_ = b2 in !ng := @(a12_, a1_, a2_, a_, b12_, b1_, b2_, b_); (* Be aware: this is not a tail recursion! *) recurs () end with | ~ gint_overflow () => begin (* Replace the sources with ones that return no terms. (You have to replace BOTH sources.) *) !ng := @(a12, a12, a2, a2, b12, b12, b2, b2); !xsource := neginf_source; !ysource := neginf_source; recurs () end end else begin !ng := @(a12, a12, a2, a2, b12, b12, b2, b2); recurs () end end in recurs () end end (*------------------------------------------------------------------*) (* Some basic operations on two continued fractions. *) fn cf_add_cf (x : cf, y : cf) : cf = apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y) fn cf_sub_cf (x : cf, y : cf) : cf = apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i ~1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y) fn cf_mul_cf (x : cf, y : cf) : cf = apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y) fn cf_div_cf (x : cf, y : cf) : cf = apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1, g0i2i 0), x, y) overload + with cf_add_cf overload - with cf_sub_cf overload * with cf_mul_cf overload / with cf_div_cf (*------------------------------------------------------------------*) typedef charptr =$extype"char *"

fn
show_with_note (expression : string,
cf         : cf,
note       : string)
: void =
if note = "" then
ignoret ($extfcall (int, "printf", "%19s => %s\n",$UN.cast{charptr} expression,
$UN.cast{charptr} (cf2string cf))) else ignoret ($extfcall (int, "printf", "%19s =>  %-46s  %s\n",
$UN.cast{charptr} expression,$UN.cast{charptr} (cf2string cf),
$UN.cast{charptr} note)) fn show_without_note (expression : string, cf : cf) : void = show_with_note (expression, cf, "") overload show with show_with_note overload show with show_without_note (*------------------------------------------------------------------*) val golden_ratio = cf_make (lam () =<cloref1> (g0i2i 1) : integer) val silver_ratio = cf_make (lam () =<cloref1> (g0i2i 2) : integer) val sqrt2 = silver_ratio - g0i2i 1 val frac_13_11 = r2cf @(g0i2i 13, g0i2i 11) val frac_22_7 = r2cf @(g0i2i 22, g0i2i 7) val one = i2cf (g0i2i 1) val two = i2cf (g0i2i 2) val three = i2cf (g0i2i 3) val four = i2cf (g0i2i 4) implement main (argc, argv) = begin if 2 <= argc then let val n = g0string2uint (argv_get_at (argv, 1)) : uint in (* Set the maximum number of terms to print. *) !cf2string_default_max_terms := g1u2u (g1ofg0 n) end; show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2"); show ("silver ratio", silver_ratio, "1 + sqrt(2)"); show ("sqrt(2)", sqrt2); show ("13/11", frac_13_11); show ("22/7", frac_22_7); show ("one", one); show ("two", two); show ("three", three); show ("four", four); show ("(1 + 1/sqrt(2))/2", apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 2, g0i2i 0), silver_ratio, sqrt2), "method 1"); show ("(1 + 1/sqrt(2))/2", apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 8), silver_ratio, silver_ratio), "method 2"); show ("(1 + 1/sqrt(2))/2", (g0i2i 1 + (one / sqrt2)) / two, "method 3"); show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2); show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2); show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2); show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2); 0 end (*------------------------------------------------------------------*) Output: $ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW bivariate-continued-fraction-task-memoizing.dats -lgc && ./a.out
golden ratio =>  [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]   (1 + sqrt(5))/2
silver ratio =>  [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]   1 + sqrt(2)
sqrt(2) =>  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 =>  [1;5,2]
22/7 =>  [3;7]
one =>  [1]
two =>  [2]
three =>  [3]
four =>  [4]
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 1
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 2
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 3
sqrt(2) + sqrt(2) =>  [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) =>  [0]
sqrt(2) * sqrt(2) =>  [2;7530688524100]
sqrt(2) / sqrt(2) =>  [1]


Using multiple precision numbers

Translation of: Standard ML
Library: ats2-xprelude
Library: GMP

For this program you need ats2-xprelude.

The program closely follows the Standard ML code, so one can compare the two languages with each other. They have similar syntaxes, but are very different. Notice, for instance, that ATS has overloads, whereas SML does not. (SML has signatures with respective namespaces.) In ATS, a function is not a closure unless you explicitly make it one, whereas in SML no special notation is needed. And so on.

ATS is translated to C, and its functions (except closures) are essentially just C functions. One can write Arduino code and Linux kernel modules in ATS, because ATS is, in some sense, an elaborate way to write C. Nevertheless, there is enough similarity between ATS and Standard ML to easily translate the SML code for this Rosetta Code task to ATS.

I have broken the program into three files, to demonstrate what an ATS program might look like, if it were broken into separately compiled parts.

The first file is an "interface" specification for a continued_fraction data type. The file is called continued_fraction.sats:

(* "Static" file. (Exported declarations.) *)

(* To set up a predictable name-mangling scheme:  *)
#define ATS_PACKNAME "rosetta-code.continued_fraction"

(* Load declarations from ats2-xprelude: *)
#include "xprelude/HATS/xprelude_sats.hats"

(* A term_generator thunk generates terms, which a continued_fraction
data structure memoizes. The internals of continued_fraction are
not exposed here. It is an abstract type, the size of (but not the
same type as) a pointer. SIDE NOTE: In ATS2, we get a conventional
function, rather than a closure, unless we say explicitly that we
want a closure; "cloref1" means a particular kind of closure one
often uses when linking the program with Boehm GC. *)
typedef term_generator = () -<cloref1> Option exrat
abstype continued_fraction = ptr

(* Create a continued fraction. *)
fn continued_fraction_make : term_generator -> continued_fraction

(* Does the indexed term exist? *)
fn term_exists : (continued_fraction, intGte 0) -> bool

(* Retrieve the indexed term. Raise an exception if there is no such
term. The precedence of the overload must exceed that of an
overload that is in the prelude. (To see what I mean, try removing
the "of 1".) *)
val get_term_exn : (continued_fraction, intGte 0) -> exrat
overload [] with get_term_exn of 1

(* Use a continued_fraction as a term_generator thunk. *)
fn continued_fraction_to_term_generator :
continued_fraction -> term_generator

(* Get a human-readable string. *)
val default_max_terms : ref (intGte 1)
fn continued_fraction_to_string_given_max_terms :
(continued_fraction, intGte 1) -> string
fn continued_fraction_to_string_default_max_terms :
continued_fraction -> string
continued_fraction_to_string_given_max_terms
continued_fraction_to_string_default_max_terms

(* Make a continued_fraction for an integer. *)
fn int_to_continued_fraction : int -> continued_fraction

(* Make a continued_fraction for a rational number. *)
fn exrat_to_continued_fraction : exrat -> continued_fraction
fn rational_to_continued_fraction :
(int, [d : int | d != 0] int d) -> continued_fraction

(* Make a continued_fraction with one term repeated forever. *)
fn continued_fraction_make_constant_term : int -> continued_fraction

(* Make a continued fraction via binary arithmetic operations. (I have
not bothered here to implement ng4, although one likely would wish
to have ng4 as well.)  *)
(* The @() denotes an unboxed tuple. A boxed tuple is written '() and
would be put in the heap. An unboxed tuple may also be written
without the @-sign, but then the compiler might confuse it with,
for instance, an argument list. (ATS2 has conventional argument
lists that are distinct from tuples, and supports
call-by-reference, where an argument is mutable.) *)
typedef ng8 = @(exrat, exrat, exrat, exrat,
exrat, exrat, exrat, exrat)
typedef continued_fraction_binary_op_cloref =
(continued_fraction, continued_fraction) -<cloref1> continued_fraction
(* ng8_make_int takes ONE argument, which is a tuple. *)
val ng8_make_int : @(int, int, int, int, int, int, int, int) -> ng8
val ng8_stopping_processing_threshold : ref exrat
val ng8_infinitization_threshold : ref exrat
val ng8_apply : ng8 -> continued_fraction_binary_op_cloref
val ng8_apply_sub : continued_fraction_binary_op_cloref
val ng8_apply_mul : continued_fraction_binary_op_cloref
val ng8_apply_div : continued_fraction_binary_op_cloref
(* The following two are regular functions, not closures. They are
translated by the ATS compiler into ordinary C functions. *)
fn ng8_apply_neg : continued_fraction -> continued_fraction
fn ng8_apply_pow : (continued_fraction, int) -> continued_fraction

(* Miscellanous continued fractions. *)
val zero : continued_fraction
val one : continued_fraction
val two : continued_fraction
val three : continued_fraction
val four : continued_fraction
//
val one_fourth : continued_fraction
val one_third : continued_fraction
val one_half : continued_fraction
val two_thirds : continued_fraction
val three_fourths : continued_fraction
//
val golden_ratio : continued_fraction
val silver_ratio : continued_fraction
val sqrt2 : continued_fraction
val sqrt5 : continued_fraction

The second file is an implementation of the stuff declared in the first file. The second file is called continued_fraction.dats:

(* "Dynamic" file. (Implementations.) *)

(* To set up a predictable name-mangling scheme:  *)
#define ATS_PACKNAME "rosetta-code.continued_fraction"

(* Load templates from the ATS2 prelude: *)

(* Load declarations and templates from ats2-xprelude: *)
#include "xprelude/HATS/xprelude.hats"

(* Load the declarations for this package: *)

typedef cf_record =
(* A cf_record is an unboxed record, denoted by @{}. A boxed record
would be written '{} and would be placed in the heap. Either way,
it is an immutable record. For a mutable record, we would have to
use vtypedef to make it a LINEAR type. *)
@{
terminated = bool,          (* Is the generator exhausted? *)
memo_count = size_t,        (* How many terms are memoized? *)

(* An arrszref is an array with runtime bounds checking. An
arrszref is less efficient than an arrayref, but will not force
us to use dependent types for the indices. *)
memo = arrszref exrat,      (* Memoized terms. *)

generate = term_generator   (* The source of terms. *)
}

(* The actual type of a continued_fraction is a MUTABLE reference to
the (immutable) cf_record. Within this file, we may also call the
type cf_t. *)
typedef cf_t = ref cf_record
assume continued_fraction = cf_t

implement
continued_fraction_make generator =
let
val record : cf_record =
@{
terminated = false,
memo_count = i2sz 0,
memo = arrszref_make_elt<exrat> (i2sz 32, exrat_make (0, 1)),
generate = generator
}
in
ref record
end

(* "fn" means a non-recursive function. A function that might be
recursive is written "fun" (or sometimes "fnx"). Incidentally: it
is common to see the recursions put into nested functions, with the
function a programmer is supposed to call being non-recursive. This
is often a matter of style. (By the way, in a "*.sats" file there
is no distinction between "fn" and "fun" that I know of.) *)
fn
resize_if_necessary (cf : cf_t, i : size_t) : void =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if size memo <= i then
let
val new_size = i2sz 2 * (succ i)
val new_memo =
arrszref_make_elt<exrat> (new_size, exrat_make (0, 1))
val new_record : cf_record =
@{
terminated = terminated,
memo_count = memo_count,
memo = new_memo,
generate = generate
}

var i : size_t          (* A C-style automatic variable. *)
in
(* A C-style for-loop. *)
for (i := i2sz 0; i <> memo_count; i := succ i)
new_memo[i] := memo[i];

!cf := new_record
end
end

fn
update_terms (cf : cf_t, i : size_t) : void =
let
fun
loop () : void =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if terminated then
()
else if i < memo_count then
()
else
case generate () of
| None () =>
let
val new_record =
@{
terminated = true,
memo_count = memo_count,
memo = memo,
generate = generate
}
in
!cf := new_record
end
| Some term =>
(* "begin-end" is a synonym for "()". *)
begin
memo[memo_count] := term;
let
val new_record =
@{
terminated = false,
memo_count = succ memo_count,
memo = memo,
generate = generate
}
in
!cf := new_record;
loop ()
end
end
end
in
loop ()
end

implement
term_exists (cf, i) =
let
val i = i2sz i

fun
loop () =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if i < memo_count then
true
else if terminated then
false
else
begin
resize_if_necessary (cf, i);
update_terms (cf, i);
loop ()
end
end
in
loop ()
end

implement
get_term_exn (cf, i) =
if i2sz i < (!cf).memo_count then
let
val memo = (!cf).memo
in
memo[i]
end
else
$raise IllegalArgExn "get_term_exn:out_of_bounds" implement continued_fraction_to_term_generator cf = let val i : ref (intGte 0) = ref 0 in lam () =<cloref1> let val j = !i in if term_exists (cf, j) then begin !i := succ j; Some (cf[j]) end else None () end end implement default_max_terms = ref 20 implement continued_fraction_to_string_given_max_terms (cf, max_terms) = let fun loop (i : intGte 0, accum : string) : string = if ~term_exists (cf, i) then (* The return value of string_append is a LINEAR, MUTABLE strptr, which we cast to a nonlinear, immutable string. (One could introduce one's own shorthands, though.) *) strptr2string (string_append (accum, "]")) else if i = max_terms then strptr2string (string_append (accum, ",...]")) else let val separator = if i = 0 then "" else if i = 1 then ";" else "," and term_string = tostring_val<exrat> cf[i] in loop (succ i, strptr2string (string_append (accum, separator, term_string))) end in loop (0, "[") end implement continued_fraction_to_string_default_max_terms cf = let val max_terms = !default_max_terms in continued_fraction_to_string_given_max_terms (cf, max_terms) end implement int_to_continued_fraction i = let val done : ref bool = ref false val i = (g0i2f i) : exrat in continued_fraction_make (lam () =<cloref1> if !done then None () else begin !done := true; Some i end) end implement exrat_to_continued_fraction num = let val done : ref bool = ref false val num : ref exrat = ref num in continued_fraction_make (lam () =<cloref1> if !done then None () else let val q = floor !num val r = !num - q in if iseqz r then !done := true else !num := reciprocal r; Some q end) end implement rational_to_continued_fraction (numer, denom) = exrat_to_continued_fraction (exrat_make (numer, denom)) implement continued_fraction_make_constant_term i = let val i = (g0i2f i) : exrat in continued_fraction_make (lam () =<cloref1> Some i) end implement ng8_make_int tuple = let val @(a12, a1, a2, a, b12, b1, b2, b) = tuple fn f (i : int) : exrat = exrat_make (i, 1) in @(f a12, f a1, f a2, f a, f b12, f b1, f b2, f b) end implement ng8_stopping_processing_threshold = ref (exrat_make (2, 1) ** 512) implement ng8_infinitization_threshold = ref (exrat_make (2, 1) ** 64) fn too_big (term : exrat) : bool = abs (term) >= abs (!ng8_stopping_processing_threshold) fn any_too_big (ng : ng8) : bool = (* "orelse" may also be (and usually is) written "||", as in C. The "orelse" notation resembles that of Standard ML. Non-shortcircuiting OR also exists, and can be written "+". *) case+ ng of (* <-- the + sign means all cases must have a match. *) | @(a, b, c, d, e, f, g, h) => too_big (a) orelse too_big (b) orelse too_big (c) orelse too_big (d) orelse too_big (e) orelse too_big (f) orelse too_big (g) orelse too_big (h) fn infinitize (term : exrat) : Option exrat = if abs (term) >= abs (!ng8_infinitization_threshold) then None () else Some term val no_terms_source : term_generator = lam () =<cloref1> None () fn divide (a : exrat, b : exrat) : @(exrat, exrat) = if iseqz b then @(exrat_make (0, 1), exrat_make (0, 1)) else (* Do integer division of the numerators of a and b. The following particular function does floor division if the divisor is positive, ceiling division if the divisor is negative. Thus the remainder is never negative. *) exrat_numerator_euclid_division (a, b) implement ng8_apply ng = lam (x, y) => let val ng : ref ng8 = ref ng and xsource : ref term_generator = ref (continued_fraction_to_term_generator x) and ysource : ref term_generator = ref (continued_fraction_to_term_generator y) fn all_b_are_zero () : bool = let val @(_, _, _, _, b12, b1, b2, b) = !ng in (* Instead of the Standard ML-like notation "andalso", one may (and usually does) use the C-like notation "&&". There is also non-shortcircuiting AND, written "*". *) iseqz b andalso iseqz b2 andalso iseqz b1 andalso iseqz b12 end fn all_four_equal (a : exrat, b : exrat, c : exrat, d : exrat) : bool = a = b && a = c && a = d fn absorb_x_term () = let val @(a12, a1, a2, a, b12, b1, b2, b) = !ng in case (!xsource) () of | Some term => let val new_ng = (a2 + (a12 * term), a + (a1 * term), a12, a1, b2 + (b12 * term), b + (b1 * term), b12, b1) in if any_too_big new_ng then (* Pretend all further x terms are infinite. *) (!ng := @(a12, a1, a12, a1, b12, b1, b12, b1); !xsource := no_terms_source) else !ng := new_ng end | None () => !ng := @(a12, a1, a12, a1, b12, b1, b12, b1) end fn absorb_y_term () = let val @(a12, a1, a2, a, b12, b1, b2, b) = !ng in case (!ysource) () of | Some term => let val new_ng = (a1 + (a12 * term), a12, a + (a2 * term), a2, b1 + (b12 * term), b12, b + (b2 * term), b2) in if any_too_big new_ng then (* Pretend all further y terms are infinite. *) (!ng := @(a12, a12, a2, a2, b12, b12, b2, b2); !ysource := no_terms_source) else !ng := new_ng end | None () => !ng := @(a12, a12, a2, a2, b12, b12, b2, b2) end fun loop () = (* ATS2 can do mutual recursion with proper tail calls, but, to stay closer to the Standard ML code, here I use only single tail recursion. To do mutual recursion with proper tail calls, one says "fnx" instead of "fun". *) if all_b_are_zero () then None () (* There are no more terms to output. *) else let val @(_, _, _, _, b12, b1, b2, b) = !ng in if iseqz b andalso iseqz b2 then (absorb_x_term (); loop ()) else if iseqz b orelse iseqz b2 then (absorb_y_term (); loop ()) else if iseqz b1 then (absorb_x_term (); loop ()) else let val @(a12, a1, a2, a, _, _, _, _) = !ng val @(q12, r12) = divide (a12, b12) and @(q1, r1) = divide (a1, b1) and @(q2, r2) = divide (a2, b2) and @(q, r) = divide (a, b) in if isneqz b12 andalso all_four_equal (q12, q1, q2, q) then (!ng := (b12, b1, b2, b, r12, r1, r2, r); (* Return a term--or, if a magnitude threshold is reached, return no more terms . *) infinitize q) else let (* Put a1, a2, and a over a common denominator and compare some magnitudes. (SIDE NOTE: We are representing big integers as EXACT rationals with denominator one, so in fact could have put a1, a2, and a over their respective denominators and compared the fractions. However, I have retained the phrasing of the Standard ML program.) *) val n1 = a1 * b2 * b and n2 = a2 * b1 * b and n = a * b1 * b2 in if abs (n1 - n) > abs (n2 - n) then (absorb_x_term (); loop ()) else (absorb_y_term (); loop ()) end end end in continued_fraction_make (lam () =<cloref1> loop ()) end (* A macro definition: *) macdef make_op (tuple) = ng8_apply (ng8_make_int ,(tuple)) implement ng8_apply_add = make_op @(0, 1, 1, 0, 0, 0, 0, 1) implement ng8_apply_sub = make_op @(0, 1, ~1, 0, 0, 0, 0, 1) implement ng8_apply_mul = make_op @(1, 0, 0, 0, 0, 0, 0, 1) implement ng8_apply_div = make_op @(0, 1, 0, 0, 0, 0, 1, 0) (* Here the closure is "wrapped" in an ordinary function. *) val _ng8_apply_neg = make_op @(0, 0, ~1, 0, 0, 0, 0, 1) implement ng8_apply_neg cf = _ng8_apply_neg (cf, cf) val _reciprocal = make_op @(0, 0, 0, 1, 0, 1, 0, 0) implement ng8_apply_pow (cf, i) = let macdef reciprocal cf = _reciprocal (,(cf), ,(cf)) fun loop (x : continued_fraction, n : int, accum : continued_fraction) : continued_fraction = if 1 < n then let val nhalf = n / 2 and xsquare = x * x in if nhalf + nhalf <> n then loop (xsquare, nhalf, accum * x) else loop (xsquare, nhalf, accum) end else if n = 1 then accum * x else accum in if 0 <= i then loop (cf, i, one) else reciprocal (loop (cf, ~i, one)) end implement zero = i2cf 0 implement one = i2cf 1 implement two = i2cf 2 implement three = i2cf 3 implement four = i2cf 4 implement one_fourth = r2cf (1, 4) implement one_third = r2cf (1, 3) implement one_half = r2cf (1, 2) implement two_thirds = r2cf (2, 3) implement three_fourths = r2cf (3, 4) implement golden_ratio = constant_term_cf 1 implement silver_ratio = constant_term_cf 2 implement sqrt2 = silver_ratio - one implement sqrt5 = (two * golden_ratio) - one The third file is the main program. It is called continued-fraction-task.dats: (* Main program. *) (* Install ats2-xprelude, being sure to enable GMP support: https://sourceforge.net/p/chemoelectric/ats2-xprelude If you have it installed already, there might have been bugfixes since. So try updating. Then, to compile the program: patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW \$(pkg-config --cflags ats2-xprelude) \
$(pkg-config --variable=PATSCCFLAGS ats2-xprelude) \ continued-fraction-task.dats continued_fraction.{s,d}ats \$(pkg-config --libs ats2-xprelude) -lgc -lm

*)

(* Load templates from the ATS2 prelude: *)

dynload "continued_fraction.dats" (* Initialize the "val". *)

(* Load declarations and templates from ats2-xprelude: *)
#include "xprelude/HATS/xprelude.hats"

fn
make_pad (n : size_t) : string =
let
val n = g1ofg0 n
prval () = lemma_g1uint_param n
implement string_tabulate$fopr<> _ = ' ' in strnptr2string (string_tabulate<> n) end fn show_with_note (expression : string, cf : continued_fraction, note : string) : void = let val cf_str = cf2string cf val expr_sz = strlen expression and cf_sz = strlen cf_str and note_sz = strlen note val expr_pad_sz = max (i2sz 19 - expr_sz, i2sz 0) and cf_pad_sz = if iseqz note_sz then i2sz 0 else max (i2sz 48 - cf_sz, i2sz 0) val expr_pad = make_pad expr_pad_sz and cf_pad = make_pad cf_pad_sz in println! (expr_pad, expression, " => ", cf_str, cf_pad, note) end fn show_without_note (expression : string, cf : continued_fraction) : void = show_with_note (expression, cf, "") overload show with show_with_note overload show with show_without_note implement main0 () = (* A main that takes no arguments and returns 0. *) begin show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2"); show ("silver ratio", silver_ratio, "(1 + sqrt(2))"); show ("sqrt2", sqrt2); show ("sqrt5", sqrt5); show ("1/4", one_fourth); show ("1/3", one_third); show ("1/2", one_half); show ("2/3", two_thirds); show ("3/4", three_fourths); show ("13/11", r2cf (13, 11)); show ("22/7", r2cf (22, 7), "approximately pi"); show ("0", zero); show ("1", one); show ("2", two); show ("3", three); show ("4", four); show ("4 + 3", four + three); show ("4 - 3", four - three); show ("4 * 3", four * three); show ("4 / 3", four / three); show ("4 ** 3", four ** 3); show ("4 ** (-3)", four ** (~3)); show ("negative 4", ~four); show ("(1 + 1/sqrt(2))/2", (one + (one / sqrt2)) / two, "method 1"); show ("(1 + 1/sqrt(2))/2", silver_ratio * (sqrt2 ** (~3)), "method 2"); show ("(1 + 1/sqrt(2))/2", ((silver_ratio ** 2) + one) / (four * two), "method 3"); show ("sqrt2 + sqrt2", sqrt2 + sqrt2); show ("sqrt2 - sqrt2", sqrt2 - sqrt2); show ("sqrt2 * sqrt2", sqrt2 * sqrt2); show ("sqrt2 / sqrt2", sqrt2 / sqrt2); end Output: To compile the program, you might try something like the following (assuming you have Boehm GC): patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW$(pkg-config --cflags ats2-xprelude) $(pkg-config --variable=PATSCCFLAGS ats2-xprelude) continued-fraction-task.dats continued_fraction.{s,d}ats$(pkg-config --libs ats2-xprelude) -lgc -lm

You have to specify some C language standard, because patscc defaults to C99.

Then run the program by typing
./a.out

The output should resemble that of the Standard ML program from which the ATS was translated. Minus signs might look different:

       golden ratio =>  [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]   (1 + sqrt(5))/2
silver ratio =>  [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]   (1 + sqrt(2))
sqrt2 =>  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt5 =>  [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...]
1/4 =>  [0;4]
1/3 =>  [0;3]
1/2 =>  [0;2]
2/3 =>  [0;1,2]
3/4 =>  [0;1,3]
13/11 =>  [1;5,2]
22/7 =>  [3;7]                                           approximately pi
0 =>  [0]
1 =>  [1]
2 =>  [2]
3 =>  [3]
4 =>  [4]
4 + 3 =>  [7]
4 - 3 =>  [1]
4 * 3 =>  [12]
4 / 3 =>  [1;3]
4 ** 3 =>  [64]
4 ** (-3) =>  [0;64]
negative 4 =>  [-4]
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 1
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 2
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 3
sqrt2 + sqrt2 =>  [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt2 - sqrt2 =>  [0]
sqrt2 * sqrt2 =>  [2]
sqrt2 / sqrt2 =>  [1]


C

Translation of: Python
Translation of: Scheme

You will need Boehm GC and the GNU Multiple Precision Library.

(Actually you can leave out the Boehm GC parts. The program leaks memory, but harmlessly. Also, you can leave out the C23 attribute specifiers.)

/*------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <gc/gc.h>              /* Boehm GC. */
#include <gmp.h>                /* GNU Multiple Precision. */

void *
my_malloc (size_t size)
{
void *p = GC_MALLOC (size);
if (p == NULL)
{
fprintf (stderr, "Memory allocation failed.\n");
exit (1);
}
return p;
}

void *
my_realloc (void *p, size_t size)
{
void *q = GC_REALLOC (p, size);
if (q == NULL)
{
fprintf (stderr, "Memory allocation failed.\n");
exit (1);
}
return q;
}

void
my_free (void *p)
{
GC_FREE (p);
}

/*------------------------------------------------------------------*/
/* Some helper functions. */

#define MIN(x, y) (((x) < (y)) ? (x) : (y))
#define MAX(x, y) (((x) > (y)) ? (x) : (y))

char *
string_append1 (const char *s1)
{
size_t n1 = strlen (s1);
char *s = my_malloc ((n1 + 1) * sizeof (char));
s[n1] = '\0';
memcpy (s, s1, n1);
return s;
}

char *
string_append2 (const char *s1, const char *s2)
{
size_t n1 = strlen (s1);
size_t n2 = strlen (s2);
char *s = my_malloc ((n1 + n2 + 1) * sizeof (char));
s[n1 + n2] = '\0';
memcpy (s, s1, n1);
memcpy (s + n1, s2, n2);
return s;
}

char *
string_append3 (const char *s1, const char *s2, const char *s3)
{
size_t n1 = strlen (s1);
size_t n2 = strlen (s2);
size_t n3 = strlen (s3);
char *s = my_malloc ((n1 + n2 + n3 + 1) * sizeof (char));
s[n1 + n2 + n3] = '\0';
memcpy (s, s1, n1);
memcpy (s + n1, s2, n2);
memcpy (s + n1 + n2, s3, n3);
return s;
}

char *
string_repeat (size_t n, const char *s)
{
/* This brute force implementation will suffice. */
char *t = "";
for (size_t i = 0; i != n; i += 1)
t = string_append2 (t, s);
return t;
}

/*------------------------------------------------------------------*/

typedef mpz_t *cf_func (size_t i, void *env);

struct cf
{
bool terminated;              /* Are there no more terms? */
size_t m;                     /* The number of terms memoized. */
size_t n;                     /* The size of memoization storage. */
mpz_t **memo;                 /* Memoization storage. */
cf_func *func;                /* A function that produces terms. */
void *env;                    /* An environment for func. */
};

typedef struct cf *cf_t;

cf_t
make_cf (cf_func *func, void *env)
{
cf_t cf = my_malloc (sizeof (struct cf));
cf->terminated = false;
cf->m = 0;
cf->n = 32;
cf->memo = my_malloc (cf->n * sizeof (mpz_t *));
cf->func = func;
cf->env = env;
return cf;
}

void
resize_cf (cf_t cf, size_t minimum)
{
/* Ensure there is at least twice the minimum storage. */
size_t size = 2 * minimum;
if (cf->n < size)
{
cf->memo = my_realloc (cf->memo, size * sizeof (mpz_t *));
cf->n = size;
}
}

void
update_cf (cf_t cf, size_t needed)
{
/* Ensure there are at least a certain number of finite terms
memoized (or else that all of them are memoized). */
if (!cf->terminated && cf->m < needed)
{
if (cf->n < needed)
resize_cf (cf, needed);
while (!cf->terminated && cf->m != needed)
{
cf->memo[cf->m] = cf->func (cf->m, cf->env);
cf->m += 1;
}
}
}

mpz_t *
cf_ref (cf_t cf, size_t i)
{
/* Get the ith term, or a NULL pointer if there is no finite ith
term. */
update_cf (cf, i + 1);
return (i < cf->m) ? cf->memo[i] : NULL;
}

size_t default_max_terms = 20;

char *
cf2string (cf_t cf, size_t max_terms)
{
if (max_terms == 0)
max_terms = default_max_terms;

size_t i = 0;
char *s = string_append1 ("[");
bool done = false;
while (!done)
{
mpz_t *term = cf_ref (cf, i);
if (term == NULL)
{
s = string_append2 (s, "]");
done = true;
}
else if (i == max_terms)
{
s = string_append2 (s, ",...]");
done = true;
}
else
{
static const char *separators[3] = { "", ";", "," };
const char *separator = separators[(i <= 1) ? i : 2];
const char *term_str = mpz_get_str (NULL, 10, *term);
s = string_append3 (s, separator, term_str);
i += 1;
}
}
return s;
}

/*------------------------------------------------------------------*/

cf_t golden_ratio;
cf_t silver_ratio;

mpz_t *
return_constant ([[maybe_unused]] size_t i, void *env)
{
mpz_t *term = my_malloc (sizeof (mpz_t));
mpz_init_set (*term, *((mpz_t *) env));
return term;
}

cf_t
make_cf_with_constant_terms (int term_si)
{
mpz_t *env = my_malloc (sizeof (mpz_t));
mpz_init_set_si (*env, term_si);
return make_cf (return_constant, env);
}

/*------------------------------------------------------------------*/

cf_t sqrt2;

mpz_t *
return_sqrt2_term (size_t i, [[maybe_unused]] void *env)
{
mpz_t *term = my_malloc (sizeof (mpz_t));
mpz_init_set_si (*term, (i == 0) ? 1 : 2);
return term;
}

cf_t
make_cf_sqrt2 (void)
{
return make_cf (return_sqrt2_term, NULL);
}

/*------------------------------------------------------------------*/

cf_t frac_13_11;
cf_t frac_22_7;
cf_t one;
cf_t two;
cf_t three;
cf_t four;

mpz_t *
return_rational_term ([[maybe_unused]] size_t i, void *env)
{
mpz_t *frac = env;
mpz_t *term = NULL;
if (mpz_sgn (frac[1]) != 0)
{
term = my_malloc (sizeof (mpz_t));
mpz_init (*term);
mpz_t r;
mpz_init (r);
mpz_fdiv_qr (*term, r, frac[0], frac[1]);
mpz_set (frac[0], frac[1]);
mpz_set (frac[1], r);
}
return term;
}

cf_t
make_cf_rational (int numerator_si, int denominator_si)
{
mpz_t *env = my_malloc (2 * sizeof (mpz_t));
mpz_init_set_si (env[0], numerator_si);
mpz_init_set_si (env[1], denominator_si);
return make_cf (return_rational_term, env);
}

cf_t
make_cf_integer (int integer_si)
{
return make_cf_rational (integer_si, 1);
}

/*------------------------------------------------------------------*/

/* Thresholds. */
mpz_t number_that_is_too_big;
mpz_t practically_infinite;

struct ng8_env
{
mpz_t ng[8];
cf_t x;
cf_t y;
size_t ix;
size_t iy;
bool xoverflow;
bool yoverflow;
};

typedef struct ng8_env *ng8_env_t;

enum ng8_index
{
ng8a12 = 0,
ng8a1  = 1,
ng8a2  = 2,
ng8a   = 3,
ng8b12 = 4,
ng8b1  = 5,
ng8b2  = 6,
ng8b   = 7
};

static bool
ng8_too_big (const mpz_t ng[8])
{
bool too_big = false;
int i = 0;
while (!too_big && i != 8)
{
too_big = (0 <= mpz_cmpabs (ng[i], number_that_is_too_big));
i += 1;
}
}

static bool
treat_ng8_term_as_infinite (const mpz_t term)
{
return (0 <= mpz_cmpabs (term, practically_infinite));
}

static void
a_plus_bc (mpz_t result, const mpz_t a,  const mpz_t b,
const mpz_t c)
{
mpz_set (result, a);
}

static void
abc (mpz_t result, const mpz_t a,  const mpz_t b, const mpz_t c)
{
mpz_mul (result, a, b);
mpz_mul (result, result, c);
}

static void
absorb_x_term (ng8_env_t env)
{
mpz_t tmp[8];
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->xoverflow) ? cf_ref (env->x, env->ix) : NULL;
env->ix += 1;
mpz_set (env->ng[ng8a2], tmp[ng8a12]);
mpz_set (env->ng[ng8a], tmp[ng8a1]);
mpz_set (env->ng[ng8b2], tmp[ng8b12]);
mpz_set (env->ng[ng8b], tmp[ng8b1]);
if (term != NULL)
{
a_plus_bc (env->ng[ng8a12], tmp[ng8a2], tmp[ng8a12], *term);
a_plus_bc (env->ng[ng8a1], tmp[ng8a], tmp[ng8a1], *term);
a_plus_bc (env->ng[ng8b12], tmp[ng8b2], tmp[ng8b12], *term);
a_plus_bc (env->ng[ng8b1], tmp[ng8b], tmp[ng8b1], *term);
if (ng8_too_big (env->ng))
{
env->xoverflow = true;
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a1], tmp[ng8a1]);
mpz_set (env->ng[ng8b12], tmp[ng8b12]);
mpz_set (env->ng[ng8b1], tmp[ng8b1]);
}
}
}

static void
absorb_y_term (ng8_env_t env)
{
mpz_t tmp[8];
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->yoverflow) ? cf_ref (env->y, env->iy) : NULL;
env->iy += 1;
mpz_set (env->ng[ng8a1], tmp[ng8a12]);
mpz_set (env->ng[ng8a], tmp[ng8a2]);
mpz_set (env->ng[ng8b1], tmp[ng8b12]);
mpz_set (env->ng[ng8b], tmp[ng8b2]);
if (term != NULL)
{
a_plus_bc (env->ng[ng8a12], tmp[ng8a1], tmp[ng8a12], *term);
a_plus_bc (env->ng[ng8a2], tmp[ng8a], tmp[ng8a2], *term);
a_plus_bc (env->ng[ng8b12], tmp[ng8b1], tmp[ng8b12], *term);
a_plus_bc (env->ng[ng8b2], tmp[ng8b], tmp[ng8b2], *term);
if (ng8_too_big (env->ng))
{
env->yoverflow = true;
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a2], tmp[ng8a2]);
mpz_set (env->ng[ng8b12], tmp[ng8b12]);
mpz_set (env->ng[ng8b2], tmp[ng8b2]);
}
}
}

mpz_t *
return_ng8_term ([[maybe_unused]] size_t i, void *env)
{
/* The McCabe complexity of this function is high. Please be careful
if modifying the code. */

ng8_env_t p = env;

mpz_t *term = NULL;

bool done = false;
while (!done)
{
const bool b12_zero = (mpz_sgn (p->ng[ng8b12]) == 0);
const bool b1_zero = (mpz_sgn (p->ng[ng8b1]) == 0);
const bool b2_zero = (mpz_sgn (p->ng[ng8b2]) == 0);
const bool b_zero = (mpz_sgn (p->ng[ng8b]) == 0);

if (b_zero && b1_zero && b2_zero && b12_zero)
done = true;            /* There are no more terms. */
else if (b_zero && b2_zero)
absorb_x_term (p);
else if (b_zero || b2_zero)
absorb_y_term (p);
else if (b1_zero)
absorb_x_term (p);
else
{
mpz_t q, r;
mpz_inits (q, r, NULL);
mpz_t q1, r1;
mpz_inits (q1, r1, NULL);
mpz_t q2, r2;
mpz_inits (q2, r2, NULL);
mpz_t q12, r12;
mpz_inits (q12, r12, NULL);

mpz_fdiv_qr (q, r, p->ng[ng8a], p->ng[ng8b]);
mpz_fdiv_qr (q1, r1, p->ng[ng8a1], p->ng[ng8b1]);
mpz_fdiv_qr (q2, r2, p->ng[ng8a2], p->ng[ng8b2]);
if (!b12_zero)
mpz_fdiv_qr (q12, r12, p->ng[ng8a12], p->ng[ng8b12]);

if (!b12_zero
&& mpz_cmp (q, q1) == 0
&& mpz_cmp (q, q2) == 0
&& mpz_cmp (q, q12) == 0)
{
// Output a term.
mpz_set (p->ng[ng8a12], p->ng[ng8b12]);
mpz_set (p->ng[ng8a1],  p->ng[ng8b1]);
mpz_set (p->ng[ng8a2],  p->ng[ng8b2]);
mpz_set (p->ng[ng8a],   p->ng[ng8b]);
mpz_set (p->ng[ng8b12], r12);
mpz_set (p->ng[ng8b1],  r1);
mpz_set (p->ng[ng8b2],  r2);
mpz_set (p->ng[ng8b],   r);
if (!treat_ng8_term_as_infinite (q))
{
term = my_malloc (sizeof (mpz_t));
mpz_init_set (*term, q);
}
done = true;
}
else
{
/* Rather than compare fractions, we will put the
numerators over a common denominator of b*b1*b2, and
then compare the new numerators. */
mpz_t n, n1, n2, n1_diff, n2_diff;
mpz_inits (n, n1, n2, n1_diff, n2_diff, NULL);
abc (n, p->ng[ng8a], p->ng[ng8b1], p->ng[ng8b2]);
abc (n1, p->ng[ng8a1], p->ng[ng8b], p->ng[ng8b2]);
abc (n2, p->ng[ng8a2], p->ng[ng8b], p->ng[ng8b1]);
mpz_sub (n1_diff, n1, n);
mpz_sub (n2_diff, n2, n);
if (mpz_cmpabs (n1_diff, n2_diff) > 0)
absorb_x_term (p);
else
absorb_y_term (p);
}
}
}

return term;
}

cf_t
make_cf_ng8 (int ng[8], cf_t x, cf_t y)
{
ng8_env_t env = my_malloc (sizeof (struct ng8_env));
for (int i = 0; i != 8; i += 1)
mpz_init_set_si (env->ng[i], ng[i]);
env->x = x;
env->y = y;
env->ix = 0;
env->iy = 0;
env->xoverflow = false;
env->yoverflow = false;
return make_cf (return_ng8_term, env);
}

/*------------------------------------------------------------------*/

static void *
gmp_malloc (size_t alloc_size)
{
return my_malloc (alloc_size);
}

static void *
gmp_realloc (void *p,
[[maybe_unused]] size_t old_size,
size_t new_size)
{
return my_realloc (p, new_size);
}

static void
gmp_free (void *p, [[maybe_unused]] size_t size)
{
/* There is no need for us to explicitly free memory, and
performance might even suffer if we do. On the other hand, maybe
GMP will free memory that otherwise would have been passed over
for collection. */
my_free (p);                  /* <-- optional */
}

void
show (const char *expression, cf_t cf, const char *note)
{
size_t nexpr = strlen (expression);
char *padding = string_repeat (MAX (19, nexpr + 1) - nexpr, " ");
char *line = string_append3 (padding, expression, " =>  ");
char *cfstr = cf2string (cf, 0);
line = string_append2 (line, cfstr);
if (note != NULL)
{
size_t ncfstr = strlen (cfstr);
padding = string_repeat (MAX (48, ncfstr + 1) - ncfstr, " ");
line = string_append3 (line, padding, note);
}
puts (line);
}

int ng8_add[8] = { 0, 1, 1, 0, 0, 0, 0, 1 };
int ng8_sub[8] = { 0, 1, -1, 0, 0, 0, 0, 1 };
int ng8_mul[8] = { 1, 0, 0, 0, 0, 0, 0, 1 };
int ng8_div[8] = { 0, 1, 0, 0, 0, 0, 1, 0 };

int
main (void)
{
GC_INIT ();

/* GMP has to be told to use Boehm GC as its allocator. */
mp_set_memory_functions (gmp_malloc, gmp_realloc, gmp_free);

/* Initialize thresholds, to values chosen merely for
demonstration. */
mpz_init_set_si (number_that_is_too_big, 1);
mpz_mul_2exp (number_that_is_too_big, number_that_is_too_big,
512);           /* 2**512 */
mpz_init_set_si (practically_infinite, 1);
mpz_mul_2exp (practically_infinite, practically_infinite,
64);            /* 2**64 */

/* Initialize global continued fractions. */
golden_ratio = make_cf_with_constant_terms (1);
silver_ratio = make_cf_with_constant_terms (2);
sqrt2 = make_cf_sqrt2 ();
frac_13_11 = make_cf_rational (13, 11);
frac_22_7 = make_cf_rational (22, 7);
one = make_cf_integer (1);
two = make_cf_integer (2);
three = make_cf_integer (3);
four = make_cf_integer (4);

/* Divide the silver ratio by 2 times the square root of 2. */
int ng8_method1[8] = { 0, 1, 0, 0, 0, 0, 2, 0 };
cf_t method1 = make_cf_ng8 (ng8_method1, silver_ratio, sqrt2);

/* Add 1/8 to 1/8th of the square of the silver ratio. */
int ng8_method2[8] = { 1, 0, 0, 1, 0, 0, 0, 8 };
cf_t method2 = make_cf_ng8 (ng8_method2, silver_ratio,
silver_ratio);

/* Thrice divide the silver ratio by the square root of 2. */
cf_t method3 = make_cf_ng8 (ng8_div, silver_ratio, sqrt2);
method3 = make_cf_ng8 (ng8_div, method3, sqrt2);
method3 = make_cf_ng8 (ng8_div, method3, sqrt2);

show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "1 + sqrt(2)");
show ("sqrt(2)", sqrt2, NULL);
show ("13/11", frac_13_11, NULL);
show ("22/7", frac_22_7, NULL);
show ("one", one, NULL);
show ("two", two, NULL);
show ("three", three, NULL);
show ("four", four, NULL);
show ("(1 + 1/sqrt(2))/2", method1, "method 1");
show ("(1 + 1/sqrt(2))/2", method2, "method 2");
show ("(1 + 1/sqrt(2))/2", method3, "method 3");
show ("sqrt(2) + sqrt(2)", make_cf_ng8 (ng8_add, sqrt2, sqrt2),
NULL);
show ("sqrt(2) - sqrt(2)", make_cf_ng8 (ng8_sub, sqrt2, sqrt2),
NULL);
show ("sqrt(2) * sqrt(2)", make_cf_ng8 (ng8_mul, sqrt2, sqrt2),
NULL);
show ("sqrt(2) / sqrt(2)", make_cf_ng8 (ng8_div, sqrt2, sqrt2),
NULL);

return 0;
}

/*------------------------------------------------------------------*/

Output:
$gcc -std=gnu2x -Wall -Wextra -g bivariate-continued-fraction-task-gmp.c -lgmp -lgc && ./a.out golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2 silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 1 + sqrt(2) sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 13/11 => [1;5,2] 22/7 => [3;7] one => [1] two => [2] three => [3] four => [4] (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3 sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] sqrt(2) - sqrt(2) => [0] sqrt(2) * sqrt(2) => [2] sqrt(2) / sqrt(2) => [1]  C++ Uses matrixNG, NG_4 and NG from Continued_fraction/Arithmetic/G(matrix_NG,_Contined_Fraction_N)#C++, and r2cf from Continued_fraction/Arithmetic/Construct_from_rational_number#C++ /* Implement matrix NG Nigel Galloway, February 12., 2013 */ class NG_8 : public matrixNG { private: int a12, a1, a2, a, b12, b1, b2, b, t; double ab, a1b1, a2b2, a12b12; const int chooseCFN(){return fabs(a1b1-ab) > fabs(a2b2-ab)? 0 : 1;} const bool needTerm() { if (b1==0 and b==0 and b2==0 and b12==0) return false; if (b==0){cfn = b2==0? 0:1; return true;} else ab = ((double)a)/b; if (b2==0){cfn = 1; return true;} else a2b2 = ((double)a2)/b2; if (b1==0){cfn = 0; return true;} else a1b1 = ((double)a1)/b1; if (b12==0){cfn = chooseCFN(); return true;} else a12b12 = ((double)a12)/b12; thisTerm = (int)ab; if (thisTerm==(int)a1b1 and thisTerm==(int)a2b2 and thisTerm==(int)a12b12){ t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm; t=a2; a2=b2; b2=t-b2*thisTerm; t=a12; a12=b12; b12=t-b12*thisTerm; haveTerm = true; return false; } cfn = chooseCFN(); return true; } void consumeTerm(){if(cfn==0){a=a1; a2=a12; b=b1; b2=b12;} else{a=a2; a1=a12; b=b2; b1=b12;}} void consumeTerm(int n){ if(cfn==0){t=a; a=a1; a1=t+a1*n; t=a2; a2=a12; a12=t+a12*n; t=b; b=b1; b1=t+b1*n; t=b2; b2=b12; b12=t+b12*n;} else{t=a; a=a2; a2=t+a2*n; t=a1; a1=a12; a12=t+a12*n; t=b; b=b2; b2=t+b2*n; t=b1; b1=b12; b12=t+b12*n;} } public: NG_8(int a12, int a1, int a2, int a, int b12, int b1, int b2, int b): a12(a12), a1(a1), a2(a2), a(a), b12(b12), b1(b1), b2(b2), b(b){ }};  Testing [3;7] + [0;2] int main() { NG_8 a(0,1,1,0,0,0,0,1); r2cf n2(22,7); r2cf n1(1,2); for(NG n(&a, &n1, &n2); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; NG_4 a3(2,1,0,2); r2cf n3(22,7); for(NG n(&a3, &n3); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; return 0; }  Output: 3 1 1 1 4 3 1 1 1 4  [1:5,2] * [3;7] int main() { NG_8 b(1,0,0,0,0,0,0,1); r2cf b1(13,11); r2cf b2(22,7); for(NG n(&b, &b1, &b2); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; for(NG n(&a, &b2, &b1); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; for(r2cf cf(286,77); cf.moreTerms(); std::cout << cf.nextTerm() << " "); std::cout << std::endl; return 0; }  Output: 3 1 2 2 3 1 2 2  [1:5,2] - [3;7] int main() { NG_8 c(0,1,-1,0,0,0,0,1); r2cf c1(13,11); r2cf c2(22,7); for(NG n(&c, &c1, &c2); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; for(r2cf cf(-151,77); cf.moreTerms(); std::cout << cf.nextTerm() << " "); std::cout << std::endl; return 0; }  Output: -1 -1 -24 -1 -2 -1 -1 -24 -1 -2  Divide [] by [3;7] int main() { NG_8 d(0,1,0,0,0,0,1,0); r2cf d1(22*22,7*7); r2cf d2(22,7); for(NG n(&d, &d1, &d2); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; return 0; }  Output: 3 7  ([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) int main() { r2cf a1(2,7); r2cf a2(13,11); NG_8 na(0,1,1,0,0,0,0,1); NG A(&na, &a1, &a2); //[0;3,2] + [1;5,2] r2cf b1(2,7); r2cf b2(13,11); NG_8 nb(0,1,-1,0,0,0,0,1); NG B(&nb, &b1, &b2); //[0;3,2] - [1;5,2] NG_8 nc(1,0,0,0,0,0,0,1); //A*B for(NG n(&nc, &A, &B); n.moreTerms(); std::cout << n.nextTerm() << " "); std::cout << std::endl; for(r2cf cf(2,7); cf.moreTerms(); std::cout << cf.nextTerm() << " "); std::cout << std::endl; for(r2cf cf(13,11); cf.moreTerms(); std::cout << cf.nextTerm() << " "); std::cout << std::endl; for(r2cf cf(-7797,5929); cf.moreTerms(); std::cout << cf.nextTerm() << " "); std::cout << std::endl; return 0; }  Common Lisp Translation of: Python (defstruct (cf (:conc-name %%cf-) (:constructor make-cf (generator))) "continued fraction" (generator nil :type function) (terminated-p nil :type boolean) (memo (make-array '(32)) :type (array integer)) (memo-count 0 :type fixnum)) (defstruct (ng8 (:constructor ng8 (a12 a1 a2 a b12 b1 b2 b))) "coefficients of a bihomographic function" (a12 0 :type integer) (a1 0 :type integer) (a2 0 :type integer) (a 0 :type integer) (b12 0 :type integer) (b1 0 :type integer) (b2 0 :type integer) (b 0 :type integer)) (defun cf-ref (cf i) "Return the ith term, or nil if there is none." (declare (cf cf) (fixnum i)) (defun get-more-terms (needed) (declare (fixnum needed)) (loop while (not (%%cf-terminated-p cf)) while (< (%%cf-memo-count cf) needed) do (let ((term (funcall (%%cf-generator cf)))) (cond (term (let ((memo (%%cf-memo cf)) (m (%%cf-memo-count cf))) (setf (aref memo m) term) (setf (%%cf-memo-count cf) (1+ m)))) (t (setf (%%cf-terminated-p cf) t)))))) (defun update (needed) (declare (fixnum needed)) (cond ((%%cf-terminated-p cf) (progn)) ((<= needed (%%cf-memo-count cf)) (progn)) ((<= needed (array-dimension (%%cf-memo cf) 0)) (get-more-terms needed)) (t (let* ((n1 (+ needed needed)) (memo1 (make-array (list n1)))) (loop for i from 0 upto (1- (%%cf-memo-count cf)) do (setf (aref memo1 i) (aref (%%cf-memo cf) i))) (setf (%%cf-memo cf) memo1) (get-more-terms needed))))) (update (1+ i)) (the (or integer null) (and (< i (%%cf-memo-count cf)) (aref (%%cf-memo cf) i)))) (defparameter *cf-max-terms* 20 "Default term-count limit for cf->string.") (defun cf->string (cf &optional (max-terms *cf-max-terms*)) "Make a readable string from a continued fraction." (declare (cf cf)) (loop with i = 0 with s = "[" do (let ((term (cf-ref cf i))) (cond ((not term) (return (concatenate 'string s "]"))) ((= i max-terms) (return (concatenate 'string s ",...]"))) (t (let ((separator (case i ((0) "") ((1) ";") (t ","))) (term-str (format nil "~A" term))) (setf i (1+ i)) (setf s (concatenate 'string s separator term-str)))))))) (defun integer->cf (num) "Transform an integer to a continued fraction." (declare (integer num)) (let ((terminated-p nil)) (declare (boolean terminated-p)) (make-cf #'(lambda () (and (not terminated-p) (progn (setf terminated-p t) num)))))) (defun ratio->cf (num) "Transform a ratio to a continued fraction." (declare (ratio num)) (let ((n (the integer (numerator num))) (d (the integer (denominator num)))) (make-cf #'(lambda () (and (not (zerop d)) (multiple-value-bind (q r) (floor n d) (setf n d) (setf d r) q)))))) ;; Thresholds chosen merely for demonstration. (defparameter number-that-is-too-big (expt 2 512)) (defparameter practically-infinite (expt 2 64)) (defun num-too-big-p (num) (declare (integer num)) (>= (abs num) (abs number-that-is-too-big))) (defun ng8-too-big-p (ng) (declare (ng8 ng)) (or (num-too-big-p (ng8-a12 ng)) (num-too-big-p (ng8-a1 ng)) (num-too-big-p (ng8-a2 ng)) (num-too-big-p (ng8-a ng)) (num-too-big-p (ng8-b12 ng)) (num-too-big-p (ng8-b1 ng)) (num-too-big-p (ng8-b2 ng)) (num-too-big-p (ng8-b ng)))) (defun treat-as-infinite-p (term) (declare (integer term)) (>= (abs term) (abs practically-infinite))) (defun quotient (u v) (declare (integer u v)) (if (zerop v) (list nil nil) (multiple-value-list (floor u v)))) (defmacro absorb-x-term (ng xsource) (let ((a12 (ng8-a12 ,ng)) (a1 (ng8-a1 ,ng)) (a2 (ng8-a2 ,ng)) (a (ng8-a ,ng)) (b12 (ng8-b12 ,ng)) (b1 (ng8-b1 ,ng)) (b2 (ng8-b2 ,ng)) (b (ng8-b ,ng)) (term (funcall ,xsource))) (if term (let ((ng^ (ng8 (+ a2 (* a12 term)) (+ a (* a1 term)) a12 a1 (+ b2 (* b12 term)) (+ b (* b1 term)) b12 b1))) (if (not (ng8-too-big-p ng^)) (setf ,ng ng^) (progn (setf ,ng (ng8 a12 a1 a12 a1 b12 b1 b12 b1)) ;; Replace the x source with one that never ;; returns a term. (setf ,xsource #'no-terms-thunk)))) (setf ,ng (ng8 a12 a1 a12 a1 b12 b1 b12 b1))))) (defmacro absorb-y-term (ng ysource) (let ((a12 (ng8-a12 ,ng)) (a1 (ng8-a1 ,ng)) (a2 (ng8-a2 ,ng)) (a (ng8-a ,ng)) (b12 (ng8-b12 ,ng)) (b1 (ng8-b1 ,ng)) (b2 (ng8-b2 ,ng)) (b (ng8-b ,ng)) (term (funcall ,ysource))) (if term (let ((ng^ (ng8 (+ a1 (* a12 term)) a12 (+ a (* a2 term)) a2 (+ b1 (* b12 term)) b12 (+ b (* b2 term)) b2))) (if (not (ng8-too-big-p ng^)) (setf ,ng ng^) (progn (setf ,ng (ng8 a12 a12 a2 a2 b12 b12 b2 b2)) ;; Replace the y source with one that never ;; returns a term. (setf ysource #'no-terms-thunk)))) (setf ,ng (ng8 a12 a12 a2 a2 b12 b12 b2 b2))))) (defun cf->thunk (cf) (let ((i 0)) #'(lambda () (let ((term (cf-ref cf i))) (setf i (1+ i)) term)))) (defun no-terms-thunk () nil) (defun apply-ng8 (ng8 x y) (declare (ng8 ng8)) (let ((ng ng8) (xsource (cf->thunk x)) (ysource (cf->thunk y))) (flet ((main () (loop with absorb for bzero = (zerop (ng8-b ng)) for b1zero = (zerop (ng8-b1 ng)) for b2zero = (zerop (ng8-b2 ng)) for b12zero = (zerop (ng8-b12 ng)) do (multiple-value-bind (q r q1 r1 q2 r2 q12 r12) (values-list (,@(quotient (ng8-a ng) (ng8-b ng)) ,@(quotient (ng8-a1 ng) (ng8-b1 ng)) ,@(quotient (ng8-a2 ng) (ng8-b2 ng)) ,@(quotient (ng8-a12 ng) (ng8-b12 ng)))) (cond ((and bzero b1zero b2zero b12zero) (return nil)) ((and bzero b2zero) (setf absorb 'x)) ((or bzero b2zero) (setf absorb 'y)) (b1zero (setf absorb 'x)) ((and (not b12zero) (= q q1 q2 q12)) ;; ;; Output a term. ;; (setf ng (ng8 (ng8-b12 ng) (ng8-b1 ng) (ng8-b2 ng) (ng8-b ng) r12 r1 r2 r)) (return (and (not (treat-as-infinite-p q)) q))) (t ;; ;; Rather than compare fractions, we will put the ;; numerators over a common denominator of ;; b*b1*b2, and then compare the new numerators. ;; (let ((n (* (ng8-a ng) (ng8-b1 ng) (ng8-b2 ng))) (n1 (* (ng8-a1 ng) (ng8-b ng) (ng8-b2 ng))) (n2 (* (ng8-a2 ng) (ng8-b ng) (ng8-b1 ng)))) (if (> (abs (- n1 n)) (abs (- n2 n))) (setf absorb 'x) (setf absorb 'y)))))) when (eq absorb 'x) do (absorb-x-term ng xsource) when (eq absorb 'y) do (absorb-y-term ng ysource)))) (make-cf #'main)))) (defun show (expression cf &optional (note "")) (format t "~A => ~A~A~%" expression (cf->string cf) note)) (defvar golden-ratio (make-cf #'(lambda () 1))) (defvar silver-ratio (make-cf #'(lambda () 2))) (defvar sqrt2 (make-cf (let ((next-term 1)) #'(lambda () (let ((term next-term)) (setf next-term 2) term))))) (defvar frac13/11 (ratio->cf 13/11)) (defvar frac22/7 (ratio->cf 22/7)) (defvar one (integer->cf 1)) (defvar two (integer->cf 2)) (defvar three (integer->cf 3)) (defvar four (integer->cf 4)) (defun cf+ (x y) (apply-ng8 (ng8 0 1 1 0 0 0 0 1) x y)) (defun cf- (x y) (apply-ng8 (ng8 0 1 -1 0 0 0 0 1) x y)) (defun cf* (x y) (apply-ng8 (ng8 1 0 0 0 0 0 0 1) x y)) (defun cf/ (x y) (apply-ng8 (ng8 0 1 0 0 0 0 1 0) x y)) (show " golden ratio" golden-ratio) (show " silver ratio" silver-ratio) (show " sqrt(2)" sqrt2) (show " 13/11" frac13/11) (show " 22/7" frac22/7) (show " 1" one) (show " 2" two) (show " 3" three) (show " 4" four) (show " (1 + 1/sqrt(2))/2" (cf/ silver-ratio (cf* sqrt2 (cf* sqrt2 sqrt2))) " method 1") (show " (1 + 1/sqrt(2))/2" (apply-ng8 (ng8 1 0 0 1 0 0 0 8) silver-ratio silver-ratio) " method 2") (show " (1 + 1/sqrt(2))/2" (cf/ (cf/ (cf/ silver-ratio sqrt2) sqrt2) sqrt2) " method 3") (show " sqrt(2) + sqrt(2)" (cf+ sqrt2 sqrt2)) (show " sqrt(2) - sqrt(2)" (cf- sqrt2 sqrt2)) (show " sqrt(2) * sqrt(2)" (cf* sqrt2 sqrt2)) (show " sqrt(2) / sqrt(2)" (cf/ sqrt2 sqrt2))  Output: $ sbcl --script bivariate-continued-fraction-task.lisp
golden ratio =>  [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]
silver ratio =>  [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2) =>  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 =>  [1;5,2]
22/7 =>  [3;7]
1 =>  [1]
2 =>  [2]
3 =>  [3]
4 =>  [4]
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]  method 1
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]  method 2
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]  method 3
sqrt(2) + sqrt(2) =>  [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) =>  [0]
sqrt(2) * sqrt(2) =>  [2]
sqrt(2) / sqrt(2) =>  [1]


D

Translation of: Python
Works with: Digital Mars D version 2.099.1
Works with: GCC version 12.2.1
//--------------------------------------------------------------------

import std.algorithm;
import std.bigint;
import std.conv;
import std.math;
import std.stdio;
import std.string;
import std.typecons;

//--------------------------------------------------------------------

class CF                        // Continued fraction.
{
alias Term = BigInt;          // The type for terms.
alias Index = size_t;         // The type for indexing terms.

protected bool terminated;    // Are there no more terms?
protected size_t m;           // The number of terms memoized.
private Term[] memo;          // Memoization storage.

static Index maxTerms = 20;   // Maximum number of terms in the
// string representation.

this ()
{
terminated = false;
m = 0;
memo.length = 32;
}

protected Nullable!Term generate ()
{
// Return terms for zero. To get different terms, override this
// method.
auto retval = (Term(0)).nullable;
if (m != 0)
retval.nullify();
return retval;
}

public Nullable!Term opIndex (Index i)
{
void update (size_t needed)
{
// Ensure all finite terms with indices 0 <= i < needed are
// memoized.
if (!terminated && m < needed)
{
if (memo.length < needed)
// To reduce the frequency of reallocation, increase the
// space to twice what might be needed right now.
memo.length = 2 * needed;

while (m != needed && !terminated)
{
auto term = generate ();
if (!term.isNull())
{
memo[m] = term.get();
m += 1;
}
else
terminated = true;
}
}
}

update (i + 1);

Nullable!Term retval;
if (i < m)
retval = memo[i].nullable;
return retval;
}

public override string toString ()
{
static string[3] separators = ["", ";", ","];

string s = "[";
Index i = 0;
bool done = false;
while (!done)
{
auto term = this[i];
if (term.isNull())
{
s ~= "]";
done = true;
}
else if (i == maxTerms)
{
s ~= ",...]";
done = true;
}
else
{
s ~= separators[(i <= 1) ? i : 2];
s ~= to!string (term.get());
i += 1;
}
}
return s;
}

public CF opBinary(string op : "+") (CF other)
{
return new cfNG8 (ng8_add, this, other);
}

public CF opBinary(string op : "-") (CF other)
{
return new cfNG8 (ng8_sub, this, other);
}

public CF opBinary(string op : "*") (CF other)
{
return new cfNG8 (ng8_mul, this, other);
}

public CF opBinary(string op : "/") (CF other)
{
return new cfNG8 (ng8_div, this, other);
}

};

//--------------------------------------------------------------------

class cfIndexed : CF  // Continued fraction with an index-to-term map.
{
alias Mapper = Nullable!Term delegate (Index);

protected Mapper map;

this (Mapper map)
{
this.map = map;
}

protected override Nullable!Term generate ()
{
return map (m);
}
}

__gshared goldenRatio =
new cfIndexed ((i) => CF.Term(1).nullable);

__gshared silverRatio =
new cfIndexed ((i) => CF.Term(2).nullable);

__gshared sqrt2 =
new cfIndexed ((i) => CF.Term(min (i + 1, 2)).nullable);

//--------------------------------------------------------------------

class cfRational : CF           // CF for a rational number.
{
private Term n;
private Term d;

this (Term numer, Term denom = Term(1))
{
n = numer;
d = denom;
}

protected override Nullable!Term generate ()
{
Nullable!Term term;
if (d != 0)
{
auto q = n / d;
auto r = n % d;
n = d;
d = r;
term = q.nullable;
}
return term;
}
}

__gshared frac_13_11 = new cfRational (CF.Term(13), CF.Term(11));
__gshared frac_22_7 = new cfRational (CF.Term(22), CF.Term(7));
__gshared one = new cfRational (CF.Term(1));
__gshared two = new cfRational (CF.Term(2));
__gshared three = new cfRational (CF.Term(3));
__gshared four = new cfRational (CF.Term(4));

//--------------------------------------------------------------------

class NG8                       // Bihomographic function.
{
public CF.Term a12, a1, a2, a;
public CF.Term b12, b1, b2, b;

this (CF.Term a12, CF.Term a1, CF.Term a2, CF.Term a,
CF.Term b12, CF.Term b1, CF.Term b2, CF.Term b)
{
this.a12 = a12;
this.a1 = a1;
this.a2 = a2;
this.a = a;
this.b12 = b12;
this.b1 = b1;
this.b2 = b2;
this.b = b;
}

this (long a12, long a1, long a2, long a,
long b12, long b1, long b2, long b)
{
this.a12 = a12;
this.a1 = a1;
this.a2 = a2;
this.a = a;
this.b12 = b12;
this.b1 = b1;
this.b2 = b2;
this.b = b;
}

this (NG8 other)
{
this.a12 = other.a12;
this.a1 = other.a1;
this.a2 = other.a2;
this.a = other.a;
this.b12 = other.b12;
this.b1 = other.b1;
this.b2 = other.b2;
this.b = other.b;
}
}

class cfNG8 : CF   // CF that is a bihomographic function of other CF.
{
private NG8 ng;
private CF x;
private CF y;
private Index ix;
private Index iy;
private bool xoverflow;
private bool yoverflow;

//
// Thresholds chosen merely for demonstration.
//
static number_that_is_too_big = // 2 ** 512
BigInt ("13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096");
static practically_infinite = // 2 ** 64
BigInt ("18446744073709551616");

this (NG8 ng, CF x, CF y)
{
this.ng = new NG8 (ng);
this.x = x;
this.y = y;
ix = 0;
iy = 0;
xoverflow = false;
yoverflow = false;
}

protected override Nullable!Term generate ()
{
// The McCabe complexity of this function is high. Please be
// careful if modifying the code.

Nullable!Term term;

bool done = false;
while (!done)
{
bool bz = (ng.b == 0);
bool b1z = (ng.b1 == 0);
bool b2z = (ng.b2 == 0);
bool b12z = (ng.b12 == 0);
if (bz && b1z && b2z && b12z)
done = true;          // There are no more terms.
else if (bz && b2z)
absorb_x_term ();
else if (bz || b2z)
absorb_y_term ();
else if (b1z)
absorb_x_term ();
else
{
Term q, r;
Term q1, r1;
Term q2, r2;
Term q12, r12;
divMod (ng.a, ng.b, q, r);
divMod (ng.a1, ng.b1, q1, r1);
divMod (ng.a2, ng.b2, q2, r2);
if (ng.b12 != 0)
divMod (ng.a12, ng.b12, q12, r12);
if (!b12z && q == q1 && q == q2 && q == q12)
{
// Output a term.
ng = new NG8 (ng.b12, ng.b1, ng.b2, ng.b,
r12, r1, r2, r);
if (!treat_as_infinite (q))
term = q.nullable;
done = true;
}
else
{
//
// Rather than compare fractions, we will put the
// numerators over a common denominator of b*b1*b2,
// and then compare the new numerators.
//
Term n = ng.a * ng.b1 * ng.b2;
Term n1 = ng.a1 * ng.b * ng.b2;
Term n2 = ng.a2 * ng.b * ng.b1;
if (abs (n1 - n) > abs (n2 - n))
absorb_x_term ();
else
absorb_y_term ();
}
}
}

return term;
}

private void absorb_x_term ()
{
Nullable!Term term;
if (!xoverflow)
term = x[ix];
ix += 1;
if (!term.isNull())
{
auto t = term.get();
auto new_ng = new NG8 (ng.a2 + (ng.a12 * t),
ng.a + (ng.a1 * t),
ng.a12, ng.a1,
ng.b2 + (ng.b12 * t),
ng.b + (ng.b1 * t),
ng.b12, ng.b1);
if (!too_big (new_ng))
ng = new_ng;
else
{
ng = new NG8 (ng.a12, ng.a1, ng.a12, ng.a1,
ng.b12, ng.b1, ng.b12, ng.b1);
xoverflow = true;
}
}
else
ng = new NG8 (ng.a12, ng.a1, ng.a12, ng.a1,
ng.b12, ng.b1, ng.b12, ng.b1);
}

private void absorb_y_term ()
{
Nullable!Term term;
if (!yoverflow)
term = y[iy];
iy += 1;
if (!term.isNull())
{
auto t = term.get();
auto new_ng = new NG8 (ng.a1 + (ng.a12 * t), ng.a12,
ng.a + (ng.a2 * t), ng.a2,
ng.b1 + (ng.b12 * t), ng.b12,
ng.b + (ng.b2 * t), ng.b2);
if (!too_big (new_ng))
ng = new_ng;
else
{
ng = new NG8 (ng.a12, ng.a12, ng.a2, ng.a2,
ng.b12, ng.b12, ng.b2, ng.b2);
yoverflow = true;
}
}
else
ng = new NG8 (ng.a12, ng.a12, ng.a2, ng.a2,
ng.b12, ng.b12, ng.b2, ng.b2);
}

private bool too_big (NG8 ng)
{
// Stop computing if a number reaches the threshold.
return (too_big (ng.a12) || too_big (ng.a1) ||
too_big (ng.a2) || too_big (ng.a) ||
too_big (ng.b12) || too_big (ng.b1) ||
too_big (ng.b2) || too_big (ng.b));
}

private bool too_big (Term u)
{
return (abs (u) >= abs (number_that_is_too_big));
}

private bool treat_as_infinite (Term u)
{
return (abs(u) >= abs (practically_infinite));
}
}

__gshared NG8 ng8_add = new NG8 (0, 1, 1, 0, 0, 0, 0, 1);
__gshared NG8 ng8_sub = new NG8 (0, 1, -1, 0, 0, 0, 0, 1);
__gshared NG8 ng8_mul = new NG8 (1, 0, 0, 0, 0, 0, 0, 1 );
__gshared NG8 ng8_div = new NG8 (0, 1, 0, 0, 0, 0, 1, 0);

//--------------------------------------------------------------------

void
show (string expression, CF cf, string note = "")
{
auto line = rightJustify (expression, 19) ~ " =>  ";
auto cf_str = to!string (cf);
if (note == "")
line ~= cf_str;
else
line ~= leftJustify (cf_str, 48) ~ note;
writeln (line);
}

int
main (char[][] args)
{
show ("golden ratio", goldenRatio, "(1 + sqrt(5))/2");
show ("silver ratio", silverRatio, "1 + sqrt(2)");
show ("sqrt(2)", sqrt2);
show ("13/11", frac_13_11);
show ("22/7", frac_22_7);
show ("one", one);
show ("two", two);
show ("three", three);
show ("four", four);
show ("(1 + 1/sqrt(2))/2",
new cfNG8 (new NG8 (0, 1, 0, 0, 0, 0, 2, 0),
silverRatio, sqrt2),
"method 1");
show ("(1 + 1/sqrt(2))/2",
new cfNG8 (new NG8 (1, 0, 0, 1, 0, 0, 0, 8),
silverRatio, silverRatio),
"method 2");
show ("(1 + 1/sqrt(2))/2", (one + (one / sqrt2)) / two,
"method 3");
show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2);
show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2);
show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2);
show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2);

return 0;
}

//--------------------------------------------------------------------

Output:
$gdc -g -Wall -Wextra bivariate_continued_fraction_task_dlang.d && ./a.out golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2 silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 1 + sqrt(2) sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 13/11 => [1;5,2] 22/7 => [3;7] one => [1] two => [2] three => [3] four => [4] (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2 (1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3 sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] sqrt(2) - sqrt(2) => [0] sqrt(2) * sqrt(2) => [2] sqrt(2) / sqrt(2) => [1]  Fortran Translation of: Python Translation of: ObjectIcon This program includes a primitive module for multiple-precision integer arithmetic. It is adequate for the task. Parts of the program might assume two's-complement representation of signed integers. The requirement that integers be two's-complement seems to me unlikely ever to become a part of Fortran standards (even though it will be required in future C standards). !--------------------------------------------------------------------- module big_integers ! Big (but not very big) integers. ! NOTE: I assume that iachar and achar do not alter the most ! significant bit. use, intrinsic :: iso_fortran_env, only: int16 implicit none private public :: big_integer public :: integer2big public :: string2big public :: big2string public :: big_sgn public :: big_cmp, big_cmpabs public :: big_neg, big_abs public :: big_addabs, big_add public :: big_subabs, big_sub public :: big_mul ! One might also include a big_muladd. public :: big_divrem ! One could also include big_div and big_rem. public :: operator(+) public :: operator(-) public :: operator(*) type :: big_integer ! The representation is sign-magnitude. The radix is 256, which ! is not speed-efficient, but which seemed relatively easy to ! work with if one were writing in standard Fortran (and assuming ! iachar and achar were "8-bit clean"). logical :: sign = .false. ! .false. => +sign, .true. => -sign. character, allocatable :: bytes(:) end type big_integer character, parameter :: zero = achar (0) character, parameter :: one = achar (1) ! An integer type capable of holding an unsigned 8-bit value. integer, parameter :: bytekind = int16 interface operator(+) module procedure big_add end interface interface operator(-) module procedure big_neg module procedure big_sub end interface interface operator(*) module procedure big_mul end interface contains elemental function logical2byte (bool) result (byte) logical, intent(in) :: bool character :: byte if (bool) then byte = one else byte = zero end if end function logical2byte elemental function logical2i (bool) result (i) logical, intent(in) :: bool integer :: i if (bool) then i = 1 else i = 0 end if end function logical2i elemental function byte2i (c) result (i) character, intent(in) :: c integer :: i i = iachar (c) end function byte2i elemental function i2byte (i) result (c) integer, intent(in) :: i character :: c c = achar (i) end function i2byte elemental function byte2bk (c) result (i) character, intent(in) :: c integer(bytekind) :: i i = iachar (c, kind = bytekind) end function byte2bk elemental function bk2byte (i) result (c) integer(bytekind), intent(in) :: i character :: c c = achar (i) end function bk2byte elemental function bk2i (i) result (j) integer(bytekind), intent(in) :: i integer :: j j = int (i) end function bk2i elemental function i2bk (i) result (j) integer, intent(in) :: i integer(bytekind) :: j j = int (iand (i, 255), kind = bytekind) end function i2bk ! Left shift of the least significant 8 bits of a bytekind integer. elemental function lshftbk (a, i) result (c) integer(bytekind), intent(in) :: a integer, intent(in) :: i integer(bytekind) :: c c = ishft (ibits (a, 0, 8 - i), i) end function lshftbk ! Right shift of the least significant 8 bits of a bytekind integer. elemental function rshftbk (a, i) result (c) integer(bytekind), intent(in) :: a integer, intent(in) :: i integer(bytekind) :: c c = ibits (a, i, 8 - i) end function rshftbk ! Left shift an integer. elemental function lshfti (a, i) result (c) integer, intent(in) :: a integer, intent(in) :: i integer :: c c = ishft (a, i) end function lshfti ! Right shift an integer. elemental function rshfti (a, i) result (c) integer, intent(in) :: a integer, intent(in) :: i integer :: c c = ishft (a, -i) end function rshfti function integer2big (i) result (a) integer, intent(in) :: i type(big_integer), allocatable :: a ! ! To write a more efficient implementation of this procedure is ! left as an exercise for the reader. ! character(len = 100) :: buffer write (buffer, '(I0)') i a = string2big (trim (buffer)) end function integer2big function string2big (s) result (a) character(len = *), intent(in) :: s type(big_integer), allocatable :: a integer :: n, i, istart, iend integer :: digit if ((s(1:1) == '-') .or. s(1:1) == '+') then istart = 2 else istart = 1 end if iend = len (s) n = (iend - istart + 2) / 2 allocate (a) allocate (a%bytes(n)) a%bytes = zero do i = istart, iend digit = ichar (s(i:i)) - ichar ('0') if (digit < 0 .or. 9 < digit) error stop a = short_multiplication (a, 10) a = short_addition (a, digit) end do a%sign = (s(1:1) == '-') call normalize (a) end function string2big function big2string (a) result (s) type(big_integer), intent(in) :: a character(len = :), allocatable :: s type(big_integer), allocatable :: q integer :: r integer :: sgn sgn = big_sgn (a) if (sgn == 0) then s = '0' else q = a s = '' do while (big_sgn (q) /= 0) call short_division (q, 10, q, r) s = achar (r + ichar ('0')) // s end do if (sgn < 0) s = '-' // s end if end function big2string function big_sgn (a) result (sgn) type(big_integer), intent(in) :: a integer :: sgn integer :: n, i n = size (a%bytes) i = 1 sgn = 1234 do while (sgn == 1234) if (i == n + 1) then sgn = 0 else if (a%bytes(i) /= zero) then if (a%sign) then sgn = -1 else sgn = 1 end if else i = i + 1 end if end do end function big_sgn function big_cmp (a, b) result (cmp) type(big_integer(*)), intent(in) :: a, b integer :: cmp if (a%sign) then if (b%sign) then cmp = -big_cmpabs (a, b) else cmp = -1 end if else if (b%sign) then cmp = 1 else cmp = big_cmpabs (a, b) end if end if end function big_cmp function big_cmpabs (a, b) result (cmp) type(big_integer(*)), intent(in) :: a, b integer :: cmp integer :: n, i integer :: ia, ib cmp = 1234 n = max (size (a%bytes), size (b%bytes)) i = n do while (cmp == 1234) if (i == 0) then cmp = 0 else ia = byteval (a, i) ib = byteval (b, i) if (ia < ib) then cmp = -1 else if (ia > ib) then cmp = 1 else i = i - 1 end if end if end do end function big_cmpabs function big_neg (a) result (c) type(big_integer), intent(in) :: a type(big_integer), allocatable :: c c = a c%sign = .not. c%sign end function big_neg function big_abs (a) result (c) type(big_integer), intent(in) :: a type(big_integer), allocatable :: c c = a c%sign = .false. end function big_abs function big_add (a, b) result (c) type(big_integer), intent(in) :: a type(big_integer), intent(in) :: b type(big_integer), allocatable :: c logical :: sign if (a%sign) then if (b%sign) then ! a <= 0, b <= 0 c = big_addabs (a, b) sign = .true. else ! a <= 0, b >= 0 c = big_subabs (a, b) sign = .not. c%sign end if else if (b%sign) then ! a >= 0, b <= 0 c = big_subabs (a, b) sign = c%sign else ! a >= 0, b >= 0 c = big_addabs (a, b) sign = .false. end if end if c%sign = sign end function big_add function big_sub (a, b) result (c) type(big_integer), intent(in) :: a type(big_integer), intent(in) :: b type(big_integer), allocatable :: c logical :: sign if (a%sign) then if (b%sign) then ! a <= 0, b <= 0 c = big_subabs (a, b) sign = .not. c%sign else ! a <= 0, b >= 0 c = big_addabs (a, b) sign = .true. end if else if (b%sign) then ! a >= 0, b <= 0 c = big_addabs (a, b) sign = .false. else ! a >= 0, b >= 0 c = big_subabs (a, b) sign = c%sign end if end if c%sign = sign end function big_sub function big_addabs (a, b) result (c) type(big_integer), intent(in) :: a, b type(big_integer), allocatable :: c ! Compute abs(a) + abs(b). integer :: n, nc, i logical :: carry type(big_integer), allocatable :: tmp n = max (size (a%bytes), size (b%bytes)) nc = n + 1 allocate(tmp) allocate(tmp%bytes(nc)) call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry) do i = 2, n call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry) end do tmp%bytes(nc) = logical2byte (carry) call normalize (tmp) c = tmp end function big_addabs function big_subabs (a, b) result (c) type(big_integer), intent(in) :: a, b type(big_integer), allocatable :: c ! Compute abs(a) - abs(b). The result is signed. integer :: n, i logical :: carry type(big_integer), allocatable :: tmp n = max (size (a%bytes), size (b%bytes)) allocate(tmp) allocate(tmp%bytes(n)) if (big_cmpabs (a, b) >= 0) then tmp%sign = .false. call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry) do i = 2, n call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry) end do else tmp%sign = .true. call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry) do i = 2, n call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry) end do end if call normalize (tmp) c = tmp end function big_subabs function big_mul (a, b) result (c) type(big_integer), intent(in) :: a, b type(big_integer), allocatable :: c ! ! This is Knuth, Volume 2, Algorithm 4.3.1M. ! integer :: na, nb, nc integer :: i, j integer :: ia, ib, ic integer :: carry type(big_integer), allocatable :: tmp na = size (a%bytes) nb = size (b%bytes) nc = na + nb + 1 allocate (tmp) allocate (tmp%bytes(nc)) tmp%bytes = zero j = 1 do j = 1, nb ib = byte2i (b%bytes(j)) if (ib /= 0) then carry = 0 do i = 1, na ia = byte2i (a%bytes(i)) ic = byte2i (tmp%bytes(i + j - 1)) ic = (ia * ib) + ic + carry tmp%bytes(i + j - 1) = i2byte (iand (ic, 255)) carry = ishft (ic, -8) end do tmp%bytes(na + j) = i2byte (carry) end if end do tmp%sign = (a%sign .neqv. b%sign) call normalize (tmp) c = tmp end function big_mul subroutine big_divrem (a, b, q, r) type(big_integer), intent(in) :: a, b type(big_integer), allocatable, intent(inout) :: q, r ! ! Division with a remainder that is never negative. Equivalently, ! this is floor division if the divisor is positive, and ceiling ! division if the divisor is negative. ! ! See Raymond T. Boute, "The Euclidean definition of the functions ! div and mod", ACM Transactions on Programming Languages and ! Systems, Volume 14, Issue 2, pp. 127-144. ! https://doi.org/10.1145/128861.128862 ! call nonnegative_division (a, b, .true., .true., q, r) if (a%sign) then if (big_sgn (r) /= 0) then q = short_addition (q, 1) r = big_sub (big_abs (b), r) end if q%sign = .not. b%sign else q%sign = b%sign end if end subroutine big_divrem function short_addition (a, b) result (c) type(big_integer), intent(in) :: a integer, intent(in) :: b type(big_integer), allocatable :: c ! Compute abs(a) + b. integer :: na, nc, i logical :: carry type(big_integer), allocatable :: tmp na = size (a%bytes) nc = na + 1 allocate(tmp) allocate(tmp%bytes(nc)) call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry) do i = 2, na call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry) end do tmp%bytes(nc) = logical2byte (carry) call normalize (tmp) c = tmp end function short_addition function short_multiplication (a, b) result (c) type(big_integer), intent(in) :: a integer, intent(in) :: b type(big_integer), allocatable :: c integer :: i, na, nc integer :: ia, ic integer :: carry type(big_integer), allocatable :: tmp na = size (a%bytes) nc = na + 1 allocate (tmp) allocate (tmp%bytes(nc)) tmp%sign = a%sign carry = 0 do i = 1, na ia = byte2i (a%bytes(i)) ic = (ia * b) + carry tmp%bytes(i) = i2byte (iand (ic, 255)) carry = ishft (ic, -8) end do tmp%bytes(nc) = i2byte (carry) call normalize (tmp) c = tmp end function short_multiplication ! Division without regard to signs. subroutine nonnegative_division (a, b, want_q, want_r, q, r) type(big_integer), intent(in) :: a, b logical, intent(in) :: want_q, want_r type(big_integer), intent(inout), allocatable :: q, r integer :: na, nb integer :: remainder na = size (a%bytes) nb = size (b%bytes) ! It is an error if b has "significant" zero-bytes or is equal to ! zero. if (b%bytes(nb) == zero) error stop if (nb == 1) then if (want_q) then call short_division (a, byte2i (b%bytes(1)), q, remainder) else block type(big_integer), allocatable :: bit_bucket call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder) end block end if if (want_r) then if (allocated (r)) deallocate (r) allocate (r) allocate (r%bytes(1)) r%bytes(1) = i2byte (remainder) end if else if (na >= nb) then call long_division (a, b, want_q, want_r, q, r) else if (want_q) q = string2big ("0") if (want_r) r = a end if end if end subroutine nonnegative_division subroutine short_division (a, b, q, r) type(big_integer), intent(in) :: a integer, intent(in) :: b type(big_integer), intent(inout), allocatable :: q integer, intent(inout) :: r ! ! This is Knuth, Volume 2, Exercise 4.3.1.16. ! ! The divisor is assumed to be positive. ! integer :: n, i integer :: ia, ib, iq type(big_integer), allocatable :: tmp ib = b n = size (a%bytes) allocate (tmp) allocate (tmp%bytes(n)) r = 0 do i = n, 1, -1 ia = (256 * r) + byte2i (a%bytes(i)) iq = ia / ib r = mod (ia, ib) tmp%bytes(i) = i2byte (iq) end do tmp%sign = a%sign call normalize (tmp) q = tmp end subroutine short_division subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder) type(big_integer), intent(in) :: a, b logical, intent(in) :: want_quotient, want_remainder type(big_integer), intent(inout), allocatable :: quotient type(big_integer), intent(inout), allocatable :: remainder ! ! This is Knuth, Volume 2, Algorithm 4.3.1D. ! ! We do not deal here with the signs of the inputs and outputs. ! ! It is assumed size(a%bytes) >= size(b%bytes), and that b has no ! leading zero-bytes and is at least two bytes long. If b is one ! byte long and nonzero, use short division. ! integer :: na, nb, m, n integer :: num_lz, num_nonlz integer :: j integer :: qhat logical :: carry ! ! We will NOT be working with VERY large numbers, and so it will ! be safe to put temporary storage on the stack. (Note: your ! Fortran might put this storage in a heap instead of the stack.) ! ! v = b, normalized to put its most significant 1-bit all the ! way left. ! ! u = a, shifted left by the same amount as b. ! ! q = the quotient. ! ! The remainder, although shifted left, will end up in u. ! integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes)) integer(bytekind) :: v(0:size (b%bytes) - 1) integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes)) na = size (a%bytes) nb = size (b%bytes) n = nb m = na - nb ! In the most significant byte of the divisor, find the number of ! leading zero bits, and the number of bits after that. block integer(bytekind) :: tmp tmp = byte2bk (b%bytes(n)) num_nonlz = bit_size (tmp) - leadz (tmp) num_lz = 8 - num_nonlz end block call normalize_v (b%bytes) ! Make the most significant bit of v be one. call normalize_u (a%bytes) ! Shifted by the same amount as v. ! Assure ourselves that the most significant bit of v is a one. if (.not. btest (v(n - 1), 7)) error stop do j = m, 0, -1 call calculate_qhat (qhat) call multiply_and_subtract (carry) q(j) = i2bk (qhat) if (carry) call add_back end do if (want_quotient) then if (allocated (quotient)) deallocate (quotient) allocate (quotient) allocate (quotient%bytes(m + 1)) quotient%bytes = bk2byte (q) call normalize (quotient) end if if (want_remainder) then if (allocated (remainder)) deallocate (remainder) allocate (remainder) allocate (remainder%bytes(n)) call unnormalize_u (remainder%bytes) call normalize (remainder) end if contains subroutine normalize_v (b_bytes) character, intent(in) :: b_bytes(n) ! ! Normalize v so its most significant bit is a one. Any ! normalization factor that achieves this goal will suffice; we ! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).) ! ! Strictly for readability, we use linear stack space for an ! intermediate result. ! integer :: i integer(bytekind) :: btmp(0:n - 1) btmp = byte2bk (b_bytes) v(0) = lshftbk (btmp(0), num_lz) do i = 1, n - 1 v(i) = ior (lshftbk (btmp(i), num_lz), & & rshftbk (btmp(i - 1), num_nonlz)) end do end subroutine normalize_v subroutine normalize_u (a_bytes) character, intent(in) :: a_bytes(m + n) ! ! Shift a leftwards to get u. Shift by as much as b was shifted ! to get v. ! ! Strictly for readability, we use linear stack space for an ! intermediate result. ! integer :: i integer(bytekind) :: atmp(0:m + n - 1) atmp = byte2bk (a_bytes) u(0) = lshftbk (atmp(0), num_lz) do i = 1, m + n - 1 u(i) = ior (lshftbk (atmp(i), num_lz), & & rshftbk (atmp(i - 1), num_nonlz)) end do u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz) end subroutine normalize_u subroutine unnormalize_u (r_bytes) character, intent(out) :: r_bytes(n) ! ! Strictly for readability, we use linear stack space for an ! intermediate result. ! integer :: i integer(bytekind) :: rtmp(0:n - 1) do i = 0, n - 1 rtmp(i) = ior (rshftbk (u(i), num_lz), & & lshftbk (u(i + 1), num_nonlz)) end do rtmp(n - 1) = rshftbk (u(n - 1), num_lz) r_bytes = bk2byte (rtmp) end subroutine unnormalize_u subroutine calculate_qhat (qhat) integer, intent(out) :: qhat integer :: itmp, rhat logical :: adjust itmp = ior (lshfti (bk2i (u(j + n)), 8), & & bk2i (u(j + n - 1))) qhat = itmp / bk2i (v(n - 1)) rhat = mod (itmp, bk2i (v(n - 1))) adjust = .true. do while (adjust) if (rshfti (qhat, 8) /= 0) then continue else if (qhat * bk2i (v(n - 2)) & & > ior (lshfti (rhat, 8), & & bk2i (u(j + n - 2)))) then continue else adjust = .false. end if if (adjust) then qhat = qhat - 1 rhat = rhat + bk2i (v(n - 1)) if (rshfti (rhat, 8) == 0) then adjust = .false. end if end if end do end subroutine calculate_qhat subroutine multiply_and_subtract (carry) logical, intent(out) :: carry integer :: i integer :: qhat_v integer :: mul_carry, sub_carry integer :: diff mul_carry = 0 sub_carry = 0 do i = 0, n ! Multiplication. qhat_v = mul_carry if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i))) mul_carry = rshfti (qhat_v, 8) qhat_v = iand (qhat_v, 255) ! Subtraction. diff = bk2i (u(j + i)) - qhat_v + sub_carry sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1. u(j + i) = i2bk (diff) end do carry = (sub_carry /= 0) end subroutine multiply_and_subtract subroutine add_back integer :: i, carry, sum q(j) = q(j) - 1_bytekind carry = 0 do i = 0, n - 1 sum = bk2i (u(j + i)) + bk2i (v(i)) + carry carry = ishft (sum, -8) u(j + i) = i2bk (sum) end do end subroutine add_back end subroutine long_division subroutine add_bytes (a, b, carry_in, c, carry_out) character, intent(in) :: a, b logical, value :: carry_in character, intent(inout) :: c logical, intent(inout) :: carry_out integer :: ia, ib, ic ia = byte2i (a) if (carry_in) ia = ia + 1 ib = byte2i (b) ic = ia + ib c = i2byte (iand (ic, 255)) carry_out = (ic >= 256) end subroutine add_bytes subroutine sub_bytes (a, b, carry_in, c, carry_out) character, intent(in) :: a, b logical, value :: carry_in character, intent(inout) :: c logical, intent(inout) :: carry_out integer :: ia, ib, ic ia = byte2i (a) ib = byte2i (b) if (carry_in) ib = ib + 1 ic = ia - ib carry_out = (ic < 0) if (carry_out) ic = ic + 256 c = i2byte (iand (ic, 255)) end subroutine sub_bytes function get_byte (a, i) result (byte) type(big_integer), intent(in) :: a integer, intent(in) :: i character :: byte if (size (a%bytes) < i) then byte = zero else byte = a%bytes(i) end if end function get_byte function byteval (a, i) result (v) type(big_integer), intent(in) :: a integer, intent(in) :: i integer :: v if (size (a%bytes) < i) then v = 0 else v = byte2i (a%bytes(i)) end if end function byteval subroutine normalize (a) type(big_integer), intent(inout) :: a logical :: done integer :: i character, allocatable :: fewer_bytes(:) ! Shorten to the minimum number of bytes. i = size (a%bytes) done = .false. do while (.not. done) if (i == 1) then done = .true. else if (a%bytes(i) /= zero) then done = .true. else i = i - 1 end if end do if (i /= size (a%bytes)) then allocate (fewer_bytes (i)) fewer_bytes = a%bytes(1:i) call move_alloc (fewer_bytes, a%bytes) end if ! If the magnitude is zero, then clear the sign bit. if (size (a%bytes) == 1) then if (a%bytes(1) == zero) then a%sign = .false. end if end if end subroutine normalize end module big_integers !--------------------------------------------------------------------- module continued_fractions use, non_intrinsic :: big_integers implicit none private public :: continued_fraction public :: term_generator public :: term_generator_procedure public :: make_continued_fraction public :: i2cf public :: make_integer_continued_fraction public :: make_integer_continued_fraction_from_integer public :: constant_term_cf public :: make_constant_term_continued_fraction public :: make_constant_term_continued_fraction_from_integer public :: apply_ng8 public :: apply_ng8_big_integers public :: apply_ng8_integers public :: ng8_coefficient_threshold public :: ng8_term_threshold public :: add_continued_fractions public :: subtract_continued_fractions public :: multiply_continued_fractions public :: divide_continued_fractions public :: cf2string public :: continued_fraction_to_string_given_max_terms public :: continued_fraction_to_string_with_default_max_terms public :: default_continued_fraction_max_terms type :: continued_fraction class(continued_fraction_record), pointer, private :: p => null () contains procedure, pass :: get_term => get_continued_fraction_term procedure, pass :: term_exists => continued_fraction_term_exists procedure, pass :: term => continued_fraction_term procedure, pass :: to_string => continued_fraction_to_string_with_default_max_terms procedure, pass :: add => add_continued_fractions generic :: operator(+) => add procedure, pass :: subtract => subtract_continued_fractions generic :: operator(-) => subtract procedure, pass :: multiply => multiply_continued_fractions generic :: operator(*) => multiply procedure, pass :: divide => divide_continued_fractions generic :: operator(/) => divide procedure, pass, private :: continued_fraction_make_new_ref generic :: assignment(=) => continued_fraction_make_new_ref final :: continued_fraction_final end type continued_fraction type :: continued_fraction_record logical, private :: terminated = .false. ! No more terms? integer, private :: m = 0 ! No. of terms memoized. type(big_integer), private, allocatable :: memo(:) ! Memoized terms. class(term_generator), pointer :: gen ! Where terms come from. integer :: refcount = 0 contains procedure, pass :: get_term => get_continued_fraction_record_term procedure, pass :: term_exists => continued_fraction_record_term_exists procedure, pass :: term => continued_fraction_record_term final :: continued_fraction_record_final end type continued_fraction_record type, abstract :: term_generator contains procedure(term_generator_procedure), pass, deferred :: generate end type term_generator interface subroutine term_generator_procedure (gen, term_exists, term) import term_generator import big_integer class(term_generator), intent(inout) :: gen logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term end subroutine term_generator_procedure end interface type, extends (term_generator) :: integer_term_generator type(big_integer), allocatable :: term logical :: no_more_terms = .false. contains procedure, pass :: generate => integer_term_generator_generate end type integer_term_generator type, extends (term_generator) :: constant_term_generator type(big_integer), allocatable :: term contains procedure, pass :: generate => constant_term_generator_generate end type constant_term_generator type, extends (term_generator) :: ng8_term_generator type(big_integer), allocatable :: a12, a1, a2, a type(big_integer), allocatable :: b12, b1, b2, b type(continued_fraction) :: x, y integer :: ix = 0 integer :: iy = 0 logical :: x_overflow = .false. logical :: y_overflow = .false. contains procedure, pass :: generate => ng8_term_generator_generate end type ng8_term_generator interface i2cf module procedure make_integer_continued_fraction module procedure make_integer_continued_fraction_from_integer end interface i2cf interface constant_term_cf module procedure make_constant_term_continued_fraction module procedure make_constant_term_continued_fraction_from_integer end interface constant_term_cf interface apply_ng8 module procedure apply_ng8_big_integers module procedure apply_ng8_integers end interface apply_ng8 interface cf2string module procedure continued_fraction_to_string_given_max_terms module procedure continued_fraction_to_string_with_default_max_terms end interface cf2string integer :: default_continued_fraction_max_terms = 20 type(big_integer), allocatable :: ng8_coefficient_threshold type(big_integer), allocatable :: ng8_term_threshold contains subroutine continued_fraction_make_new_ref (dst, src) class(continued_fraction), intent(inout) :: dst class(continued_fraction), intent(in) :: src if (associated (dst%p)) deallocate (dst%p) dst%p => src%p dst%p%refcount = dst%p%refcount + 1 end subroutine continued_fraction_make_new_ref subroutine continued_fraction_final (cf) type(continued_fraction), intent(inout) :: cf cf%p%refcount = cf%p%refcount - 1 if (cf%p%refcount == 0) deallocate (cf%p) end subroutine continued_fraction_final function make_continued_fraction (gen) result (cf) class(term_generator), pointer, intent(in) :: gen type(continued_fraction) :: cf allocate (cf%p) allocate (cf%p%memo(0:31)) ! The starting size is arbitrary. cf%p%gen => gen cf%p%refcount = cf%p%refcount + 1 end function make_continued_fraction subroutine continued_fraction_record_final (cfrec) type(continued_fraction_record), intent(inout) :: cfrec deallocate (cfrec%gen) end subroutine continued_fraction_record_final function make_integer_continued_fraction (bigint) result (cf) type(big_integer), intent(in) :: bigint type(continued_fraction) :: cf class(integer_term_generator), pointer :: gen allocate (gen) gen%term = bigint cf = make_continued_fraction (gen) end function make_integer_continued_fraction function make_integer_continued_fraction_from_integer (i) result (cf) integer, intent(in) :: i type(continued_fraction) :: cf cf = make_integer_continued_fraction (integer2big (i)) end function make_integer_continued_fraction_from_integer subroutine integer_term_generator_generate (gen, term_exists, term) class(integer_term_generator), intent(inout) :: gen logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term term_exists = (.not. gen%no_more_terms) if (term_exists) term = gen%term gen%no_more_terms = .true. end subroutine integer_term_generator_generate function make_constant_term_continued_fraction (bigint) result (cf) type(big_integer), intent(in) :: bigint type(continued_fraction) :: cf class(constant_term_generator), pointer :: gen allocate (gen) gen%term = bigint cf = make_continued_fraction (gen) end function make_constant_term_continued_fraction function make_constant_term_continued_fraction_from_integer (i) result (cf) integer, intent(in) :: i type(continued_fraction) :: cf cf = make_constant_term_continued_fraction (integer2big (i)) end function make_constant_term_continued_fraction_from_integer subroutine constant_term_generator_generate (gen, term_exists, term) class(constant_term_generator), intent(inout) :: gen logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term term_exists = .true. if (term_exists) term = gen%term end subroutine constant_term_generator_generate function apply_ng8_big_integers (a12, a1, a2, a, & & b12, b1, b2, b, x, y) result (cf) type(big_integer), intent(in) :: a12, a1, a2, a type(big_integer), intent(in) :: b12, b1, b2, b class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf class(ng8_term_generator), pointer :: gen allocate (gen) gen%a12 = a12; gen%a1 = a1; gen%a2 = a2; gen%a = a gen%b12 = b12; gen%b1 = b1; gen%b2 = b2; gen%b = b gen%x = x gen%y = y cf = make_continued_fraction (gen) end function apply_ng8_big_integers function apply_ng8_integers (a12, a1, a2, a, & & b12, b1, b2, b, x, y) result (cf) integer, intent(in) :: a12, a1, a2, a integer, intent(in) :: b12, b1, b2, b class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf cf = apply_ng8_big_integers (integer2big (a12), & & integer2big (a1), & & integer2big (a2), & & integer2big (a), & & integer2big (b12), & & integer2big (b1), & & integer2big (b2), & & integer2big (b), x, y) end function apply_ng8_integers subroutine ng8_term_generator_generate (gen, term_exists, term) class(ng8_term_generator), intent(inout) :: gen logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term logical :: done logical :: b12z, b1z, b2z, bz type(big_integer), allocatable :: q12, r12 type(big_integer), allocatable :: q1, r1 type(big_integer), allocatable :: q2, r2 type(big_integer), allocatable :: q, r logical :: absorb_x, absorb_y, compare_fracs done = .false. do while (.not. done) absorb_x = .false. absorb_y = .false. compare_fracs = .false. b12z = (big_sgn (gen%b12) == 0) b1z = (big_sgn (gen%b1) == 0) b2z = (big_sgn (gen%b2) == 0) bz = (big_sgn (gen%b) == 0) if (b12z .and. b1z .and. b2z .and. bz) then ! There are no more terms. term_exists = .false. done = .true. else if (b2z .and. bz) then absorb_x = .true. else if (b2z .or. bz) then absorb_y = .true. else if (b1z) then absorb_x = .true. else call big_divrem (gen%a1, gen%b1, q1, r1) call big_divrem (gen%a2, gen%b2, q2, r2) call big_divrem (gen%a, gen%b, q, r) if (.not. b12z) then call big_divrem (gen%a12, gen%b12, q12, r12) if (big_cmp (q, q1) /= 0) then compare_fracs = .true. else if (big_cmp (q, q2) /= 0) then compare_fracs = .true. else if (big_cmp (q, q12) /= 0) then compare_fracs = .true. else call output_term done = .true. end if end if end if if (compare_fracs) call compare_fractions (absorb_x, absorb_y) if (absorb_x) call absorb_x_term if (absorb_y) call absorb_y_term end do contains subroutine output_term gen%a12 = gen%b12; gen%a1 = gen%b1; gen%a2 = gen%b2; gen%a = gen%b gen%b12 = r12; gen%b1 = r1; gen%b2 = r2; gen%b = r term_exists = (.not. treat_as_infinite (q)) if (term_exists) term = q end subroutine output_term subroutine compare_fractions (absorb_x, absorb_y) logical, intent(inout) :: absorb_x, absorb_y ! Rather than compare fractions, we will put the numerators over ! a common denominator of b1*b2*b, and then compare the new ! numerators. type(big_integer), allocatable :: n1, n2, n n1 = gen%a1 * gen%b2 * gen%b n2 = gen%a2 * gen%b1 * gen%b n = gen%a * gen%b1 * gen%b2 if (big_cmpabs (n1 - n, n2 - n) > 0) then absorb_x = .true. else absorb_y = .true. end if end subroutine compare_fractions subroutine absorb_x_term logical :: term_exists type(big_integer), allocatable :: term type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b if (gen%x_overflow) then term_exists = .false. else term_exists = gen%x%term_exists(gen%ix) end if new_a2 = gen%a12 new_a = gen%a1 new_b2 = gen%b12 new_b = gen%b1 if (term_exists) then term = gen%x%term(gen%ix) new_a12 = gen%a2 + (gen%a12 * term) new_a1 = gen%a + (gen%a1 * term) new_b12 = gen%b2 + (gen%b12 * term) new_b1 = gen%b + (gen%b1 * term) if (any_too_big (new_a12, new_a1, new_a2, new_a, & & new_b12, new_b1, new_b2, new_b)) then gen%x_overflow = .true. new_a12 = gen%a12 new_a1 = gen%a1 new_b12 = gen%b12 new_b1 = gen%b1 end if else new_a12 = gen%a12 new_a1 = gen%a1 new_b12 = gen%b12 new_b1 = gen%b1 end if gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b gen%ix = gen%ix + 1 end subroutine absorb_x_term subroutine absorb_y_term logical :: term_exists type(big_integer), allocatable :: term type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b if (gen%y_overflow) then term_exists = .false. else term_exists = gen%y%term_exists(gen%iy) end if new_a1 = gen%a12 new_a = gen%a2 new_b1 = gen%b12 new_b = gen%b2 if (term_exists) then term = gen%y%term(gen%iy) new_a12 = gen%a1 + (gen%a12 * term) new_a2 = gen%a + (gen%a2 * term) new_b12 = gen%b1 + (gen%b12 * term) new_b2 = gen%b + (gen%b2 * term) if (any_too_big (new_a12, new_a1, new_a2, new_a, & & new_b12, new_b1, new_b2, new_b)) then gen%y_overflow = .true. new_a12 = gen%a12 new_a2 = gen%a2 new_b12 = gen%b12 new_b2 = gen%b2 end if else new_a12 = gen%a12 new_a2 = gen%a2 new_b12 = gen%b12 new_b2 = gen%b2 end if gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b gen%iy = gen%iy + 1 end subroutine absorb_y_term function any_too_big (a, b, c, d, e, f, g, h) result (yes) type(big_integer), intent(in) :: a, b, c, d, e, f, g, h logical :: yes if (too_big (a)) then yes = .true. else if (too_big (b)) then yes = .true. else if (too_big (c)) then yes = .true. else if (too_big (d)) then yes = .true. else if (too_big (e)) then yes = .true. else if (too_big (f)) then yes = .true. else if (too_big (g)) then yes = .true. else if (too_big (h)) then yes = .true. else yes = .false. end if end function any_too_big function too_big (coef) result (yes) type(big_integer), intent(in) :: coef logical :: yes if (.not. allocated (ng8_coefficient_threshold)) then ! 2**512 ng8_coefficient_threshold = string2big ('1340780792994259709957& &402499820584612747936582059239337772356144372176403007354& &697680187429816690342769003185818648605085375388281194656& &9946433649006084096') end if yes = (big_cmpabs (coef, ng8_coefficient_threshold) >= 0) end function too_big function treat_as_infinite (term) result (yes) type(big_integer), intent(in) :: term logical :: yes if (.not. allocated (ng8_term_threshold)) then ! 2**64 ng8_term_threshold = string2big ('18446744073709551616') end if yes = (big_cmpabs (term, ng8_term_threshold) >= 0) end function treat_as_infinite end subroutine ng8_term_generator_generate function add_continued_fractions (x, y) result (cf) class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf cf = apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1, x, y) end function add_continued_fractions function subtract_continued_fractions (x, y) result (cf) class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf cf = apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1, x, y) end function subtract_continued_fractions function multiply_continued_fractions (x, y) result (cf) class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf cf = apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1, x, y) end function multiply_continued_fractions function divide_continued_fractions (x, y) result (cf) class(continued_fraction), intent(in) :: x, y type(continued_fraction) :: cf cf = apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0, x, y) end function divide_continued_fractions subroutine get_continued_fraction_term (cf, i, term_exists, term) class(continued_fraction), intent(in) :: cf integer, intent(in) :: i logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term call get_continued_fraction_record_term (cf%p, i, term_exists, term) end subroutine get_continued_fraction_term subroutine get_continued_fraction_record_term (cfrec, i, term_exists, term) class(continued_fraction_record), intent(inout) :: cfrec integer, intent(in) :: i logical, intent(out) :: term_exists type(big_integer), allocatable, intent(out) :: term if (i < 0) error stop call update (i + 1) term_exists = (i < cfrec%m) if (term_exists) term = cfrec%memo(i) contains subroutine update (needed) integer :: needed logical :: term_exists1 type(big_integer), allocatable :: term1 if (.not. cfrec%terminated .and. cfrec%m < needed) then if (size (cfrec%memo) < needed) then block ! Allocate more storage. type(big_integer), allocatable :: memo1(:) allocate (memo1(0 : (2 * needed) - 1)) memo1(0:cfrec%m - 1) = cfrec%memo(0:cfrec%m - 1) call move_alloc (memo1, cfrec%memo) end block end if do while (.not. cfrec%terminated .and. cfrec%m < needed) call cfrec%gen%generate (term_exists1, term1) if (term_exists1) then cfrec%memo(cfrec%m) = term1 cfrec%m = cfrec%m + 1 else cfrec%terminated = .true. end if end do end if end subroutine update end subroutine get_continued_fraction_record_term function continued_fraction_term_exists (cf, i) result (term_exists) class(continued_fraction), intent(in) :: cf integer, intent(in) :: i logical :: term_exists term_exists = continued_fraction_record_term_exists (cf%p, i) end function continued_fraction_term_exists function continued_fraction_record_term_exists (cfrec, i) result (term_exists) class(continued_fraction_record), intent(inout) :: cfrec integer, intent(in) :: i logical :: term_exists type(big_integer), allocatable :: term call get_continued_fraction_record_term (cfrec, i, term_exists, term) end function continued_fraction_record_term_exists function continued_fraction_term (cf, i) result (term) class(continued_fraction), intent(in) :: cf integer, intent(in) :: i type(big_integer), allocatable :: term term = continued_fraction_record_term (cf%p, i) end function continued_fraction_term function continued_fraction_record_term (cfrec, i) result (term) class(continued_fraction_record), intent(inout) :: cfrec integer, intent(in) :: i type(big_integer), allocatable :: term logical :: term_exists call get_continued_fraction_record_term (cfrec, i, term_exists, term) if (.not. term_exists) error stop end function continued_fraction_record_term function continued_fraction_to_string_given_max_terms (cf, max_terms) result (s) class(continued_fraction), intent(in) :: cf integer, intent(in) :: max_terms character(len = :), allocatable :: s s = continued_fraction_record_to_string_given_max_terms (cf%p, max_terms) end function continued_fraction_to_string_given_max_terms function continued_fraction_record_to_string_given_max_terms (cfrec, max_terms) result (s) class(continued_fraction_record), intent(inout) :: cfrec integer, intent(in) :: max_terms character(len = :), allocatable :: s integer :: i logical :: done i = 0 s = '[' done = .false. do while (.not. done) if (.not. cfrec%term_exists(i)) then s = s // "]" done = .true. else if (i == max_terms) then s = s // ",...]" done = .true. else select case (i) case (0); continue case (1); s = s // ";" case default; s = s // "," end select s = s // big2string (cfrec%term(i)) i = i + 1 end if end do end function continued_fraction_record_to_string_given_max_terms function continued_fraction_to_string_with_default_max_terms (cf) result (s) class(continued_fraction), intent(in) :: cf character(len = :), allocatable :: s s = continued_fraction_record_to_string_with_default_max_terms (cf%p) end function continued_fraction_to_string_with_default_max_terms function continued_fraction_record_to_string_with_default_max_terms (cfrec) result (s) class(continued_fraction_record), intent(inout) :: cfrec character(len = :), allocatable :: s integer :: max_terms max_terms = max (default_continued_fraction_max_terms, 1) s = continued_fraction_record_to_string_given_max_terms (cfrec, max_terms) end function continued_fraction_record_to_string_with_default_max_terms end module continued_fractions !--------------------------------------------------------------------- program bivariate_continued_fraction_task use, non_intrinsic :: big_integers use, non_intrinsic :: continued_fractions implicit none type(continued_fraction) :: golden_ratio type(continued_fraction) :: silver_ratio type(continued_fraction) :: sqrt2 type(continued_fraction) :: one type(continued_fraction) :: two type(continued_fraction) :: three type(continued_fraction) :: four type(continued_fraction) :: method1 type(continued_fraction) :: method2 type(continued_fraction) :: method3 block ! ! Let us start with a test of the long division routine, with some ! numbers known to trigger a bug in the first version I ! posted. That version had a buggy "add_back" routine. ! ! (How I found such numbers is easy: I used random search.) ! type(big_integer), allocatable :: a, b, q, r a = string2big ("95292350834616415707142739736356097545484658215015733475& &690528634954280668802285176284181482202905629004984123564335019024318905") b = string2big ("63677949970178275389170357551071104191609722674550547056511830750") call big_divrem (a, b, q, r) if (big_sgn (((b * q) + r) - a) /= 0) error stop a = string2big ("5286200801181288750950358142425694618335361315503743069535407838& &1104411448878793976933118436177295215313131557463887718741957154") b = string2big ("54401097470188014066128968444633185848791550678521") call big_divrem (a, b, q, r) if (big_sgn (((b * q) + r) - a) /= 0) error stop a = string2big ("74352827755975214935544861176217106447734695144315262422& &288346418457330023596489437154599028318030933202302606937951415862696330") b = string2big ("291979433784649910583546698460221489986784915256036707914578& &892106828527219639012712") call big_divrem (a, b, q, r) if (big_sgn (((b * q) + r) - a) /= 0) error stop a = string2big ("7515839498480018152556264500298705705667515770181724145893& &9544448656273749453586884931339958104923411455488844854494605760712247") b = string2big ("8600698996698965932302079501896131441135807557744554902970070& &402964318496325075877027770784963001") call big_divrem (a, b, q, r) if (big_sgn (((b * q) + r) - a) /= 0) error stop a = string2big ("13370595927769020368832742717678609885835798503146654175875& &149837801359758893206045930442389897206420586502996797614097489470778") b = string2big ("871343613388") call big_divrem (a, b, q, r) if (big_sgn (((b * q) + r) - a) /= 0) error stop end block golden_ratio = constant_term_cf (1) silver_ratio = constant_term_cf (2) one = i2cf (1) two = i2cf (2) three = i2cf (3) four = i2cf (4) sqrt2 = silver_ratio - one method1 = apply_ng8 (0, 1, 0, 0, 0, 0, 2, 0, silver_ratio, sqrt2) method2 = apply_ng8 (1, 0, 0, 1, 0, 0, 0, 8, silver_ratio, silver_ratio) method3 = (one / two) * (one + (one / sqrt2)) call show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2") call show ("silver ratio", silver_ratio, "(1 + sqrt(2))") call show ("sqrt(2)", sqrt2, "silver ratio minus 1") call show ("13/11", i2cf (13) / i2cf (11), "") call show ("22/7", i2cf (22) / i2cf (7), "") call show ("1", one, "") call show ("2", two, "") call show ("3", three, "") call show ("4", four, "") call show ("(1 + 1/sqrt(2))/2", method1, "method 1") call show ("(1 + 1/sqrt(2))/2", method2, "method 2") call show ("(1 + 1/sqrt(2))/2", method3, "method 3") call show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2, "") call show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2, "") call show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2, "") call show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2, "") contains subroutine show (expression, cf, note) character(len = *), intent(in) :: expression class(continued_fraction), intent(in) :: cf character(len = *), intent(in) :: note write (*, '(A19, " => ", A, T73, A)') & & expression, cf%to_string(), note end subroutine show end program bivariate_continued_fraction_task !---------------------------------------------------------------------  Output: $ gfortran -O2 -g -fbounds-check -Wall -Wextra bivariate_continued_fraction_task.f90 && ./a.out
golden ratio =>  [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]   (1 + sqrt(5))/2
silver ratio =>  [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]   (1 + sqrt(2))
sqrt(2) =>  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]   silver ratio minus 1
13/11 =>  [1;5,2]
22/7 =>  [3;7]
1 =>  [1]
2 =>  [2]
3 =>  [3]
4 =>  [4]
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 1
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 2
(1 + 1/sqrt(2))/2 =>  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]   method 3
sqrt(2) + sqrt(2) =>  [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) =>  [0]
sqrt(2) * sqrt(2) =>  [2]
sqrt(2) / sqrt(2) =>  [1]


Go

Adding to the existing package from the Continued_fraction/Arithmetic/Construct_from_rational_number#Go task, re-uses cf.go and rat.go as given in that task.

File ng8.go:

package cf

import "math"

// A 2×4 matix:
//     [ a₁₂   a₁   a₂   a ]
//     [ b₁₂   b₁   b₂   b ]
//
// which when "applied" to two continued fractions N1 and N2
// gives a new continued fraction z such that:
//
//         a₁₂ * N1 * N2  +  a₁ * N1  +  a₂ * N2  +  a
//     z = -------------------------------------------
//         b₁₂ * N1 * N2  +  b₁ * N1  +  b₂ * N2  +  b
//
// Examples:
//      NG8{0,1,1,0,  0,0,0,1} gives N1 + N2
//      NG8{0,1,-1,0, 0,0,0,1} gives N1 - N2
//      NG8{1,0,0,0,  0,0,0,1} gives N1 * N2
//      NG8{0,1,0,0,  0,0,1,0} gives N1 / N2
//      NG8{21,-15,28,-20, 0,0,0,1} gives 21*N1*N2 -15*N1 +28*N2 -20
//                               which is (3*N1 + 4) * (7*N2 - 5)
type NG8 struct {
A12, A1, A2, A int64
B12, B1, B2, B int64
}

// Basic identities as NG8 matrices.
var (
NG8Add = NG8{0, 1, 1, 0, 0, 0, 0, 1}
NG8Sub = NG8{0, 1, -1, 0, 0, 0, 0, 1}
NG8Mul = NG8{1, 0, 0, 0, 0, 0, 0, 1}
NG8Div = NG8{0, 1, 0, 0, 0, 0, 1, 0}
)

func (ng *NG8) needsIngest() bool {
if ng.B12 == 0 || ng.B1 == 0 || ng.B2 == 0 || ng.B == 0 {
return true
}
x := ng.A / ng.B
return ng.A1/ng.B1 != x || ng.A2/ng.B2 != x && ng.A12/ng.B12 != x
}

func (ng *NG8) isDone() bool {
return ng.B12 == 0 && ng.B1 == 0 && ng.B2 == 0 && ng.B == 0
}

func (ng *NG8) ingestWhich() bool { // true for N1, false for N2
if ng.B == 0 && ng.B2 == 0 {
return true
}
if ng.B == 0 || ng.B2 == 0 {
return false
}
x1 := float64(ng.A1) / float64(ng.B1)
x2 := float64(ng.A2) / float64(ng.B2)
x := float64(ng.A) / float64(ng.B)
return math.Abs(x1-x) > math.Abs(x2-x)
}

func (ng *NG8) ingest(isN1 bool, t int64) {
if isN1 {
// [ a₁₂   a₁   a₂   a ] becomes [ a₂+a₁₂*t  a+a₁*t  a₁₂  a₁]
// [ b₁₂   b₁   b₂   b ]         [ b₂+b₁₂*t  b+b₁*t  b₁₂  b₁]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.A2+ng.A12*t, ng.A+ng.A1*t, ng.A12, ng.A1,
ng.B2+ng.B12*t, ng.B+ng.B1*t, ng.B12, ng.B1
} else {
// [ a₁₂   a₁   a₂   a ] becomes [ a₁+a₁₂*t  a₁₂  a+a₂*t  a₂]
// [ b₁₂   b₁   b₂   b ]         [ b₁+b₁₂*t  b₁₂  b+b₂*t  b₂]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.A1+ng.A12*t, ng.A12, ng.A+ng.A2*t, ng.A2,
ng.B1+ng.B12*t, ng.B12, ng.B+ng.B2*t, ng.B2
}
}

func (ng *NG8) ingestInfinite(isN1 bool) {
if isN1 {
// [ a₁₂   a₁   a₂   a ] becomes [ a₁₂  a₁  a₁₂  a₁ ]
// [ b₁₂   b₁   b₂   b ]         [ b₁₂  b₁  b₁₂  b₁ ]
ng.A2, ng.A, ng.B2, ng.B =
ng.A12, ng.A1,
ng.B12, ng.B1
} else {
// [ a₁₂   a₁   a₂   a ] becomes [ a₁₂  a₁₂  a₂  a₂ ]
// [ b₁₂   b₁   b₂   b ]         [ b₁₂  b₁₂  b₂  b₂ ]
ng.A1, ng.A, ng.B1, ng.B =
ng.A12, ng.A2,
ng.B12, ng.B2
}
}

func (ng *NG8) egest(t int64) {
// [ a₁₂   a₁   a₂   a ] becomes [     b₁₂       b₁       b₂      b   ]
// [ b₁₂   b₁   b₂   b ]         [ a₁₂-b₁₂*t  a₁-b₁*t  a₂-b₂*t  a-b*t ]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.B12, ng.B1, ng.B2, ng.B,
ng.A12-ng.B12*t, ng.A1-ng.B1*t, ng.A2-ng.B2*t, ng.A-ng.B*t
}

// ApplyTo "applies" the matrix ng to the continued fractions
// N1 and N2, returning the resulting continued fraction.
// After ingesting limit terms without any output terms the resulting
// continued fraction is terminated.
func (ng NG8) ApplyTo(N1, N2 ContinuedFraction, limit int) ContinuedFraction {
return func() NextFn {
next1, next2 := N1(), N2()
done := false
sinceEgest := 0
return func() (int64, bool) {
if done {
return 0, false
}
for ng.needsIngest() {
sinceEgest++
if sinceEgest > limit {
done = true
return 0, false
}
isN1 := ng.ingestWhich()
next := next2
if isN1 {
next = next1
}
if t, ok := next(); ok {
ng.ingest(isN1, t)
} else {
ng.ingestInfinite(isN1)
}
}
sinceEgest = 0
t := ng.A / ng.B
ng.egest(t)
done = ng.isDone()
return t, true
}
}
}


File ng8_test.go:

package cf

import "fmt"

func ExampleNG8() {
cases := [...]struct {
op     string
r1, r2 Rat
ng     NG8
}{
{"+", Rat{22, 7}, Rat{1, 2}, NG8Add},
{"*", Rat{13, 11}, Rat{22, 7}, NG8Mul},
{"-", Rat{13, 11}, Rat{22, 7}, NG8Sub},
{"/", Rat{22 * 22, 7 * 7}, Rat{22, 7}, NG8Div},
}
for _, tc := range cases {
n1 := tc.r1.AsContinuedFraction()
n2 := tc.r2.AsContinuedFraction()
z := tc.ng.ApplyTo(n1, n2, 1000)
fmt.Printf("%v %s %v is %v %v %v gives %v\n",
tc.r1, tc.op, tc.r2,
tc.ng, n1, n2, z,
)
}

z := NG8Mul.ApplyTo(Sqrt2, Sqrt2, 1000)
fmt.Println("√2 * √2 =", z)

// Output:
// 22/7 + 1/2 is {0 1 1 0 0 0 0 1} [3; 7] [0; 2] gives [3; 1, 1, 1, 4]
// 13/11 * 22/7 is {1 0 0 0 0 0 0 1} [1; 5, 2] [3; 7] gives [3; 1, 2, 2]
// 13/11 - 22/7 is {0 1 -1 0 0 0 0 1} [1; 5, 2] [3; 7] gives [-1; -1, -24, -1, -2]
// 484/49 / 22/7 is {0 1 0 0 0 0 1 0} [9; 1, 7, 6] [3; 7] gives [3; 7]
// √2 * √2 = [1; 0, 1]
}

Output:

(Note, [1; 0, 1] = 1 + 1 / (0 + 1/1 ) = 2, so the answer is correct, it should however be normalised to the more reasonable form of [2].)

22/7 + 1/2 is {0 1 1 0 0 0 0 1} [3; 7] [0; 2] gives [3; 1, 1, 1, 4]
13/11 * 22/7 is {1 0 0 0 0 0 0 1} [1; 5, 2] [3; 7] gives [3; 1, 2, 2]
13/11 - 22/7 is {0 1 -1 0 0 0 0 1} [1; 5, 2] [3; 7] gives [-1; -1, -24, -1, -2]
484/49 / 22/7 is {0 1 0 0 0 0 1 0} [9; 1, 7, 6] [3; 7] gives [3; 7]
√2 * √2 = [1; 0, 1]


Translation of: Mercury

This Haskell follows the Mercury, in using infinitely long lazy lists to represent continued fractions. There are two kinds of terms: "infinite" and "finite integer".

----------------------------------------------------------------------

data Term = InfiniteTerm | IntegerTerm Integer
type ContinuedFraction = [Term] -- The list should be infinitely long.

type NG8 = (Integer, Integer, Integer, Integer,
Integer, Integer, Integer, Integer)

----------------------------------------------------------------------

cf2string (cf :: ContinuedFraction) =
loop 0 "[" cf
where loop i s lst =
case lst of {
(InfiniteTerm : _) -> s ++ "]" ;
(IntegerTerm value : tail) ->
(if i == 20 then
s ++ ",...]"
else
let {
sepStr =
case i of {
0 -> "";
1 -> ";";
_ -> ","
};
termStr = show value;
s1 = s ++ sepStr ++ termStr
}
in loop (i + 1) s1 tail)
}

----------------------------------------------------------------------

repeatingTerm (term :: Term) =
term : repeatingTerm term

infiniteContinuedFraction = repeatingTerm InfiniteTerm

i2cf (i :: Integer) =
-- Continued fraction representing an integer.
IntegerTerm i : infiniteContinuedFraction

r2cf (n :: Integer) (d :: Integer) =
-- Continued fraction representing a rational number.
let (q, r) = divMod n d in
(if r == 0 then
(IntegerTerm q : infiniteContinuedFraction)
else
(IntegerTerm q : r2cf d r))

----------------------------------------------------------------------

add_cf = apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1)
sub_cf = apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1)
mul_cf = apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1)
div_cf = apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0)

apply_ng8
(ng :: NG8)
(x :: ContinuedFraction)
(y :: ContinuedFraction) =
--
let (a12, a1, a2, a, b12, b1, b2, b) = ng in
if iseqz [b12, b1, b2, b] then
infiniteContinuedFraction -- No more finite terms to output.
else if iseqz [b2, b] then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else if atLeastOne_iseqz [b2, b] then
let (ng1, x1, y1) = absorb_y_term ng x y in
apply_ng8 ng1 x1 y1
else if iseqz [b1] then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else
let {
(q12, r12) = maybeDivide a12 b12;
(q1, r1) = maybeDivide a1 b1;
(q2, r2) = maybeDivide a2 b2;
(q, r) = maybeDivide a b
}
in
if not (iseqz [b12]) && q == q12 && q == q1 && q == q2 then
-- Output a term.
(if integerExceedsInfinitizingThreshold q then
infiniteContinuedFraction
else
let new_ng = (b12, b1, b2, b, r12, r1, r2, r) in
(IntegerTerm q : apply_ng8 new_ng x y))
else
-- Put a1, a2, and a over a common denominator and compare
-- some magnitudes.
let {
n1 = a1 * b2 * b;
n2 = a2 * b1 * b;
n = a * b1 * b2
}
in
(if abs (n1 - n) > abs (n2 - n) then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else
let (ng1, x1, y1) = absorb_y_term ng x y in
apply_ng8 ng1 x1 y1)

absorb_x_term
((a12, a1, a2, a, b12, b1, b2, b) :: NG8)
(x :: ContinuedFraction)
(y :: ContinuedFraction) =
--
case x of {
(IntegerTerm n : xtail) -> (
let new_ng = (a2 + (a12 *`