Continued fraction/Arithmetic/G(matrix ng, continued fraction n)

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

This task investigates mathmatical operations that can be performed on a single continued fraction. This requires only a baby version of NG:

I may perform perform the following operations:

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

I output a term if the integer parts of and are equal. Otherwise I input a term from N. If I need a term from N but N has no more terms I inject .

When I input a term t my internal state: is transposed thus

When I output a term t my internal state: is transposed thus

When I need a term t but there are no more my internal state: is transposed thus

I am done when b1 and b are zero.

Demonstrate your solution by calculating:

[1;5,2] + 1/2
[3;7] + 1/2
[3;7] divided by 4

Using a generator for (e.g., from Continued fraction) calculate . You are now at the starting line for using Continued Fractions to implement Arithmetic-geometric mean without ulps and epsilons.

The first step in implementing Arithmetic-geometric mean is to calculate do this now to cross the starting line and begin the race.

11l

Translation of: Python
T NG
   Int a1, a, b1, b
   F (a1, a, b1, b)
      .a1 = a1
      .a  = a
      .b1 = b1
      .b  = b

   F ingress(n)
      (.a, .a1) = (.a1, .a + .a1 * n)
      (.b, .b1) = (.b1, .b + .b1 * n)

   F needterm()
      R (.b == 0 | .b1 == 0) | !(.a I/ .b == .a1 I/ .b1)

   F egress()
      V n = .a I/ .b
      (.a, .b) = (.b, .a - .b * n)
      (.a1, .b1) = (.b1, .a1 - .b1 * n)
      R n

   F egress_done()
      I .needterm()
         (.a, .b) = (.a1, .b1)
      R .egress()

   F done()
      R .b == 0 & .b1 == 0

F r2cf(=n1, =n2)
   [Int] r
   L n2 != 0
      (n1, V t1_n2) = (n2, divmod(n1, n2))
      n2 = t1_n2[1]
      r [+]= t1_n2[0]
   R r

V data = [(‘[1;5,2] + 1/2’,      2,1,0,2, (13, 11)),
          (‘[3;7] + 1/2’,        2,1,0,2, (22,  7)),
          (‘[3;7] divided by 4’, 1,0,0,4, (22,  7))]

L(string, a1, a, b1, b, r) data
   print(‘#<20->’.format(string), end' ‘’)
   V op = NG(a1, a, b1, b)
   L(n) r2cf(r[0], r[1])
      I !op.needterm()
         print(‘ ’op.egress(), end' ‘’)
      op.ingress(n)
   L
      print(‘ ’op.egress_done(), end' ‘’)
      I op.done()
         L.break
   print()
Output:
[1;5,2] + 1/2       -> 1 1 2 7
[3;7] + 1/2         -> 3 1 1 1 4
[3;7] divided by 4  -> 0 1 3 1 2

Ada

Translation of: ATS
Translation of: C
----------------------------------------------------------------------

with ada.text_io; use ada.text_io;
with ada.strings; use ada.strings;
with ada.strings.fixed; use ada.strings.fixed;

with ada.unchecked_deallocation;

procedure univariate_continued_fraction_task is

--  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
-- A generator is a tree structure, accessed by a generator_t.

type generator_record;
type generator_t is access generator_record;

type generator_list_node;
type generator_list_t is access generator_list_node;

type generator_list_node is
  record
    car : generator_t;
    cdr : generator_list_t;
  end record;

type generator_workspace is array (natural range <>) of integer;
type generator_worksp_t is access generator_workspace;

type generator_proc_t is
  access procedure (workspace   : in generator_worksp_t;
                    sources     : in generator_list_t;
                    term_exists : out boolean;
                    term        : out integer);

type generator_record is
  record
    run       : generator_proc_t;   -- What does the work.
    worksize  : natural;            -- The size of workspace.
    initial   : generator_worksp_t; -- The initial value of workspace.
    workspace : generator_worksp_t; -- The workspace.
    sources   : generator_list_t;   -- The sources of input terms.
  end record;

procedure free_generator_workspace is
  new ada.unchecked_deallocation (generator_workspace,
                                  generator_worksp_t);

procedure free_generator_list_node is
  new ada.unchecked_deallocation (generator_list_node,
                                  generator_list_t);

procedure free_generator_record is
  new ada.unchecked_deallocation (generator_record,
                                  generator_t);

procedure initialize_workspace (gen : in generator_t) is
begin
  for i in 0 .. gen.worksize - 1 loop
    gen.workspace(i) := gen.initial(i);
  end loop;
end initialize_workspace;

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

-- Re-initialize a generator for a computation.
procedure initialize_generator_t (gen : in generator_t) is
  p : generator_list_t;
begin
  if gen /= null then
    initialize_workspace (gen);

    p := gen.sources;
    while p /= null loop
      initialize_generator_t (p.car);
      p := p.cdr;
    end loop;
  end if;
end initialize_generator_t;

-- Free the storage of a generator.
procedure free_generator_t (gen : in out generator_t) is
  p, q : generator_list_t;
begin
  if gen /= null then
    free_generator_workspace (gen.initial);
    free_generator_workspace (gen.workspace);

    p := gen.sources;
    while p /= null loop
      q := p.cdr;
      free_generator_t (p.car);
      free_generator_list_node (p);
      p := q;
    end loop;

    free_generator_record (gen);
  end if;
end free_generator_t;

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

-- Run a generator and print its output.
procedure put_generator_output (gen       : in generator_t;
                                max_terms : in positive) is
  sep         : integer range 0 .. 2;
  terms_count : natural;
  term_exists : boolean;
  term        : integer;
  done        : boolean;
begin
  initialize_generator_t (gen);

  terms_count := 0;
  sep := 0;
  done := false;
  while not done loop
    if terms_count = max_terms then
      put (",...]");
      done := true;
    else
      gen.run (gen.workspace, gen.sources, term_exists, term);
      if term_exists then
        case sep is
          when 0 =>
            put ("[");
            sep := 1;
          when 1 =>
            put (";");
            sep := 2;
          when 2 =>
            put (",");
        end case;
        put (trim (integer'image (term), left));
        terms_count := terms_count + 1;
      else
        put ("]");
        done := true;
      end if;
    end if;
  end loop;

  initialize_generator_t (gen);
end put_generator_output;

--  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
-- Generators for continued fraction terms of rational numbers.

procedure r2cf_run (workspace   : in generator_worksp_t;
                    sources     : in generator_list_t;
                    term_exists : out boolean;
                    term        : out integer) is
  n, d, q, r : integer;
begin
  d := workspace(1);
  term_exists := (d /= 0);
  if term_exists then
    n := workspace(0);

    -- We shall use the kind of integer division that is most
    -- "natural" in Ada: truncation towards zero.
    q := n / d;                 -- Truncation towards zero.
    r := n rem d;               -- The remainder may be negative.

    workspace(0) := d;
    workspace(1) := r;

    term := q;
  end if;
end r2cf_run;

-- Make a generator for the fraction n/d.
function r2cf_make (n : in integer;
                    d : in integer)
return generator_t is
  gen : generator_t := new generator_record;
begin
  gen.run := r2cf_run'access;
  gen.worksize := 2;
  gen.initial := new generator_workspace(0 .. gen.worksize - 1);
  gen.workspace := new generator_workspace(0 .. gen.worksize - 1);
  gen.initial(0) := n;
  gen.initial(1) := d;
  initialize_workspace (gen);
  gen.sources := null;
  return gen;
end r2cf_make;

--  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
-- A generator for continued fraction terms of sqrt(2).

procedure sqrt2_run (workspace   : in generator_worksp_t;
                     sources     : in generator_list_t;
                     term_exists : out boolean;
                     term        : out integer) is
begin
  term_exists := true;
  term := workspace(0);
  workspace(0) := 2;
end sqrt2_run;

-- Make a generator for the fraction n/d.
function sqrt2_make
return generator_t is
  gen : generator_t := new generator_record;
begin
  gen.run := sqrt2_run'access;
  gen.worksize := 1;
  gen.initial := new generator_workspace(0 .. gen.worksize - 1);
  gen.workspace := new generator_workspace(0 .. gen.worksize - 1);
  gen.initial(0) := 1;
  initialize_workspace (gen);
  gen.sources := null;
  return gen;
end sqrt2_make;

--  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
-- Generators for the application of homographic functions to other
-- generators.

procedure hfunc_take_term (workspace : in generator_worksp_t;
                           sources   : in generator_list_t) is
  term_exists1 : boolean;
  term1        : integer;
  src          : generator_t := sources.car;
  a1, a, b1, b : integer;
begin
  src.run (src.workspace, src.sources, term_exists1, term1);
  a1 := workspace(0);
  b1 := workspace(2);
  if term_exists1 then
    a := workspace(1);
    b := workspace(3);
    workspace(0) := a + (a1 * term1);
    workspace(1) := a1;
    workspace(2) := b + (b1 * term1);
    workspace(3) := b1;
  else
    workspace(1) := a1;
    workspace(3) := b1;
  end if;
end hfunc_take_term;

procedure hfunc_run (workspace   : in generator_worksp_t;
                     sources     : in generator_list_t;
                     term_exists : out boolean;
                     term        : out integer) is
  done         : boolean;
  a1, a, b1, b : integer;
  q1, q        : integer;
begin
  done := false;
  while not done loop
    b1 := workspace(2);
    b := workspace(3);
    if b1 = 0 and b = 0 then
      term_exists := false;
      done := true;
    else
      a1 := workspace(0);
      a := workspace(1);
      if b1 /= 0 and b /= 0 then
        q1 := a1 / b1;
        q := a / b;
        if q1 = q then
          workspace(0) := b1;
          workspace(1) := b;
          workspace(2) := a1 - (b1 * q);
          workspace(3) := a - (b * q);
          term_exists := true;
          term := q;
          done := true;
        else
          hfunc_take_term (workspace, sources);
        end if;
      else
        hfunc_take_term (workspace, sources);
      end if;
    end if;
  end loop;
end hfunc_run;

function hfunc_make (a1     : in integer;
                     a      : in integer;
                     b1     : in integer;
                     b      : in integer;
                     source : in generator_t)
return generator_t is
  gen : generator_t := new generator_record;
begin
  gen.run := hfunc_run'access;
  gen.worksize := 4;
  gen.initial := new generator_workspace(0 .. gen.worksize - 1);
  gen.workspace := new generator_workspace(0 .. gen.worksize - 1);
  gen.initial(0) := a1;
  gen.initial(1) := a;
  gen.initial(2) := b1;
  gen.initial(3) := b;
  initialize_workspace (gen);
  gen.sources := new generator_list_node;
  gen.sources.car := source;
  gen.sources.cdr := null;
  return gen;
end hfunc_make;

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

max_terms : constant positive := 20;

procedure run_gen (expr : in string;
                   gen  : in generator_t) is
begin
  put (expr);
  put (" => ");
  put_generator_output (gen, max_terms);
  put_line ("");
end run_gen;

gen : generator_t;

begin

  gen := r2cf_make (13, 11);
  run_gen ("13/11", gen);
  free_generator_t (gen);

  gen := r2cf_make (22, 7);
  run_gen ("22/7", gen);
  free_generator_t (gen);

  gen := sqrt2_make;
  run_gen ("sqrt(2)", gen);
  free_generator_t (gen);

  gen := hfunc_make (2, 1, 0, 2, r2cf_make (13, 11));
  run_gen ("13/11 + 1/2", gen);
  free_generator_t (gen);

  gen := hfunc_make (2, 1, 0, 2, r2cf_make (22, 7));
  run_gen ("22/7 + 1/2", gen);
  free_generator_t (gen);

  gen := hfunc_make (1, 0, 0, 4, r2cf_make (22, 7));
  run_gen ("(22/7)/4", gen);
  free_generator_t (gen);

  gen := hfunc_make (1, 0, 0, 2, sqrt2_make);
  run_gen ("sqrt(2)/2", gen);
  free_generator_t (gen);

  gen := hfunc_make (0, 1, 1, 0, sqrt2_make);
  run_gen ("1/sqrt(2)", gen);
  free_generator_t (gen);

  gen := hfunc_make (1, 2, 0, 4, sqrt2_make);
  run_gen ("(2 + sqrt(2))/4", gen);
  free_generator_t (gen);

  -- Demonstrate that you can compose homographic functions:
  gen := hfunc_make (1, 0, 0, 2,
                     hfunc_make (1, 1, 0, 1,
                                 hfunc_make (0, 1, 1, 0,
                                             sqrt2_make)));
  run_gen ("(1 + 1/sqrt(2))/2", gen);
  free_generator_t (gen);

end univariate_continued_fraction_task;

----------------------------------------------------------------------
Output:
$ gnatmake -q univariate_continued_fraction_task.adb && ./univariate_continued_fraction_task
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

ATS

Using non-linear types

The approach used here leaks memory freely, and so a program using it may require a garbage collector. An advantage, however, is that it memoizes results.

For no reason except my curiosity, this demo program outputs its results in MathML.

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

#include "share/atspre_staload.hats"

(* We need consistent definitions of division and remainder. Let us
   set those here. For convenience (because the prelude provides it),
   we will use truncation towards zero. *)
infixl ( / ) div
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod

(* Some useful math Unicode, so we do not need entities in markup
   using these characters. *)
val plus_sign = "&#x002B;"
val dot_operator = "&#x22C5;"
val centered_ellipsis = "&#x22EF;"
val right_arrow = "&#x2192;"

(*------------------------------------------------------------------*)
(* Continued fractions as processes for generating terms. The terms
   are memoized and are accessed by their zero-based index. The terms
   are represented as any one of the signed integer types for which
   there is a typekind. *)

abstype cf (tk : tkind) = ptr

(* A ref of a cf has the advantage, over a cf itself, that it can be
   more safely used in a closure. *)
typedef cfref (tk : tkind) = ref (cf tk)

typedef cf_generator (tk : tkind) =
  () -<cloref1> Option (g0int tk)

extern fn {tk : tkind}
cf_make :
  cf_generator tk -> cf tk

extern fn {tk : tkind}
cfref_make :
  cf_generator tk -> cfref tk

extern fn {tk  : tkind}
          {tki : tkind}
cf_get_at_guint :
  {i : int}
  (&cf tk >> _, g1uint (tki, i)) -> Option (g0int tk)

extern fn {tk  : tkind}
          {tki : tkind}
cf_get_at_gint :
  {i : nat}
  (&cf tk >> _, g1int (tki, i)) -> Option (g0int tk)

overload cf_get_at with cf_get_at_guint
overload cf_get_at with cf_get_at_gint
overload [] with cf_get_at

macdef cfref_get_at (cfref, i) =
  let
    val @(pf, fpf | p) = ref_vtakeout ,(cfref)
    val retval = cf_get_at (!p, ,(i))
    prval () = fpf pf
  in
    retval
  end

extern fn {tk : tkind}
cf2mathml_max_terms
          (cf        : &cf tk >> _,
           max_terms : size_t)
    : string

extern fn {tk : tkind}
cf2mathml_default_max_terms
          (cf : &cf tk >> _)
    : string

extern fn {tk : tkind}
cfref2mathml_max_terms
          (cfref     : cfref tk,
           max_terms : size_t)
    : string

extern fn {tk : tkind}
cfref2mathml_default_max_terms
          (cfref : cfref tk)
    : string

overload cf2mathml with cf2mathml_max_terms
overload cf2mathml with cf2mathml_default_max_terms
overload cfref2mathml with cfref2mathml_max_terms
overload cfref2mathml with cfref2mathml_default_max_terms

(* To use a cfref as a generator, you need cfref2generator. It would
   do no good to use the cf object's internal generator directly,
   because its state would be wrong. In any case, the internals of a
   cf are hidden from the programmer. *)
extern fn {tk : tkind}
cfref2generator :
  cfref tk -> cf_generator tk

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

local

  typedef _cf (tk         : tkind,
               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 (g0int tk, n), (* Memoized terms. *)
      gen = cf_generator tk          (* A thunk to generate terms. *)
    }
  typedef _cf (tk : tkind, m : int) =
    [terminated : bool]
    [n : int | m <= n]
    _cf (tk, terminated, m, n)
  typedef _cf (tk : tkind) =
    [m : int]
    _cf (tk, m)

  fn {tk : tkind}
  _cf_get_more_terms
            {terminated : bool}
            {m          : int}
            {n          : int}
            {needed     : int | m <= needed; needed <= n}
            (cf         : _cf (tk, terminated, m, n),
             needed     : size_t needed)
      : [m1 : int | m <= m1; m1 <= needed]
        [n1 : int | m1 <= n1]
        _cf (tk, 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 (tk, m1 < needed, m1, n1) =
        if i = needed then
          '{
            terminated = false,
            m = needed,
            n = (cf.n),
            memo = memo,
            gen = (cf.gen)
          }
        else
          begin
            case+ (cf.gen) () of
            | None () =>
              '{
                terminated = true,
                m = i,
                n = (cf.n),
                memo = memo,
                gen = (cf.gen)
              }
            | Some term =>
              begin
                memo[i] := term;
                loop (succ i)
              end
          end
    in
      loop (cf.m)
    end

  fn {tk : tkind}
  _cf_update {terminated : bool}
             {m          : int}
             {n          : int | m <= n}
             {needed     : int}
             (cf         : _cf (tk, terminated, m, n),
              needed     : size_t needed)
      : [terminated1 : bool]
        [m1 : int | m <= m1]
        [n1 : int | m1 <= n1]
        _cf (tk, 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<tk> (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<tk> (cf1, needed)
        end
    end

in (* local *)

  assume cf tk = _cf tk

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

  implement {tk}
  cfref_make gen =
    ref (cf_make gen)

  implement {tk} {tki}
  cf_get_at_guint {i} (cf, i) =
    let
      prval () = lemma_g1uint_param i
      val i : size_t i = g1u2u i
      val cf1 = _cf_update<tk> (cf, succ i)
    in
      cf := cf1;
      if i < (cf1.m) then
        Some (arrayref_get_at<g0int tk> (cf1.memo, i))
      else
        None ()
    end

  implement {tk} {tki}
  cf_get_at_gint (cf, i) =
    cf_get_at_guint<tk><sizeknd> (cf, g1i2u i)

end (* local *)

implement {tk}
cf2mathml_max_terms (cf, max_terms) =
  let
    fun
    loop (i     : Size_t,
          cf    : &cf tk >> _,
          sep   : string,
          accum : string)
        : string =
      if i = max_terms then
        strptr2string
          (string_append
            (accum, "<mo>,</mo><mo>",
             centered_ellipsis, "</mo><mo>]</mo>"))
      else
        begin
          case+ cf[i] of
          | None () =>
            strptr2string (string_append (accum, "<mo>]</mo>"))
          | Some term =>
            let
              val term_str = tostring_val<g0int tk> term
              val accum =
                strptr2string (string_append (accum, sep, "<mn>",
                                              term_str, "</mn>"))
              val sep =
                if sep = "<mo>[</mo>" then
                  "<mo>;</mo>"
                else
                  "<mo>,</mo>"
            in
              loop (succ i, cf, sep, accum)
            end
        end
  in
    loop (i2sz 0, cf, "<mo>[</mo>", "")
  end

implement {tk}
cf2mathml_default_max_terms cf =
  let
    #ifndef DEFAULT_CF_MAX_TERMS #then
      #define DEFAULT_CF_MAX_TERMS 20
    #endif
  in
    cf2mathml_max_terms (cf, i2sz DEFAULT_CF_MAX_TERMS)
  end

implement {tk}
cfref2mathml_max_terms (cfref, max_terms) =
  let
    val @(pf, fpf | p) = ref_vtakeout cfref
    val retval = cf2mathml_max_terms (!p, max_terms)
    prval () = fpf pf
  in
    retval
  end

implement {tk}
cfref2mathml_default_max_terms cfref =
  let
    val @(pf, fpf | p) = ref_vtakeout cfref
    val retval = cf2mathml_default_max_terms !p
    prval () = fpf pf
  in
    retval
  end

implement {tk}
cfref2generator cfref =
  let
    val index : ref Size_t = ref (i2sz 0)
  in
    lam () =>
      let
        val i = !index
        val retval = cfref_get_at (cfref, i)
      in
        !index := succ i;
        retval
      end
  end

(*------------------------------------------------------------------*)
(* A homographic function. *)

typedef hfunc (tk : tkind) =
  @{
    a1 = g0int tk,
    a = g0int tk,
    b1 = g0int tk,
    b = g0int tk
  }

extern fn {tk : tkind}
hfunc_make :
  (g0int tk, g0int tk, g0int tk, g0int tk) -<> hfunc tk

extern fn {tk : tkind}
hfunc_apply_generator2generator :
  (hfunc tk, cf_generator tk) -> cf_generator tk

extern fn {tk : tkind}
hfunc_apply_cfref2cfref :
  (hfunc tk, cfref tk) -> cfref tk

overload hfunc_apply with hfunc_apply_generator2generator
overload hfunc_apply with hfunc_apply_cfref2cfref

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

implement {tk}
hfunc_make (a1, a, b1, b) =
  @{
    a1 = a1,
    a = a,
    b1 = b1,
    b = b
  }

fn {tk : tkind}
take_term_from_ngen
          (state : ref (hfunc tk),
           ngen  : cf_generator tk)
    : void =
  let
    val @{
          a1 = a1,
          a = a,
          b1 = b1,
          b = b
        } = !state
  in
    case+ ngen () of
    | Some term =>
      !state :=
        @{
          a1 = a + (a1 * term),
          a = a1,
          b1 = b + (b1 * term),
          b = b1
        }
    | None () =>
      !state :=
        @{
          a1 = a1,
          a = a1,
          b1 = b1,
          b = b1
        }
  end

fn {tk : tkind}
adjust_state_for_term_output
          (state : ref (hfunc tk),
           term  : g0int tk)
    : void =
  let
    val @{
          a1 = a1,
          a = a,
          b1 = b1,
          b = b
        } = !state
  in
    !state :=
      @{
        a1 = b1,
        a = b,
        b1 = a1 - (b1 * term),
        b = a - (b * term)
      }
  end

implement {tk}
hfunc_apply_generator2generator (f, ngen) =
  let
    val state : ref (hfunc tk) = ref f

    val hgen =
      lam () =<cloref1>
        let
          fun
          loop () : Option (g0int tk) =
            let
              val b1_iseqz = iseqz (!state.b1)
              and b_iseqz = iseqz (!state.b)
            in
              if b1_iseqz * b_iseqz then
                None ()
              else if (~b1_iseqz) * (~b_iseqz) then
                let
                  val q1 = (!state.a1) div (!state.b1)
                  and q = (!state.a) div (!state.b)
                in
                  if q1 = q then
                    begin
                      adjust_state_for_term_output<tk> (state, q);
                      Some q
                    end
                  else
                    begin
                      take_term_from_ngen (state, ngen);
                      loop ()
                    end
                end
              else
                begin
                  take_term_from_ngen (state, ngen);
                  loop ()
                end
            end
        in
          loop ()
        end
  in
    hgen : cf_generator tk
  end

implement {tk}
hfunc_apply_cfref2cfref (f, cfref) =
  let
    val gen1 = cfref2generator<tk> cfref
    val gen2 = hfunc_apply_generator2generator<tk> (f, gen1)
  in
    cfref_make<tk> gen2
  end

(*------------------------------------------------------------------*)
(* Let us create some continued fractions. *)

extern fn {tk : tkind}
r2cf :
  (g0int tk, g0int tk) -> cf tk

implement {tk}
r2cf (n, d) =
  let
    val n = ref<g0int tk> n
    val d = ref<g0int tk> d

    fn
    gen () :<cloref1> Option (g0int tk) =
      if iseqz !d then
        None ()
      else
        let
          val @(numer, denom) = @(!n, !d)
          val q = numer div denom
          and r = numer rem denom
        in
          !n := denom;
          !d := r;
          Some q
        end
  in
    cf_make gen
  end

val cfref_13_11 = ref (r2cf (13LL, 11LL)) (* 13/11 = [1;5,2] *)
val cfref_22_7 = ref (r2cf (22LL, 7LL))   (* 22/7 = [3;7] *)
val cfref_sqrt2 =               (* sqrt(2) = [1;2,2,2,...] *)
  let
    val term : ref llint = ref 1LL
    val gen =
      lam () =<cloref1>
        let
          val retval = !term
        in
          if retval = 1LL then
            !term := 2LL;
          Some retval
        end
  in
    cfref_make (gen : cf_generator llintknd)
  end

(*------------------------------------------------------------------*)
(* Let us create some homographic functions that correspond to unary
   arithmetic operations. *)

val add_one_half = hfunc_make (2LL, 1LL, 0LL, 2LL)
val add_one = hfunc_make (1LL, 1LL, 0LL, 1LL)
val divide_by_two = hfunc_make (1LL, 0LL, 0LL, 2LL)
val divide_by_four = hfunc_make (1LL, 0LL, 0LL, 4LL)
val take_reciprocal = hfunc_make (0LL, 1LL, 1LL, 0LL)
val add_one_then_div_two = hfunc_make (1LL, 1LL, 0LL, 2LL)
val add_two_then_div_four = hfunc_make (1LL, 2LL, 0LL, 4LL)

(*------------------------------------------------------------------*)
(* Now let us derive some continued fractions. *)

local
  macdef apply = hfunc_apply<llintknd>
  macdef cfref2ml = cfref2mathml<llintknd>
in
  val cfref_13_11_plus_1_2 = apply (add_one_half, cfref_13_11)
  val cfref_22_7_plus_1_2 = apply (add_one_half, cfref_22_7)
  val cfref_22_7_div_4 = apply (divide_by_four, cfref_22_7)

  (* The following two give the same result: *)
  val cfref_sqrt2_div_2 = apply (divide_by_two, cfref_sqrt2)
  val cfref_1_div_sqrt2 = apply (take_reciprocal, cfref_sqrt2)
  val () = assertloc (cfref2ml cfref_sqrt2_div_2
                          = cfref2ml cfref_1_div_sqrt2)

  (* The following three give the same result: *)
  val cfref_2_plus_sqrt2_grouped_div_4 =
    apply (add_two_then_div_four, cfref_sqrt2)
  val cfref_half_of_1_plus_half_sqrt2 =
    apply (add_one_then_div_two, cfref_sqrt2_div_2)
  val cfref_half_of_1_plus_1_div_sqrt2 =
    apply (divide_by_two, (apply (add_one, cfref_sqrt2_div_2)))
  val () = assertloc (cfref2ml cfref_2_plus_sqrt2_grouped_div_4
                          = cfref2ml cfref_half_of_1_plus_half_sqrt2)
  val () = assertloc (cfref2ml cfref_half_of_1_plus_half_sqrt2
                          = cfref2ml cfref_half_of_1_plus_1_div_sqrt2)
end

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

implement
main () =
  let
    macdef cfref2ml = cfref2mathml<llintknd>
    macdef apply = hfunc_apply<llintknd>

    macdef text (s) =
      strptr2string (string_append ("<p>", ,(s), "</p>"))

    macdef becomes =
      strptr2string (string_append ("<mo>", right_arrow, "</mo>"))

    macdef start_math =
      "<math xmlns='http://www.w3.org/1998/Math/MathML'>"
    macdef stop_math = "</math>"

    macdef start_table = "<mtable>"
    macdef stop_table = "</mtable>"

    macdef left_side (s) = 
      strptr2string
        (string_append
          ("<mtd columnalign='right'>", ,(s), "</mtd>"))
    macdef middle (s) = 
      strptr2string
        (string_append
          ("<mtd columnalign='center'>", ,(s), "</mtd>"))
    macdef right_side (s) = 
      strptr2string
        (string_append
          ("<mtd columnalign='left'>", ,(s), "</mtd>"))
    macdef entry (left, right) =
      strptr2string
        (string_append
          ("<mtr>",
           left_side (,(left)),
           middle becomes,
           right_side (,(right)),
           "</mtr>"))

    macdef num s =
      strptr2string (string_append ("<mn>", ,(s), "</mn>"))
    macdef id s =
      strptr2string (string_append ("<mi>", ,(s), "</mi>"))
    macdef oper s =
      strptr2string (string_append ("<mo>", ,(s), "</mo>"))

    macdef frac (n, d) = 
      strptr2string (string_append ("<mfrac>", ,(n), ,(d),
                                    "</mfrac>"))
    macdef numfrac (n, d) = frac (num ,(n), num ,(d))

    macdef sqrt_of (s) =
      strptr2string (string_append ("<msqrt>", ,(s), "</msqrt>"))
  in
    println! (start_math);
    println! (start_table);

    println! (entry (numfrac ("13", "11"), cfref2ml cfref_13_11));
    println! (entry (numfrac ("22", "7"), cfref2ml cfref_22_7));
    println! (entry (sqrt_of (num "2"), cfref2ml cfref_sqrt2));

    println! (entry (strptr2string
                        (string_append (numfrac ("13", "11"),
                                        oper plus_sign,
                                        numfrac ("1", "2"))),
                     cfref2ml cfref_13_11_plus_1_2));
    println! (entry (strptr2string
                        (string_append (numfrac ("22", "7"),
                                        oper plus_sign,
                                        numfrac ("1", "2"))),
                     cfref2ml cfref_22_7_plus_1_2));
    println! (entry (frac (numfrac ("22", "7"), num ("4")),
                     cfref2ml cfref_22_7_div_4));
    println! (entry (frac (sqrt_of (num "2"), num ("2")),
                     cfref2ml cfref_sqrt2_div_2));
    println! (entry (frac (num ("1"), sqrt_of (num "2")),
                     cfref2ml cfref_1_div_sqrt2));
    println! (entry (strptr2string
                        (string_append
                            (numfrac ("1", "4"),
                             oper dot_operator,
                             strptr2string
                                (string_append
                                  (oper "(", num "2", oper plus_sign,
                                   sqrt_of (num "2"), oper ")")))),
                     cfref2ml cfref_2_plus_sqrt2_grouped_div_4));
    println! (entry (strptr2string
                        (string_append
                            (numfrac ("1", "2"),
                             oper dot_operator,
                             strptr2string
                                (string_append
                                  (oper "(", num "1", oper plus_sign,
                                   frac (sqrt_of (num "2"), num "2"),
                                   oper ")")))),
                     cfref2ml cfref_half_of_1_plus_half_sqrt2));
    println! (entry (strptr2string
                        (string_append
                            (numfrac ("1", "2"),
                             oper dot_operator,
                             strptr2string
                                (string_append
                                  (oper "(", num "1", oper plus_sign,
                                   frac (num "1", sqrt_of (num "2")),
                                   oper ")")))),
                     cfref2ml cfref_half_of_1_plus_1_div_sqrt2));

    println! (stop_table);
    println! (stop_math);
    0
  end

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

Run something such as:

$ patscc -DATS_MEMALLOC_GCBDW -O2 -std=gnu2x univariate-continued-fraction-task.dats -lgc
$ ./a.out > foo.html
$ firefox foo.html

You can use a different browser, but it might not render the MathML in the way Firefox does here:

Output:

The output from the program, as rendered by Firefox.

Using linear types

This method of implementation purposely avoids the need for a garbage collector, and should be safe against the possibility of a memory leak. It does not memoize results, however. Memoization could be added, but effective safe use of it might require a host of other features, such as "generator splitters" or reference counting. By safe I mean safe against memory leaks and double-freeing. These are issues that do not arise in a program with garbage collection.

The demo program outputs LuaTeX macro code (for plain TeX, not LaTeX).

One thing a person may notice is the opt_some/opt_none/opt_unsome/opt_unnone: this is compile-time safety against the possibility of using an uninitialized return value, when the return is by reference parameter.

(*------------------------------------------------------------------*)
(* In this implementation, memory is allocated only for linear
   types. Thus discipline is maintained in the freeing of allocated
   space. There is, however, no memoization. *)
(*------------------------------------------------------------------*)

#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"

(* We need consistent definitions of division and remainder. Let us
   set those here. For convenience (because the prelude provides it),
   we will use truncation towards zero. *)
infixl ( / ) div
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod

(* We will be using linear lists. Define a convenient notation. *)
#define NIL list_vt_nil ()
#define ::  list_vt_cons

(*------------------------------------------------------------------*)
(* Something we will use: copy the contents of one arrayptr to
   another arrayptr. *)

extern fn {a : t@ype}
arrayptr_copy_over
          {n : int}
          (n   : int n,
           src : !arrayptr (a, n),
           dst : !arrayptr (a, n))
    : void

implement {a}
arrayptr_copy_over {n} (n, src, dst) =
  let
    fun
    loop (i   : intGte 0,
          src : !arrayptr (a, n),
          dst : !arrayptr (a, n))
        : void =
      if i < n then
        begin
          dst[i] := src[i];
          loop (succ i, src, dst)
        end
  in
    loop (0, src, dst)
  end

overload copy_over with arrayptr_copy_over

(*------------------------------------------------------------------*)
(* The basics. *)

(* For this pedagogical example, let us choose a particular integer
   type, thus avoiding the clutter of template notation. *)
typedef integer = int

(* A generator is a recursive type that forms a tree. *)
datavtype generator_vt =
| generator_vt_nil of ()        (* The nil generator. *)
| {n : int}
  generator_vt_cons of          (* A non-nil generator. *)
    (generator_func_vt n,       (* What does the work. *)
     int n,                     (* The size of workspace. *)
     arrayptr (integer, n),     (* The initial value of workspace. *)
     arrayptr (integer, n),     (* The workspace. *)
     List_vt generator_vt)      (* The sources. *)
where generator_func_vt (n : int) =
  (int n,                         (* The size of workspace. *)
   !arrayptr (integer, n),        (* The workspace. *)
   !List_vt generator_vt,         (* The sources. *)
   &bool? >> bool b,              (* Is there a term? *)
   &integer? >> opt (integer, b)) (* The term, if any. *)
    -> #[b : bool] void

(* Reinitializes a generator. (Needed because there is no
   memoization.) *)
extern fn generator_vt_initialize : (&generator_vt) -> void
overload initialize with generator_vt_initialize

(* Frees a generator. *)
extern fn generator_vt_free : generator_vt -> void
overload free with generator_vt_free

(*------------------------------------------------------------------*)
(* A function to print the output of a generator as Plain TeX. *)

extern fn
fprinttex_generator_output
          (outf      : FILEref,
           gen       : &generator_vt,
           max_terms : intGte 1)
    : void

(*------------------------------------------------------------------*)
(* Some functions to make generators. *)

extern fn                       (* For a rational number. *)
r2cf_make (n : integer,
           d : integer)
    : generator_vt

extern fn                       (* For the square root of 2. *)
sqrt2_make ()
    : generator_vt

extern fn                       (* For a homographic function. *)
hfunc_make (a1      : integer,
            a       : integer,
            b1      : integer,
            b       : integer,
            sources : List1_vt generator_vt)
    : generator_vt

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

implement
generator_vt_initialize gen =
  let
    fun
    recurs (gen : &generator_vt) : void =
      case+ gen of
      | generator_vt_nil () => ()
      | @ generator_vt_cons (_, worksize, initial, workspace,
                             sources) =>
        let
          fun
          initialize_recursively
                    (p : !List_vt generator_vt)
              : void =
            case+ p of
            | NIL => ()
            | @ (gen :: tail) =>
              begin
                recurs gen;
                initialize_recursively tail;
                fold@ p
              end
        in
          copy_over (worksize, initial, workspace);
          initialize_recursively sources;
          fold@ gen
        end
  in
    recurs gen
  end

implement
generator_vt_free gen =
  let
    fun
    recurs (gen : generator_vt) : void =
      case+ gen of
      | ~ generator_vt_nil () => ()
      | ~ generator_vt_cons (_, _, initial, workspace, sources) =>
        begin
          free initial;
          free workspace;
          list_vt_freelin_fun (sources, lam x =<fun1> recurs x)
        end
  in
    recurs gen
  end

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

implement
fprinttex_generator_output (outf, gen, max_terms) =
  let
    fun
    loop (gen         : &generator_vt >> _,
          sep         : int,
          terms_count : intGte 0)
        : void =
      if terms_count = max_terms then
        fprint! (outf, ",\\cdots\\,]")
      else
        let
          var term_exists : bool?
          var term : integer?
        in
          case+ gen of
          | generator_vt_nil () => ()
          | @ generator_vt_cons (run, worksize, _, workspace,
                                 sources) =>
            let
              var term_exists : bool?
              var term : integer?
            in
              run (worksize, workspace, sources, term_exists, term);
              if term_exists then
                let
                  prval () = opt_unsome term
                  prval () = fold@ gen
                in
                  case+ sep of
                  | 0 => fprint! (outf, "[\\,")
                  | 1 => fprint! (outf, ";")
                  | _ => fprint! (outf, ",");
                  fprint! (outf, term);
                  loop (gen, if sep = 0 then 1 else 2,
                        succ terms_count)
                end
              else
                let
                  prval () = opt_unnone term
                  prval () = fold@ gen
                in
                  fprint! (outf, "\\,]")              
                end
            end
        end
  in
    initialize gen;
    loop (gen, 0, 0);
    initialize gen
  end

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

fn
r2cf_run : generator_func_vt 2 =
  lam (worksize, workspace, _sources, term_exists, term) =>
    let
      prval () = lemma_arrayptr_param workspace
      val () = assertloc (2 <= worksize)
      val d = arrayptr_get_at<integer> (workspace, 1)
    in
      if d = 0 then
        begin
          term_exists := false;
          {prval () = opt_none term}
        end
      else
        let
          val n = workspace[0]
          val @(q, r) = @(n div d, n rem d)
        in
          workspace[0] := d;
          workspace[1] := r;
          term_exists := true;
          term := q;
          {prval () = opt_some term}
        end
    end

implement
r2cf_make (n, d) =
  let
    val workspace = arrayptr_make_elt (i2sz 2, 0)
    val initial = arrayptr_make_elt (i2sz 2, 0)
    val () = initial[0] := n
    and () = initial[1] := d
  in
    copy_over (2, initial, workspace);
    generator_vt_cons (r2cf_run, 2, initial, workspace, NIL)
  end

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

fn
sqrt2_run : generator_func_vt 1 =
  lam (worksize, workspace, _sources, term_exists, term) =>
    let
      prval () = lemma_arrayptr_param workspace
      val () = assertloc (1 <= worksize)
    in
      term_exists := true;
      term := arrayptr_get_at<integer> (workspace, 0);
      {prval () = opt_some term};
      arrayptr_set_at<integer> (workspace, 0, 2)
    end

implement
sqrt2_make () =
  let
    val workspace = arrayptr_make_elt (i2sz 1, 0)
    val initial = arrayptr_make_elt (i2sz 1, 0)
    val () = initial[0] := 1
  in
    copy_over (1, initial, workspace);
    generator_vt_cons (sqrt2_run, 1, initial, workspace, NIL)
  end

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

fn
hfunc_run : generator_func_vt 4 =
  lam (worksize, workspace, sources, term_exists, term) =>
    let
      prval () = lemma_arrayptr_param workspace
      val () = assertloc (4 <= worksize)

      fun
      loop (workspace   : !arrayptr (integer, 4),
            sources     : !List_vt generator_vt,
            term_exists : &bool? >> bool b,
            term        : &integer? >> opt (integer, b))
          : #[b : bool] void =
        let
          val b1 = arrayptr_get_at<integer> (workspace, 2)
          and b = arrayptr_get_at<integer> (workspace, 3)
        in
          if b1 = 0 && b = 0 then
            begin
              term_exists := false;
              {prval () = opt_none term}
            end
          else
            let
              val a1 = workspace[0]
              and a = workspace[1]

              fn
              take_term (workspace : !arrayptr (integer, 4),
                         sources   : !List_vt generator_vt)
                  : void =
                let
                  val- @ (source :: _) = sources
                  val- @ generator_vt_cons
                          (run1, worksize1, _, workspace1,
                           sources1) = source

                  var term_exists1 : bool?
                  var term1 : integer?
                in
                  run1 (worksize1, workspace1, sources1,
                        term_exists1, term1);
                  if term_exists1 then
                    let
                      prval () = opt_unsome term1
                    in
                      workspace[0] := a + (a1 * term1);
                      workspace[1] := a1;
                      workspace[2] := b + (b1 * term1);
                      workspace[3] := b1;
                      fold@ source;
                      fold@ sources
                    end
                  else
                    let
                      prval () = opt_unnone term1
                    in
                      workspace[1] := a1;
                      workspace[3] := b1;
                      fold@ source;
                      fold@ sources
                    end
                end
            in
              if b1 <> 0 && b <> 0 then
                let
                  val q1 = a1 div b1
                  and q = a div b
                in
                  if q1 = q then
                    begin
                      workspace[0] := b1;
                      workspace[1] := b;
                      workspace[2] := a1 - (b1 * q);
                      workspace[3] := a - (b * q);
                      term_exists := true;
                      term := q;
                      {prval () = opt_some term}
                    end
                  else
                    begin
                      take_term (workspace, sources);
                      loop (workspace, sources, term_exists, term)
                    end
                end
              else
                begin
                  take_term (workspace, sources);
                  loop (workspace, sources, term_exists, term)
                end
            end
        end
    in
      loop (workspace, sources, term_exists, term)
    end

implement
hfunc_make (a1, a, b1, b, sources) =
  let
    val workspace = arrayptr_make_elt (i2sz 4, 0)
    val initial = arrayptr_make_elt (i2sz 4, 0)
    val () = initial[0] := a1
    val () = initial[1] := a
    val () = initial[2] := b1
    val () = initial[3] := b
  in
    copy_over (4, initial, workspace);
    generator_vt_cons (hfunc_run, 4, initial, workspace, sources)
  end

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

#define MAX_TERMS 20
#define GOES_TO "&\\rightarrow "
#define END_LINE "\\cr\n"

fn
fprinttex_rational_number
          (outf : FILEref,
           n    : integer,
           d    : integer)
    : void =
  let
    var gen = r2cf_make (n, d)
  in
    fprint! (outf, n, "\\over ", d, GOES_TO);
    fprinttex_generator_output (outf, gen, MAX_TERMS);
    fprint! (outf, END_LINE);
    free gen
  end

fn
fprinttex_sqrt2
          (outf : FILEref)
    : void =
  let
    var gen = sqrt2_make ()
  in
    fprint! (outf, "\\sqrt 2", GOES_TO);
    fprinttex_generator_output (outf, gen, MAX_TERMS);
    fprint! (outf, END_LINE);
    free gen
  end

fn
fprinttex_hfunc_of_rational_number
          (outf : FILEref,
           expr : string,
           a1   : integer,
           a    : integer,
           b1   : integer,
           b    : integer,
           n    : integer,
           d    : integer)
    : void =
  let
    var gen = hfunc_make (a1, a, b1, b, r2cf_make (n, d) :: NIL)
  in
    fprint! (outf, expr, GOES_TO);
    fprinttex_generator_output (outf, gen, MAX_TERMS);
    fprint! (outf, END_LINE);
    free gen
  end

fn
fprinttex_hfunc_of_sqrt2
          (outf : FILEref,
           expr : string,
           a1   : integer,
           a    : integer,
           b1   : integer,
           b    : integer)
    : void =
  let
    var gen = hfunc_make (a1, a, b1, b, sqrt2_make () :: NIL)
  in
    fprint! (outf, expr, GOES_TO);
    fprinttex_generator_output (outf, gen, MAX_TERMS);
    fprint! (outf, END_LINE);
    free gen
  end

fn
fprinttex_complicated
          (outf : FILEref)
    : void =
  (* Here we demonstrate a more complicated nesting of generators. *)
  let
    (* gen1 = 1/sqrt(2) *)
    val gen1 = hfunc_make (0, 1, 1, 0, sqrt2_make () :: NIL)

    (* gen2 = 1 + gen1 *)
    val gen2 = hfunc_make (1, 1, 0, 1, gen1 :: NIL)

    (* gen = gen2 / 2 *)
    var gen = hfunc_make (1, 0, 0, 2, gen2 :: NIL)
  in
    fprint! (outf, "{1 + {1\\over\\sqrt 2}}\\over 2", GOES_TO);
    fprinttex_generator_output (outf, gen, MAX_TERMS);
    fprint! (outf, END_LINE);
    free gen
  end

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

fn
fprint_14point (outf : FILEref) : void =
  begin
    fprintln! (outf, "%%% This file is public domain.");
    fprintln! (outf, "%%% Originally written 1992, Don Hosek.");
    fprintln! (outf, "%%% This declaration added by Clea F. Rees 2008/11/16 with the permission of Dan Hosek.");
    fprintln! (outf, "%%%");
    fprintln! (outf, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
    fprintln! (outf, "% This file sets up a fourteen point environment for TeX. It can be initialized");
    fprintln! (outf, "% with the '\\fourteenpoint' macro.");
    fprintln! (outf, "% It also sets up a '\\tenpoint' macro in case you want to go back down.");
    fprintln! (outf, "% By Don Hosek");
    fprintln! (outf, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
    fprintln! (outf, " ");
    fprintln! (outf, "\\ifx\\tenpoint\\undefined\\let\\loadedfrommacro=Y");
    fprintln! (outf, "         \\input 10point");
    fprintln! (outf, "         \\let\\loadedfrommacro=N\\fi");
    fprintln! (outf, " ");
    fprintln! (outf, "%%%");
    fprintln! (outf, "%%% Load in the fonts");
    fprintln! (outf, "%%%");
    fprintln! (outf, "\\font\\fourteenrm=cmr12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteeni=cmmi12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteensy=cmsy10 scaled \\magstep2");
    fprintln! (outf, "\\font\\fourteenex=cmex10 scaled \\magstep2");
    fprintln! (outf, "\\font\\fourteenbf=cmbx12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteensl=cmsl12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteentt=cmtt12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteenit=cmti12 scaled \\magstep1");
    fprintln! (outf, "\\font\\fourteencsc=cmcsc10 scaled \\magstep2");
    fprintln! (outf, " ");
    fprintln! (outf, "%%%");
    fprintln! (outf, "%%% Set up the fourteenpoint macros");
    fprintln! (outf, "%%%");
    fprintln! (outf, "\\ifx\\fourteenpoint\\undefined");
    fprintln! (outf, "   \\def\\fourteenpoint{\\def\\rm{\\fam0\\fourteenrm}% switch to 14-point type");
    fprintln! (outf, "       \\textfont0=\\fourteenrm \\scriptfont0=\\tenrm \\scriptscriptfont0=\\sevenrm");
    fprintln! (outf, "       \\textfont1=\\fourteeni  \\scriptfont1=\\teni  \\scriptscriptfont1=\\seveni");
    fprintln! (outf, "       \\textfont2=\\fourteensy \\scriptfont2=\\tensy \\scriptscriptfont2=\\sevensy");
    fprintln! (outf, "       \\textfont3=\\fourteenex \\scriptfont3=\\fourteenex");
    fprintln! (outf, "                              \\scriptscriptfont3=\\fourteenex");
    fprintln! (outf, "       \\textfont\\itfam=\\fourteenit  \\def\\it{\\fam\\itfam\\fourteenit}%");
    fprintln! (outf, "       \\textfont\\slfam=\\fourteensl  \\def\\sl{\\fam\\slfam\\fourteensl}%");
    fprintln! (outf, "       \\textfont\\ttfam=\\fourteentt  \\def\\tt{\\fam\\ttfam\\fourteentt}%");
    fprintln! (outf, "       \\textfont\\bffam=\\fourteenbf  \\scriptfont\\bffam=\\tenbf");
    fprintln! (outf, "        \\scriptscriptfont\\bffam=\\sevenbf  \\def\\bf{\\fam\\bffam\\fourteenbf}%");
    fprintln! (outf, "       \\textfont\\scfam=\\fourteencsc \\def\\sc{\\fam\\scfam\\fourteencsc}%");
    fprintln! (outf, "       \\normalbaselineskip=17pt");
    fprintln! (outf, "       \\setbox\\strutbox=\\hbox{\\vrule height11.9pt depth6.3pt width0pt}%");
    fprintln! (outf, "       \\normalbaselines\\rm}");
    fprintln! (outf, "   \\fi")
  end

implement
main () =
  let
    val outf = stdout_ref
  in
    (* I assume the TeX processor is LuaTeX. *)
    fprintln! (outf, "\\pagewidth 6in\\hoffset-1in\\hsize 6in");
    fprintln! (outf, "\\pageheight 6in\\voffset-1.05in\\vsize 6in");

    (* Suppress the page number. *)
    fprintln! (outf, "\\footline={}");

    (* Print large. *)
    fprint_14point (outf);
    fprintln! (outf, "\\fourteenpoint");

    fprintln! (outf, "\\normallineskip 6pt");
    fprintln! (outf, "\\normalbaselines");

    fprintln! (outf, "$$\\eqalign{");

    fprinttex_rational_number (outf, 13, 11);
    fprinttex_rational_number (outf, 22, 7);
    fprinttex_sqrt2 (outf);

    fprinttex_hfunc_of_rational_number
      (outf, "{13\\over 11} + {1\\over 2}", 2, 1, 0, 2, 13, 11);
    fprinttex_hfunc_of_rational_number
      (outf, "{22\\over 7} + {1\\over 2}", 2, 1, 0, 2, 22, 7);
    fprinttex_hfunc_of_rational_number
      (outf, "{22\\over 7}\\over 4", 1, 0, 0, 4, 22, 7);
    fprinttex_hfunc_of_sqrt2
      (outf, "{\\sqrt 2}\\over 2", 1, 0, 0, 2);
    fprinttex_hfunc_of_sqrt2
      (outf, "1\\over\\sqrt 2", 0, 1, 1, 0);
    fprinttex_hfunc_of_sqrt2
      (outf, "{2 + \\sqrt 2}\\over 4", 1, 2, 0, 4);
    fprinttex_complicated outf;

    fprintln! (outf, "}$$");
    fprintln! (outf, "\\bye");
    0
  end

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

Run something like:

$ patscc -DATS_MEMALLOC_LIBC -O2 -std=gnu2x univariate-continued-fraction-task-no-gc.dats
$ ./a.out > foo.tex
$ luatex foo.tex
$ okular foo.pdf # Or your favorite PDF viewer.
Output:

The output of the program.

Using lazy non-linear types

Translation of: Haskell
Translation of: Mercury

This method is simple, and it memoizes terms. However, the memoization is in a linked list rather than a randomly accessible array.

The recurs routines do not execute recursions, but instead (thanks to $delay) create what I will call "recursive thunks". Otherwise the stack would overflow.

The code leaks memory, so using a garbage collector may be a good idea.

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

#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"

(* For convenience (because the prelude provides it), we will use
   integer division with truncation towards zero. *)
infixl ( / ) div
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod

(*------------------------------------------------------------------*)
(* The definition of a continued fraction, and a few simple ones. *)

typedef cf (tk : tkind) = stream (g0int tk)

(* A "continued fraction" with no terms. *)
fn {tk : tkind}
cfnil ()
    : cf tk =
  stream_make_nil<g0int tk> ()

(* A continued fraction of one term followed by more terms. *)
fn {tk : tkind}
cfcons (term : g0int tk,
        more : cf tk)
    : cf tk =
  stream_make_cons<g0int tk> (term, more)

(* A continued fraction with all terms equal. *)
fn {tk : tkind}
repeat_forever (term : g0int tk)
    : cf tk =
  let
    fun recurs () : stream_con (g0int tk) =
      stream_cons (term, $delay recurs ())
  in
    $delay recurs ()
  end

(* The square root of two. *)
fn {tk : tkind}
sqrt2 ()
    : cf tk =
  cfcons<tk> (g0i2i 1, repeat_forever<tk> (g0i2i 2))

(*------------------------------------------------------------------*)
(* A continued fraction for a rational number. *)

typedef ratnum (tk : tkind) = @(g0int tk, g0int tk)

fn {tk : tkind}
r2cf_integers (n : g0int tk,
               d : g0int tk)
    : cf tk =
  let
    fun recurs (n : g0int tk,
                d : g0int tk)
        : cf tk =
      if iseqz d then
        cfnil<tk> ()
      else
        cfcons<tk> (n div d, recurs (d, n rem d))
  in
    recurs (n, d)
  end

fn {tk : tkind}
r2cf_ratnum (r : ratnum tk)
    : cf tk =
  r2cf_integers (r.0, r.1)

overload r2cf with r2cf_integers
overload r2cf with r2cf_ratnum

(*------------------------------------------------------------------*)
(* Application of a homographic function to a continued fraction. *)

typedef ng4 (tk : tkind) = @(g0int tk, g0int tk,
                             g0int tk, g0int tk)

fn {tk : tkind}
apply_ng4 (ng4      : ng4 tk,
           other_cf : cf tk)
    : cf tk =
  let
    typedef t = g0int tk

    fun
    recurs (a1       : t,
            a        : t,
            b1       : t,
            b        : t,
            other_cf : cf tk)
        : stream_con t =
      let
        fn {}
        eject_term (a1       : t,
                    a        : t,
                    b1       : t,
                    b        : t,
                    other_cf : cf tk,
                    term     : t)
            : stream_con t =
          stream_cons (term, $delay recurs (b1, b, a1 - (b1 * term),
                                            a - (b * term), other_cf))

        fn {}
        absorb_term (a1       : t,
                     a        : t,
                     b1       : t,
                     b        : t,
                     other_cf : cf tk)
            : stream_con t =
          case+ !other_cf of
          | stream_nil () =>
            recurs (a1, a1, b1, b1, other_cf)
          | stream_cons (term, rest) =>
            recurs (a + (a1 * term), a1, b + (b1 * term), b1, rest)
      in
        if iseqz b1 && iseqz b then
          stream_nil ()
        else if iseqz b1 || iseqz b then
          absorb_term (a1, a, b1, b, other_cf)
        else
          let
            val q1 = a1 div b1
            and q = a div b
          in
            if q1 = q then
              eject_term (a1, a, b1, b, other_cf, q)
            else
              absorb_term (a1, a, b1, b, other_cf)
          end
      end

    val @(a1, a, b1, b) = ng4
  in
    $delay recurs (a1, a, b1, b, other_cf)
  end

(*------------------------------------------------------------------*)
(* Some special cases of homographic functions. *)

fn {tk : tkind}
integer_add_cf (n  : g0int tk,
                cf : cf tk)
    : cf tk =
  apply_ng4 (@(g0i2i 1, n, g0i2i 0, g0i2i 1), cf)

fn {tk : tkind}
cf_add_ratnum (cf : cf tk,
               r  : ratnum tk)
    : cf tk =
  let
    val @(n, d) = r
  in
    apply_ng4 (@(d, n, g0i2i 0, d), cf)
  end

fn {tk : tkind}
cf_mul_ratnum (cf : cf tk,
               r  : ratnum tk)
    : cf tk =
  let
    val @(n, d) = r
  in
    apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), cf)
  end

fn {tk : tkind}
cf_div_integer (cf : cf tk,
                n  : g0int tk)
    : cf tk =
  apply_ng4 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i n), cf)

fn {tk : tkind}
integer_div_cf (n  : g0int tk,
                cf : cf tk)
    : cf tk =
  apply_ng4 (@(g0i2i 0, g0i2i n, g0i2i 1, g0i2i 0), cf)

overload + with integer_add_cf
overload + with cf_add_ratnum
overload * with cf_mul_ratnum
overload / with cf_div_integer
overload / with integer_div_cf

(*------------------------------------------------------------------*)
(* cf2string: convert a continued fraction to a string. *)

fn {tk : tkind}
cf2string_max_terms_given
          (cf        : cf tk,
           max_terms : intGte 1)
    : string =
  let
    fun
    loop (i     : intGte 0,
          cf    : cf tk,
          accum : List0 string)
        : List0 string =
      case+ !cf of
      | stream_nil () => list_cons ("]", accum)
      | stream_cons (term, rest) =>
        if i = max_terms then
          list_cons (",...]", accum)
        else
          let
            val accum =
              list_cons
                (tostring_val<g0int tk> term,
                 (case+ i of
                  | 0 => accum
                  | 1 => list_cons (";", accum)
                  | _ => list_cons (",", accum)) : List0 string)
          in
            loop (succ i, rest, accum)
          end

    val string_lst = list_vt2t (reverse (loop (0, cf, list_sing "[")))
  in
    strptr2string (stringlst_concat string_lst)
  end

extern fn {tk : tkind}
cf2string$max_terms :
  () -> intGte 1

implement {tk} cf2string$max_terms () = 20

fn {tk : tkind}
cf2string_max_terms_default
          (cf : cf tk)
    : string =
  cf2string_max_terms_given<tk> (cf, cf2string$max_terms<tk> ())

overload cf2string with cf2string_max_terms_given
overload cf2string with cf2string_max_terms_default

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

fn {tk : tkind}
show (expression : string,
      cf         : cf tk)
    : void =
  begin
    print! expression;
    print! " => ";
    println! (cf2string<tk> cf);
  end

implement
main () =
  let
    val cf_13_11 = r2cf (13, 11)
    val cf_22_7 = r2cf (22, 7)
    val cf_sqrt2 = sqrt2<intknd> ()
    val cf_1_sqrt2 = 1 / cf_sqrt2
  in
    show ("13/11", cf_13_11);
    show ("22/7", cf_22_7);
    show ("sqrt(2)", cf_sqrt2);
    show ("13/11 + 1/2", cf_13_11 + @(1, 2));
    show ("22/7 + 1/2", cf_22_7 + @(1, 2));
    show ("(22/7)/4", cf_22_7 * @(1, 4));
    show ("1/sqrt(2)", cf_1_sqrt2);
    show ("(2 + sqrt(2))/4", apply_ng4 (@(1, 2, 0, 4), cf_sqrt2));

    (* To show it can be done, write the following without using
       results already obtained: *)
    show ("(1 + 1/sqrt(2))/2", (1 + 1/sqrt2())/2);

    0
  end

(*------------------------------------------------------------------*)
Output:
$ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW univariate-continued-fraction-task-lazy.dats -lgc && ./a.out
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

C

With a garbage collector

Translation of: ATS
Library: Boehm GC

I use a garbage collector to make the free use of heap space both easier and more error-free. Or, in many cases, you can simply let the memory leak.

The output of a generator is memoized and can be reused without computing again. If the generator needs to produce more terms, it picks up where it left off.

/*------------------------------------------------------------------*/
/* For C with Boehm GC as garbage collector. */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <gc.h>                 /* Boehm GC. */

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

/* Let us choose an integer type. */
typedef long long int integer;

/* We need consistent definitions of division and remainder. Let us
   set those here. For convenience (because C provides it), we will
   use truncation towards zero. */
#define DIV(a, b) ((a) / (b))
#define REM(a, b) ((a) % (b))

/* Choose a memory allocator: Boehm GC. (Ideally one should check for
   NULL return values, but for this pedagogical example let us skip
   that.) */
#define MALLOC_INIT() GC_INIT ()
#define MALLOC GC_MALLOC
#define REALLOC GC_REALLOC
#define FREE GC_FREE            /* Or one could make this a no-op. */

/*------------------------------------------------------------------*/
/* Some operations on char-strings. (In practice, I would write an m4
   macro to create such repetitive C functions for me. Of course, it
   is also possible to use <stdarg.h> or some such [generally unsafe]
   mechanism.) */

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

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

static 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 = MALLOC (n1 + n2 + n3 + 1);
  s[n1 + n2 + n3] = '\0';
  memcpy (s, s1, n1);
  memcpy (s + n1, s2, n2);
  memcpy (s + n1 + n2, s3, n3);
  return s;
}

/*------------------------------------------------------------------*/
/* Continued fractions as processes for generating terms. The terms
   are memoized and are accessed by their zero-based index. */

typedef void cf_generator_func_t (void *env, bool *there_is_a_term,
                                  integer *term);

struct _cf_generator  /* For practical purposes, this is a closure. */
{
  cf_generator_func_t *func;
  void *env;
};

typedef struct _cf_generator *cf_generator_t;

struct _cf
{
  bool terminated;          /* No more terms? */
  size_t m;                 /* The number of terms computed so far. */
  size_t n;                 /* The size of memo storage. */
  integer *memo;            /* Memoized terms. */
  cf_generator_t gen;       /* A closure to generate terms. */
};

typedef struct _cf *cf_t;

cf_generator_t
cf_generator_make (cf_generator_func_t func, void *env)
{
  cf_generator_t gen = MALLOC (sizeof (struct _cf_generator));
  gen->func = func;
  gen->env = env;
  return gen;
}

cf_t
cf_make (cf_generator_t gen)
{
  const size_t start_size = 8;
  integer *memo = MALLOC (start_size * sizeof (integer));
  cf_t cf = MALLOC (sizeof (struct _cf));
  cf->terminated = false;
  cf->m = 0;
  cf->n = start_size;
  cf->memo = memo;
  cf->gen = gen;
  return cf;
}

static void
_cf_get_more_terms (cf_t cf, size_t needed)
{
  size_t term_count = cf->m;
  bool done = false;
  while (!done)
    {
      if (term_count == needed)
        {
          cf->m = needed;
          done = true;
        }
      else
        {
          bool there_is_a_term;
          integer term;
          cf->gen->func (cf->gen->env, &there_is_a_term, &term);
          if (there_is_a_term)
            {
              cf->memo[term_count] = term;
              term_count += 1;
            }
          else
            {
              cf->terminated = true;
              cf->m = term_count;
              done = true;
            }
        }
    }
}

static void
_cf_update (cf_t cf, size_t needed)
{
  if (cf->terminated || needed <= (cf->m))
    /* Do nothing. */ ;
  else if (needed <= (cf->n))
    _cf_get_more_terms (cf, needed);
  else
    {
      /* Provide twice the needed storage. */
      cf->n = 2 * needed;
      cf->memo = REALLOC (cf->memo, cf->n * sizeof (integer));
      _cf_get_more_terms (cf, needed);
    }
}

void
cf_get_at (cf_t cf, size_t i, bool *there_is_a_term,
           integer *term)
{
  _cf_update (cf, i + 1);
  *there_is_a_term = (i < (cf->m));
  if (*there_is_a_term)
    *term = cf->memo[i];
}

char *
cf2string_max_terms (cf_t cf, size_t max_terms)
{
  char *s = string_append1 ("[");
  const char *sep = "";
  size_t i = 0;
  bool done = false;
  while (!done)
    {
      if (i == max_terms)
        {
          s = string_append2 (s, ",...]");
          done = true;
        }
      else
        {
          bool there_is_a_term;
          integer term;
          cf_get_at (cf, i, &there_is_a_term, &term);
          if (there_is_a_term)
            {
              char buf1[1];
              const int numeral_len =
                snprintf (buf1, 1, "%jd", (intmax_t) term);
              char buf[numeral_len + 1];
              snprintf (buf, numeral_len + 1, "%jd", (intmax_t) term);
              s = string_append3 (s, sep, buf);
              sep = (sep[0] == '\0') ? ";" : ",";
              i += 1;
            }
          else
            {
              s = string_append2 (s, "]");
              done = true;
            }
        }
    }
  return s;
}

char *
cf2string (cf_t cf)
{
  const size_t default_max_terms = 20;
  return cf2string_max_terms (cf, default_max_terms);
}

/*------------------------------------------------------------------*/
/* Using a cf_t as a cf_generator_t. */

struct _cf_gen_env
{
  cf_t cf;
  size_t i;
};

static void
cf_gen_run (void *env, bool *there_is_a_term, integer *term)
{
  struct _cf_gen_env *state = env;
  cf_get_at (state->cf, state->i, there_is_a_term, term);
  state->i += 1;
}

cf_generator_t
cf_gen_make (cf_t cf)
{
  struct _cf_gen_env *state = MALLOC (sizeof (struct _cf_gen_env));
  state->cf = cf;
  state->i = 0;
  return cf_generator_make (cf_gen_run, state);
}

/*------------------------------------------------------------------*/
/* A homographic function. */

struct _hfunc
{
  integer a1;
  integer a;
  integer b1;
  integer b;
};

typedef struct _hfunc *hfunc_t;

struct _hfunc_gen_env
{
  struct _hfunc state;
  cf_generator_t gen;
};

hfunc_t
hfunc_make (integer a1, integer a, integer b1, integer b)
{
  hfunc_t f = MALLOC (sizeof (struct _hfunc));
  f->a1 = a1;
  f->a = a;
  f->b1 = b1;
  f->b = b;
  return f;
}

static void
_take_term_from_ngen (struct _hfunc *state,
                      cf_generator_t ngen)
{
  const integer a1 = state->a1;
  const integer b1 = state->b1;

  bool there_is_a_term;
  integer term;
  ngen->func (ngen->env, &there_is_a_term, &term);
  if (there_is_a_term)
    {
      const integer a = state->a;
      const integer b = state->b;

      state->a1 = a + (a1 * term);
      state->a = a1;
      state->b1 = b + (b1 * term);
      state->b = b1;
    }
  else
    {
      state->a = a1;
      state->b = b1;
    }
}

static void
_adjust_state_for_term_output (struct _hfunc *state,
                               integer term)
{
  const integer a1 = state->a1;
  const integer a = state->a;
  const integer b1 = state->b1;
  const integer b = state->b;

  state->a1 = b1;
  state->a = b;
  state->b1 = a1 - (b1 * term);
  state->b = a - (b * term);
}

static void
hfunc_gen_run (void *env, bool *there_is_a_term, integer *term)
{
  struct _hfunc *state = &(((struct _hfunc_gen_env *) env)->state);
  cf_generator_t ngen = ((struct _hfunc_gen_env *) env)->gen;

  bool done = false;
  while (!done)
    {
      const bool b1_iseqz = (state->b1 == 0);
      const bool b_iseqz = (state->b == 0);
      if (b1_iseqz && b_iseqz)
        {
          *there_is_a_term = false;
          done = true;
        }
      else if (!b1_iseqz && !b_iseqz)
        {
          const integer q1 = DIV (state->a1, state->b1);
          const integer q = DIV (state->a, state->b);
          if (q1 == q)
            {
              _adjust_state_for_term_output (state, q);
              *there_is_a_term = true;
              *term = q;
              done = true;
            }
          else
            _take_term_from_ngen (state, ngen);
        }
      else
        _take_term_from_ngen (state, ngen);
    }
}

/* Make a new generator that applies an hfunc_t to another
   generator. */
cf_generator_t
hfunc_gen_make (hfunc_t f, cf_generator_t gen)
{
  struct _hfunc_gen_env *env =
    MALLOC (sizeof (struct _hfunc_gen_env));
  env->state = *f;
  env->gen = gen;
  return cf_generator_make (hfunc_gen_run, env);
}

/* Make a new cf_t that applies an hfunc_t to another cf_t. */
cf_t
hfunc_apply (hfunc_t f, cf_t cf)
{
  cf_generator_t gen1 = cf_gen_make (cf);
  cf_generator_t gen2 = hfunc_gen_make (f, gen1);
  return cf_make (gen2);
}

/*------------------------------------------------------------------*/
/* Creation of a cf_t for a rational number. */

struct _r2cf_gen_env
{
  integer n;
  integer d;
};

static void
r2cf_gen_run (void *env, bool *there_is_a_term, integer *term)
{
  struct _r2cf_gen_env *state = env;
  *there_is_a_term = (state->d != 0);
  if (*there_is_a_term)
    {
      const integer q = DIV (state->n, state->d);
      const integer r = REM (state->n, state->d);
      state->n = state->d;
      state->d = r;
      *term = q;
    }
}

cf_generator_t
r2cf_gen_make (integer n, integer d)
{
  struct _r2cf_gen_env *state =
    MALLOC (sizeof (struct _r2cf_gen_env));
  state->n = n;
  state->d = d;
  return cf_generator_make (r2cf_gen_run, state);
}

cf_t
r2cf (integer n, integer d)
{
  return cf_make (r2cf_gen_make (n, d));
}

/*------------------------------------------------------------------*/
/* Creation of a cf_t for sqrt(2). */

struct _sqrt2_gen_env
{
  integer term;
};

static void
sqrt2_gen_run (void *env, bool *there_is_a_term, integer *term)
{
  struct _sqrt2_gen_env *state = env;
  *there_is_a_term = true;
  *term = state->term;
  state->term = 2;
}

cf_generator_t
sqrt2_gen_make (void)
{
  struct _sqrt2_gen_env *state =
    MALLOC (sizeof (struct _sqrt2_gen_env));
  state->term = 1;
  return cf_generator_make (sqrt2_gen_run, state);
}

cf_t
sqrt2_make (void)
{
  return cf_make (sqrt2_gen_make ());
}

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

int
main (void)
{
  MALLOC_INIT ();

  hfunc_t add_one_half = hfunc_make (2, 1, 0, 2);
  hfunc_t add_one = hfunc_make (1, 1, 0, 1);
  hfunc_t divide_by_two = hfunc_make (1, 0, 0, 2);
  hfunc_t divide_by_four = hfunc_make (1, 0, 0, 4);
  hfunc_t take_reciprocal = hfunc_make (0, 1, 1, 0);
  hfunc_t add_one_then_div_two = hfunc_make (1, 1, 0, 2);
  hfunc_t add_two_then_div_four = hfunc_make (1, 2, 0, 4);

  cf_t cf_13_11 = r2cf (13, 11);
  cf_t cf_22_7 = r2cf (22, 7);
  cf_t cf_sqrt2 = sqrt2_make ();

  cf_t cf_13_11_plus_1_2 = hfunc_apply (add_one_half, cf_13_11);
  cf_t cf_22_7_plus_1_2 = hfunc_apply (add_one_half, cf_22_7);
  cf_t cf_22_7_div_4 = hfunc_apply (divide_by_four, cf_22_7);

  /* The following two give the same result: */
  cf_t cf_sqrt2_div_2 = hfunc_apply (divide_by_two, cf_sqrt2);
  cf_t cf_1_div_sqrt2 = hfunc_apply (take_reciprocal, cf_sqrt2);
  assert (strcmp (cf2string (cf_sqrt2_div_2),
                  cf2string (cf_1_div_sqrt2)) == 0);

  /* The following three give the same result: */
  cf_t cf_2_plus_sqrt2_grouped_div_4 =
    hfunc_apply (add_two_then_div_four, cf_sqrt2);
  cf_t cf_half_of_1_plus_half_sqrt2 =
    hfunc_apply (add_one_then_div_two, cf_sqrt2_div_2);
  cf_t cf_half_of_1_plus_1_div_sqrt2 =
    hfunc_apply (divide_by_two,
                 hfunc_apply (add_one, cf_sqrt2_div_2));
  assert (strcmp (cf2string (cf_2_plus_sqrt2_grouped_div_4),
                  cf2string (cf_half_of_1_plus_half_sqrt2)) == 0);
  assert (strcmp (cf2string (cf_half_of_1_plus_half_sqrt2),
                  cf2string (cf_half_of_1_plus_1_div_sqrt2)) == 0);

  printf ("13/11 => %s\n", cf2string (cf_13_11));
  printf ("22/7 => %s\n", cf2string (cf_22_7));
  printf ("sqrt(2) => %s\n", cf2string (cf_sqrt2));
  printf ("13/11 + 1/2 => %s\n", cf2string (cf_13_11_plus_1_2));
  printf ("22/7 + 1/2 => %s\n", cf2string (cf_22_7_plus_1_2));
  printf ("(22/7)/4 => %s\n", cf2string (cf_22_7_div_4));
  printf ("sqrt(2)/2 => %s\n", cf2string (cf_sqrt2_div_2));
  printf ("1/sqrt(2) => %s\n", cf2string (cf_1_div_sqrt2));
  printf ("(2+sqrt(2))/4 => %s\n",
          cf2string (cf_2_plus_sqrt2_grouped_div_4));
  printf ("(1+sqrt(2)/2)/2 => %s\n",
          cf2string (cf_half_of_1_plus_half_sqrt2));
  printf ("(1+1/sqrt(2))/2 => %s\n",
          cf2string (cf_half_of_1_plus_1_div_sqrt2));

  return 0;
}

/*------------------------------------------------------------------*/
Output:
$ cc univariate-continued-fraction-task.c -lgc && ./a.out
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2+sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(1+sqrt(2)/2)/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

Without a garbage collector

Translation of: ATS

This code is ported from ATS code that uses linear typing to ensure good malloc/free discipline. I try to retain the memory-use pattern of the ATS, to be safe against double-free errors, and hopefully to avoid memory leaks. (Absence of leaks is harder to be sure of.) When a generator is used as source of terms for another generator, the latter generator "consumes" the former (as if the types were linear), severely limiting reusability. Because of this, reusability of a generator is limited, and so I do not bother with memoization.

I use the [[maybe_unused]] notation of C23. You can simply remove the very few instances of it, if they bother you.

For no good reason, the program writes its output as code for input to groff.

/*------------------------------------------------------------------*/
/* For C23 without need of a garbage collector. */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>

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

/* We need consistent definitions of division and remainder. Let us
   set those here. For convenience (because C provides it), we will
   use truncation towards zero. */
#define DIV(a, b) ((a) / (b))
#define REM(a, b) ((a) % (b))

/* Choose a memory allocator. (Ideally one should check for NULL
   return values, but for this pedagogical example let us skip
   that.) */
#define MALLOC_INIT() do { } while (0) /* A no-op. */
#define MALLOC malloc
#define REALLOC realloc
#define FREE free

/*------------------------------------------------------------------*/
/* The basics. */

/* The integer type. */
typedef long long int integer;

/* A generator is a recursive type that forms a tree. */
typedef struct generator *generator_t;
typedef struct generator_list *generator_list_t;
struct generator_list
{
  generator_t car;
  generator_list_t cdr;
};
typedef void generator_func_t (integer *workspace,
                               generator_list_t sources,
                               bool *term_exists,
                               integer *term);
struct generator
{
  generator_func_t *run;     /* What does the work. */
  size_t worksize;           /* The size of the workspace. */
  integer *initial;          /* The initial value of the workspace. */
  integer *workspace;        /* The workspace itself. */
  generator_list_t sources;  /* The sources of input terms. */
};

/* Reinitializes a generator. (Needed because there is no
   memoization.) */
void generator_t_initialize (generator_t);

/* Frees a generator. */
void generator_t_free (generator_t);

/*------------------------------------------------------------------*/
/* A function to print the output of a generator in a form suitable
   for eqn/troff. */

void ftroff_generator_output (FILE *, generator_t, int max_terms);

/*------------------------------------------------------------------*/
/* Some functions to make generators. */

/* For a rational number. */
generator_t r2cf_make (integer n, integer d);

/* For the square root of 2. */
generator_t sqrt2_make (void);

/* For a homographic function. */
generator_t hfunc_make (integer a1, integer a, integer b1, integer b,
                        generator_t source);

/*------------------------------------------------------------------*/
/* Implementations. */

void
generator_t_initialize (generator_t gen)
{
  if (gen != NULL)
    {
      memcpy (gen->workspace, gen->initial,
              gen->worksize * sizeof (integer));
      for (generator_list_t p = gen->sources; p != NULL; p = p->cdr)
        generator_t_initialize (p->car);
    }
}

void
generator_t_free (generator_t gen)
{
  if (gen != NULL)
    {
      FREE (gen->initial);
      FREE (gen->workspace);

      generator_list_t p = gen->sources;
      while (p != NULL)
        {
          generator_list_t q = p->cdr;
          generator_t_free (p->car);
          FREE (p);
          p = q;
        }

      FREE (gen);
    }
}

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

void
ftroff_generator_output (FILE *outf, generator_t gen, int max_terms)
{
  assert (1 <= max_terms);

  generator_t_initialize (gen);

  int terms_count = 0;
  int sep = 0;
  bool done = false;
  while (!done)
    {
      if (terms_count == max_terms)
        {
          fprintf (outf, ", ~ ... ~ ]");
          done = true;
        }
      else
        {
          bool term_exists;
          integer term;
          gen->run (gen->workspace, gen->sources, &term_exists,
                    &term);
          if (term_exists)
            {
              switch (sep)
                {
                case 0:
                  fprintf (outf, "[ ^ ");
                  sep = 1;
                  break;
                case 1:
                  fprintf (outf, "; ^ ");
                  sep = 2;
                  break;
                default:
                  fprintf (outf, " , ");
                  break;
                }
              fprintf (outf, "%jd", (intmax_t) term);
              terms_count += 1;
            }
          else
            {
              fprintf (outf, "^ ] ");
              done = true;
            }
        }
    }

  generator_t_initialize (gen);
}

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

static void
r2cf_run (integer *workspace,
          [[maybe_unused]] generator_list_t sources,
          bool *term_exists, integer *term)
{
  integer d = workspace[1];
  *term_exists = (d != 0);
  if (*term_exists)
    {
      integer n = workspace[0];
      integer q = DIV (n, d);
      integer r = REM (n, d);
      workspace[0] = d;
      workspace[1] = r;
      *term = q;
    }
}

generator_t
r2cf_make (integer n, integer d)
{
  generator_t gen = MALLOC (sizeof (*gen));
  gen->run = r2cf_run;
  gen->worksize = 2;
  gen->initial = MALLOC (gen->worksize * sizeof (integer));
  gen->workspace = MALLOC (gen->worksize * sizeof (integer));
  gen->initial[0] = n;
  gen->initial[1] = d;
  memcpy (gen->workspace, gen->initial,
          gen->worksize * sizeof (integer));
  gen->sources = NULL;
  return gen;
}

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

static void
sqrt2_run (integer *workspace,
           [[maybe_unused]] generator_list_t sources,
           bool *term_exists, integer *term)
{
  *term_exists = true;
  *term = workspace[0];
  workspace[0] = 2;
}

generator_t
sqrt2_make (void)
{
  generator_t gen = MALLOC (sizeof (*gen));
  gen->run = sqrt2_run;
  gen->worksize = 1;
  gen->initial = MALLOC (gen->worksize * sizeof (integer));
  gen->workspace = MALLOC (gen->worksize * sizeof (integer));
  gen->initial[0] = 1;
  memcpy (gen->workspace, gen->initial,
          gen->worksize * sizeof (integer));
  gen->sources = NULL;
  return gen;
}

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

static void
hfunc_take_term (integer *workspace, generator_list_t sources)
{
  generator_t src = sources->car;
  bool term_exists1;
  integer term1;
  src->run (src->workspace, src->sources, &term_exists1, &term1);
  integer a1 = workspace[0];
  integer b1 = workspace[2];
  if (term_exists1)
    {
      integer a = workspace[1];
      integer b = workspace[3];
      workspace[0] = a + (a1 * term1);
      workspace[1] = a1;
      workspace[2] = b + (b1 * term1);
      workspace[3] = b1;
    }
  else
    {
      workspace[1] = a1;
      workspace[3] = b1;
    }
}

static void
hfunc_run (integer *workspace, generator_list_t sources,
           bool *term_exists, integer *term)
{
  bool done = false;
  while (!done)
    {
      integer b1 = workspace[2];
      integer b = workspace[3];
      if (b1 == 0 && b == 0)
        {
          *term_exists = false;
          done = true;
        }
      else
        {
          integer a1 = workspace[0];
          integer a = workspace[1];
          if (b1 != 0 && b != 0)
            {
              integer q1 = DIV (a1, b1);
              integer q = DIV (a, b);
              if (q1 == q)
                {
                  workspace[0] = b1;
                  workspace[1] = b;
                  workspace[2] = a1 - (b1 * q);
                  workspace[3] = a - (b * q);
                  *term_exists = true;
                  *term = q;
                  done = true;
                }
              else
                hfunc_take_term (workspace, sources);
            }
          else
            hfunc_take_term (workspace, sources);
        }
    }
}

generator_t
hfunc_make (integer a1, integer a, integer b1, integer b,
            generator_t source)
{
  generator_t gen = MALLOC (sizeof (*gen));
  gen->run = hfunc_run;
  gen->worksize = 4;
  gen->initial = MALLOC (gen->worksize * sizeof (integer));
  gen->workspace = MALLOC (gen->worksize * sizeof (integer));
  gen->initial[0] = a1;
  gen->initial[1] = a;
  gen->initial[2] = b1;
  gen->initial[3] = b;
  memcpy (gen->workspace, gen->initial,
          gen->worksize * sizeof (integer));
  gen->sources = MALLOC (sizeof (struct generator_list));
  gen->sources->car = source;
  gen->sources->cdr = NULL;
  return gen;
}

/*------------------------------------------------------------------*/
/* Components of the demonstration. */

#define MAX_TERMS 20
#define GOES_TO " ~ -> ~ "
#define START_EQ ".EQ\n"
#define STOP_EQ "\n.EN\n"
#define NEW_LINE "\n"

void
ftroff_rational_number (FILE *outf, integer n, integer d)
{
  generator_t gen = r2cf_make (n, d);
  fprintf (outf, "%s %jd over %jd %s",
           START_EQ, (intmax_t) n, (intmax_t) d, GOES_TO);
  ftroff_generator_output (outf, gen, MAX_TERMS);
  fprintf (outf, "%s%s", STOP_EQ, NEW_LINE);
  generator_t_free (gen);
}

void
ftroff_sqrt2 (FILE *outf)
{
  generator_t gen = sqrt2_make ();
  fprintf (outf, "%s sqrt 2 %s", START_EQ, GOES_TO);
  ftroff_generator_output (outf, gen, MAX_TERMS);
  fprintf (outf, "%s%s", STOP_EQ, NEW_LINE);
  generator_t_free (gen);
}

void
ftroff_hfunc_of_rational_number (FILE *outf,
                                 const char *expr,
                                 integer a1, integer a,
                                 integer b1, integer b,
                                 integer n, integer d)
{
  generator_t gen = hfunc_make (a1, a, b1, b, r2cf_make (n, d));
  fprintf (outf, "%s %s %s", START_EQ, expr, GOES_TO);
  ftroff_generator_output (outf, gen, MAX_TERMS);
  fprintf (outf, "%s%s", STOP_EQ, NEW_LINE);
  generator_t_free (gen);
}

void
ftroff_hfunc_of_sqrt2 (FILE *outf, const char *expr,
                       integer a1, integer a, integer b1, integer b)
{
  generator_t gen = hfunc_make (a1, a, b1, b, sqrt2_make ());
  fprintf (outf, "%s %s %s", START_EQ, expr, GOES_TO);
  ftroff_generator_output (outf, gen, MAX_TERMS);
  fprintf (outf, "%s%s", STOP_EQ, NEW_LINE);
  generator_t_free (gen);
}

void
ftroff_complicated (FILE *outf)
{
  /* This function demonstrates a more complicated nesting of
     generators. */

  /* gen1 = 1/sqrt(2) */
  generator_t gen1 = hfunc_make (0, 1, 1, 0, sqrt2_make ());

  /* gen2 = 1 + gen1 */
  generator_t gen2 = hfunc_make (1, 1, 0, 1, gen1);

  /* gen = gen2 / 2 */
  generator_t gen = hfunc_make (1, 0, 0, 2, gen2);

  fprintf (outf, "%s {1 ~ + ~ { 1 over { sqrt 2 } }} over 2 %s",
           START_EQ, GOES_TO);
  ftroff_generator_output (outf, gen, MAX_TERMS);
  fprintf (outf, "%s%s", STOP_EQ, NEW_LINE);

  generator_t_free (gen);
}

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

int
main (void)
{
  MALLOC_INIT ();

  FILE *outf = stdout;

  /* Output is for "eqn -Tpdf | groff -Tpdf -P-p6i,5.5i -ms" */

  fprintf (outf, ".nr PO 0.25i\n");
  fprintf (outf, ".nr HM 0.25i\n");
  fprintf (outf, ".ps 14\n");

  ftroff_rational_number (outf, 13, 11);
  ftroff_rational_number (outf, 22, 7);
  ftroff_sqrt2 (outf);
  ftroff_hfunc_of_rational_number
    (outf, "{ 13 over 11 } ~ + ~ { 1 over 2 }",
     2, 1, 0, 2, 13, 11);
  ftroff_hfunc_of_rational_number
    (outf, "{ 22 over 7 } ~ + ~ { 1 over 2 }",
     2, 1, 0, 2, 22, 7);
  ftroff_hfunc_of_rational_number
    (outf, "{ ^ {\"\\s-3\" 22 over 7 \"\\s+3\"}} over 4",
     1, 0, 0, 4, 22, 7);
  ftroff_hfunc_of_sqrt2
    (outf, "{ sqrt 2 } over 2", 1, 0, 0, 2);
  ftroff_hfunc_of_sqrt2
    (outf, "1 over { sqrt 2 }", 0, 1, 1, 0);
  ftroff_hfunc_of_sqrt2
    (outf, "{ 2 ~ + ~ { sqrt 2 }} over 4", 1, 2, 0, 4);
  ftroff_complicated (outf);

  return 0;
}

/*------------------------------------------------------------------*/
Output:

Run something like the following:

$ gcc -Wall -Wextra -g -std=gnu2x univariate-continued-fraction-task-no-gc.c
$ ./a.out | eqn -Tpdf | groff -Tpdf -P-p6i,5.5i -ms > foo.pdf
$ okular foo.pdf

In place of okular, you can use your favorite PDF viewer. To make the PNG, I used ImageMagick and ran:

$ convert -flatten -quality 0 foo.pdf foo.png

The output of the C program that uses no garbage collector, as rendered by groff.

C++

Uses ContinuedFraction and r2cf from Continued_fraction/Arithmetic/Construct_from_rational_number#C++

/* Interface for all matrixNG classes
   Nigel Galloway, February 10th., 2013.
*/
class matrixNG {
  private:
  virtual void consumeTerm(){}
  virtual void consumeTerm(int n){}
  virtual const bool needTerm(){}
  protected: int cfn = 0, thisTerm;
             bool haveTerm = false;
  friend class NG;
};
/* Implement the babyNG matrix
   Nigel Galloway, February 10th., 2013.
*/
class NG_4 : public matrixNG {
  private: int a1, a, b1, b, t;
  const bool needTerm() {
    if (b1==0 and b==0) return false;
    if (b1==0 or b==0) return true; else thisTerm = a/b;
    if (thisTerm==(int)(a1/b1)){
      t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm;
      haveTerm=true; return false;
    }
    return true;
  }
  void consumeTerm(){a=a1; b=b1;}
  void consumeTerm(int n){t=a; a=a1; a1=t+a1*n; t=b; b=b1; b1=t+b1*n;}
  public:
  NG_4(int a1, int a, int b1, int b): a1(a1), a(a), b1(b1), b(b){}
};
/* Implement a Continued Fraction which returns the result of an arithmetic operation on
   1 or more Continued Fractions (Currently 1 or 2).
   Nigel Galloway, February 10th., 2013.
*/
class NG : public ContinuedFraction {
  private:
   matrixNG* ng;
   ContinuedFraction* n[2];
  public:
  NG(NG_4* ng, ContinuedFraction* n1): ng(ng){n[0] = n1;}
  NG(NG_8* ng, ContinuedFraction* n1, ContinuedFraction* n2): ng(ng){n[0] = n1; n[1] = n2;}
  const int nextTerm() {ng->haveTerm = false; return ng->thisTerm;}
  const bool moreTerms(){
    while(ng->needTerm()) if(n[ng->cfn]->moreTerms()) ng->consumeTerm(n[ng->cfn]->nextTerm()); else ng->consumeTerm();
    return ng->haveTerm;
  }
};

Testing

[1;5,2] + 1/2

int main() {
  NG_4 a1(2,1,0,2);
  r2cf n1(13,11);
  for(NG n(&a1, &n1); n.moreTerms(); std::cout << n.nextTerm() << " ");
  std::cout << std::endl;
  return 0;
}
Output:
1 1 2 7

[3;7] * 7/22

int main() {
  NG_4 a2(7,0,0,22);
  r2cf n2(22,7);
  for(NG n(&a2, &n2); n.moreTerms(); std::cout << n.nextTerm() << " ");
  std::cout << std::endl;
  return 0;
}
Output:
1

[3;7] + 1/22

int main() {
  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;7] divided by 4

int main() {
  NG_4 a4(1,0,0,4);
  r2cf n4(22,7);
  for(NG n(&a4, &n4); n.moreTerms(); std::cout << n.nextTerm() << " ");
  std::cout << std::endl;
  return 0;
}
Output:
0 1 3 1 2

First I generate as a continued fraction, then I obtain an approximate value using r2cf for comparison.

int main() {
  NG_4 a5(0,1,1,0);
  SQRT2 n5;
  int i = 0;
  for(NG n(&a5, &n5); n.moreTerms() and i++ < 20; std::cout << n.nextTerm() << " ");
  std::cout << "..." << std::endl;
  for(r2cf cf(10000000, 14142136); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
  std::cout << std::endl;
  return 0;
}
Output:
0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
0 1 2 2 2 2 2 2 2 2 2 6 1 2 4 1 1 2

First I generate as a continued fraction, then I obtain an approximate value using r2cf for comparison.

int main() {
  int i = 0;
  NG_4 a6(1,1,0,2);
  SQRT2 n6;
  for(NG n(&a6, &n6); n.moreTerms() and i++ < 20; std::cout << n.nextTerm() << " ");
  std::cout << "..." << std::endl;
  for(r2cf cf(24142136, 20000000); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
  std::cout << std::endl;
  return 0;
}
Output:
1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 ...
1 4 1 4 1 4 1 4 1 4 3 2 1 9 5

Common Lisp

Translation of: ATS
Translation of: Scheme

This implementation memoizes terms of continued fractions. It mostly follows the coding of the Scheme, which itself was translated from ATS. Tail recursions in the Scheme have been replaced with ordinary loops. (Common Lisp standards do not require optimization of tail calls.)

I have tested with various CL implementations, including SBCL, CLISP, ECL, Clozure CL.

(defstruct cf-record
  terminated-p                    ; Are these all the terms there are?
  m                               ; How many terms are memoized so far?
  memo                            ; Where terms are memoized.
  gen)                            ; A thunk that generates terms.

(deftype continued-fraction 'cf-record)

(defun make-continued-fraction (gen)
  (make-cf-record :terminated-p nil
                  :m 0
                  :memo (make-array '(8))
                  :gen gen))

(defun cf-get-more-terms (cf needed)
  (loop with term
        for i from (cf-record-m cf) upto needed
        do (setf term (funcall (cf-record-gen cf)))
        unless term
          do (setf (cf-record-terminated-p cf) t)
        end
        while term do (setf (aref (cf-record-memo cf) i) term)
        finally (setf (cf-record-m cf) i)))

(defun cf-update (cf needed)
  (cond ((cf-record-terminated-p cf) (progn))
        ((<= needed (cf-record-m cf)) (progn))
        ((<= needed (array-dimension (cf-record-memo cf) 0))
         (cf-get-more-terms cf needed))
        (t ;; Provide twice the room that might be needed.
         (let* ((n1 (+ needed needed))
                (memo1 (make-array (list n1))))
           (loop for i from 0 upto (1- (cf-record-m cf))
                 do (setf (aref memo1 i)
                          (aref (cf-record-memo cf) i)))
           (setf (cf-record-memo cf) memo1)
           (cf-get-more-terms cf needed)))))

(defun continued-fraction-ref (cf i)
  (cf-update cf (1+ i))
  (and (< i (cf-record-m cf))
       (aref (cf-record-memo cf) i)))

(defun continued-fraction-to-thunk (cf)
  ;; Make a generator from a continued fraction.
  (let ((i 0))
    (lambda ()
      (let ((term (continued-fraction-ref cf i)))
        (setf i (1+ i))
        term))))

(defun continued-fraction-to-string (cf &optional (max-terms 20))
  (loop with sep = 0
        with accum = "["
        with term
        for i from 0 upto (1- max-terms)
        do (setf term (continued-fraction-ref cf i))
        if term
          do (let ((sep-str (case sep
                              ((0) "")
                              ((1) ";")
                              ((2) ","))))
               (setf sep (min (1+ sep) 2))
               (setf accum (concatenate 'string accum sep-str
                                        (format nil "~A" term))))
        else
          do (setf accum (concatenate 'string accum "]"))
             (return accum)
        end
        finally (setf accum (concatenate 'string accum ",...]"))
                (return accum)))

(defun r2cf (x)
  ;; This algorithm works directly with exact rationals, rather
  ;; than numerator and denominator separately.
  (let ((ratnum (coerce x 'rational))
        (terminated-p nil))
    (make-continued-fraction
     (lambda ()
       (and (not terminated-p)
            (multiple-value-bind (q r) (floor ratnum)
              (if (zerop r)
                  (setf terminated-p t)
                  (setf ratnum (/ r)))
              q))))))

(defstruct homographic-function a1 a b1 b)

(defun apply-homographic-function (hfunc cf)
  (let* ((gen (continued-fraction-to-thunk cf))
         (state (copy-homographic-function hfunc)))
    (make-continued-fraction
     (lambda ()
       (loop
         do (let ((a1 (homographic-function-a1 state))
                  (a (homographic-function-a state))
                  (b1 (homographic-function-b1 state))
                  (b (homographic-function-b state)))
              (cond ((and (zerop b1) (zerop b)) (return nil))
                    ((and (not (zerop b1)) (not (zerop b)))
                     (let ((q1 (nth-value 0 (floor a1 b1)))
                           (q (nth-value 0 (floor a b))))
                       (when (= q1 q)
                         (setf state (make-homographic-function
                                      :a1 b1
                                      :a b
                                      :b1 (- a1 (* b1 q))
                                      :b (- a (* b q))))
                         (return q)))))
              (let ((term (funcall gen)))
                (if term
                    (setf state
                          (make-homographic-function
                           :a1 (+ a (* a1 term))
                           :a a1
                           :b1 (+ b (* b1 term))
                           :b b1))
                    (progn
                      (setf (homographic-function-a state) a1)
                      (setf (homographic-function-b state)
                            b1))))))))))

(defun make-hf (a1 a b1 b)
  (make-homographic-function :a1 a1 :a a :b1 b1 :b b))

(defun apply-hf (hfunc cf)
  (apply-homographic-function hfunc cf))

(defun cf2string (cf)
  (continued-fraction-to-string cf))

(defvar cf+1/2 (make-hf 2 1 0 2))
(defvar cf/2 (make-hf 1 0 0 2))
(defvar cf/4 (make-hf 1 0 0 4))
(defvar 1/cf (make-hf 0 1 1 0))
(defvar 2+cf./4 (make-hf 1 2 0 4))
(defvar 1+cf./2 (make-hf 1 1 0 2))

(defvar cf_13/11 (r2cf 13/11))
(defvar cf_22/7 (r2cf 22/7))
(defvar cf_sqrt2
  (let ((next-term 1))
    (make-continued-fraction
     (lambda ()
       (let ((term next-term))
         (setf next-term 2)
         term)))))

(format t "13/11 => ~A~%" (cf2string cf_13/11))
(format t "22/7 => ~A~%" (cf2string cf_22/7))
(format t "sqrt(2) => ~A~%" (cf2string cf_sqrt2))
(format t  "13/11 + 1/2 => ~A~%"
        (cf2string (apply-hf cf+1/2 cf_13/11)))
(format t  "22/7 + 1/2 => ~A~%"
        (cf2string (apply-hf cf+1/2 cf_22/7)))
(format t  "(22/7)/4 => ~A~%"
        (cf2string (apply-hf cf/4 cf_22/7)))
(format t  "sqrt(2)/2 => ~A~%"
        (cf2string (apply-hf cf/2 cf_sqrt2)))
(format t  "1/sqrt(2) => ~A~%"
        (cf2string (apply-hf 1/cf cf_sqrt2)))
(format t  "(2 + sqrt(2))/4 => ~A~%"
        (cf2string (apply-hf 2+cf./4 cf_sqrt2)))
(format t  "(1 + 1/sqrt(2))/2 => ~A~%"
        (cf2string (apply-hf 1+cf./2 (apply-hf 1/cf cf_sqrt2))))
(format t  "sqrt(2)/4 + 1/2 => ~A~%"
        (cf2string (apply-hf cf+1/2 (apply-hf cf/4 cf_sqrt2))))
(format t  "(sqrt(2)/2)/2 + 1/2 => ~A~%"
        (cf2string (apply-hf cf+1/2
                            (apply-hf cf/2
                                     (apply-hf cf/2 cf_sqrt2)))))
(format t  "(1/sqrt(2))/2 + 1/2 => ~A~%"
        (cf2string (apply-hf cf+1/2
                            (apply-hf cf/2
                                     (apply-hf 1/cf cf_sqrt2)))))
Output:

SBCL might be the most likely CL implementation to be installed:

$ sbcl --script univariate-continued-fraction-task.lisp
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]
sqrt(2)/4 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(sqrt(2)/2)/2 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(1/sqrt(2))/2 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]

D

Translation of: ATS
Works with: Digital Mars D version 2.099.1
Works with: GCC version 12.2.1

This implementation memoizes terms of a continued fraction. It leaks memory, but D provides a garbage collector.

import std.conv;
import std.stdio;

alias index_t = uint; // The type for indexing terms of a continued
                      // fraction.

alias integer = long;   // The type for terms of a continued fraction.

class cf_t    //  A continued fraction, with memoization of its terms.
{
  protected bool terminated;   // Are there more terms to be memoized?
  protected index_t m;         // How many terms are memoized so far?
  private integer[] memo;      // Memoized terms.

  public index_t maxTerms = 20; // Maximum number of terms in the
                                // string representation.

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

  protected void generate (ref bool termExists, ref integer term)
  {
    // Return terms for zero. To get different terms, override this
    // method. (I am used to using a closure or similar for the
    // generator, and not having to derive a new continued fraction
    // type, to have a new kind of generator. However, I am trying to
    // do what is more natural within the programming language.)
    termExists = (m == 0);
    term = 0;
  }

  public void getAt (index_t i, ref bool termExists, ref integer term)
  {
    void memoizeMoreTerms (index_t needed)
    {
      while (m != needed && !terminated)
        {
          bool termExists;
          integer term;
          generate (termExists, term);
          if (termExists)
            {
              memo[m] = term;
              m += 1;
            }
          else
            terminated = true;
        }
    }

    void update (index_t needed)
    {
      // If necessary, memoize more terms, perhaps increasing the
      // space in which to store them.
      if (!terminated && m < needed)
        {
          if (memo.length < needed)
            {
              // Increase the space to twice what might be needed
              // right now.
              memo.length = 2 * needed;
            }
          memoizeMoreTerms (needed);
        }
    }

    update (i + 1);
    termExists = (i < m);
    if (termExists)
      term = memo[i];
  }

  public override string toString ()
  {
    string s = "[";
    int sep = 0;
    index_t i = 0;
    bool done = false;
    while (!done)
      {
        if (i == maxTerms)
          {
            s ~= ",...]";
            done = true;
          }
        else
          {
            bool termExists;
            integer term;
            getAt (i, termExists, term);
            if (termExists)
              {
                final switch (sep)
                  {
                  case 0 :
                    sep = 1;
                    break;
                  case 1 :
                    s ~= ";";
                    sep = 2;
                    break;
                  case 2 :
                    s ~= ",";
                    break;
                  }
                s ~= to!string (term);
                i += 1;
              }
            else
              {
                s ~= "]";
                done = true;
              }
          }
      }
    return s;
  }
}

class cfSqrt2_t : cf_t          // A continued fraction for sqrt(2).
{
  override final void generate (ref bool termExists, ref integer term)
  {
    termExists = true;
    term = (m == 0 ? 1 : 2);
  }
}

class cfRational : cf_t // A continued fraction for a rational number.
{
  private integer n;            // Numerator.
  private integer d;            // Denominator.

  this (integer numerator, integer denominator)
  {
    assert (denominator != 0);
    n = numerator;
    d = denominator;
  }

  override void generate (ref bool termExists, ref integer term)
  {
    termExists = (d != 0);
    if (termExists)
      {
        auto q = n / d;
        auto r = n % d;
        n = d;
        d = r;
        term = q;
      }
  }
}

class hfunc_t                   // A homographic function.
{
  public integer a1;
  public integer a;
  public integer b1;
  public integer b;

  this (integer a1, integer a, integer b1, integer b)
  {
    this.a1 = a1;
    this.a = a;
    this.b1 = b1;
    this.b = b;
  }
}

class cfHfunc_t : cf_t // A continued fraction that is a homographic
                       // function of some other continued fraction.
{
  private integer a1;
  private integer a;
  private integer b1;
  private integer b;
  private cf_t gen;
  private index_t index;

  this (hfunc_t hfunc, cf_t gen)
  {
    a1 = hfunc.a1;
    a = hfunc.a;
    b1 = hfunc.b1;
    b = hfunc.b;
    this.gen = gen;
    index = 0;
  }

  override void generate (ref bool termExists, ref integer term)
  {
    bool done = false;
    while (!done)
      {
        if (b1 == 0 && b == 0)
          {
            termExists = false;
            done = true;
          }
        else if (b1 != 0 && b != 0)
          {
            auto q1 = a1 / b1;
            auto q = a / b;
            if (q1 == q)
              {
                const a1_ = a1;
                const a_ = a;
                const b1_ = b1;
                const b_ = b;
                a1 = b1_;
                a = b_;
                b1 = a1_ - (b1_ * q);
                b = a_ - (b_ * q);
                termExists = true;
                term = q;
                done = true;
              }
          }

        if (!done)
          {
            gen.getAt (index, termExists, term);
            index += 1;
            if (termExists)
              {
                const a1_ = a1;
                const a_ = a;
                const b1_ = b1;
                const b_ = b;
                a1 = a_ + (a1_ * term);
                a = a1_;
                b1 = b_ + (b1_ * term);
                b = b1_;
              }
            else
              {
                a = a1;
                b = b1;
              }
          }
      }
  }
}

int
main (char[][] args)
{
  auto hf_cf_add_1_2 = new hfunc_t (2, 1, 0, 2);
  auto hf_cf_div_2 = new hfunc_t (1, 0, 0, 2);
  auto hf_cf_div_4 = new hfunc_t (1, 0, 0, 4);
  auto hf_1_div_cf = new hfunc_t (0, 1, 1, 0);

  auto cf_13_11 = new cfRational (13, 11);
  auto cf_22_7 = new cfRational (22, 7);
  auto cf_sqrt2 = new cfSqrt2_t ();

  auto cf_13_11_add_1_2 = new cfHfunc_t (hf_cf_add_1_2, cf_13_11);
  auto cf_22_7_add_1_2 = new cfHfunc_t (hf_cf_add_1_2, cf_22_7);
  auto cf_22_7_div_4 = new cfHfunc_t (hf_cf_div_4, cf_22_7);
  auto cf_sqrt2_div_2 = new cfHfunc_t (hf_cf_div_2, cf_sqrt2);
  auto cf_1_div_sqrt2 = new cfHfunc_t (hf_1_div_cf, cf_sqrt2);
  auto cf_2_add_sqrt2__div_4 =
    new cfHfunc_t (new hfunc_t (1, 2, 0, 4), cf_sqrt2);
  auto cf_1_add_1_div_sqrt2__div_2 =
    new cfHfunc_t (new hfunc_t (1, 1, 0, 2), cf_1_div_sqrt2);
  auto cf_sqrt2_div_4_add_1_2 =
    new cfHfunc_t (hf_cf_add_1_2,
                   new cfHfunc_t (hf_cf_div_4, cf_sqrt2));

  void show (string expr, cf_t cf)
  {
    writeln (expr, cf.toString());
  }

  show ("13/11 => ", cf_13_11);
  show ("22/7 => ", cf_22_7);
  show ("sqrt(2) => ", cf_sqrt2);
  show ("13/11 + 1/2 => ", cf_13_11_add_1_2);
  show ("22/7 + 1/2 => ", cf_22_7_add_1_2);
  show ("(22/7)/4 => ", cf_22_7_div_4);
  show ("sqrt(2)/2 => ", cf_sqrt2_div_2);
  show ("1/sqrt(2) => ", cf_1_div_sqrt2);
  show ("(2 + sqrt(2))/4 => ", cf_2_add_sqrt2__div_4);
  show ("(1 + 1/sqrt(2))/2 => ", cf_1_add_1_div_sqrt2__div_2);
  show ("sqrt(2)/4 + 1/2 => ", cf_sqrt2_div_4_add_1_2);
  show ("(sqrt(2)/2)/2 + 1/2 => ",
           new cfHfunc_t (hf_cf_add_1_2,
                          new cfHfunc_t (hf_cf_div_2,
                                         cf_sqrt2_div_2)));

  // Demonstrate a deeper nesting of anonymous cf_t.
  show ("(1/sqrt(2))/2 + 1/2 => ",
        new cfHfunc_t (hf_cf_add_1_2,
                       new cfHfunc_t (hf_cf_div_2,
                                      new cfHfunc_t (hf_1_div_cf,
                                                     cf_sqrt2))));

  return 0;
}
Output:
$ dmd -g univariate_continued_fraction_task.d && ./univariate_continued_fraction_task
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]
sqrt(2)/4 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(sqrt(2)/2)/2 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(1/sqrt(2))/2 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]

Fortran

Translation of: ATS
Translation of: C

Unlike the ATS and C implementations upon which this Fortran code is based, there is no garbage collector. The memory management is tricky, and if you find bugs in it, please feel free to fix them.

(Fortran standards allow garbage collection, but the NAG compiler is the only Fortran compiler I know of that offers garbage collection as an option. I am using GNU Fortran.)

I have been liberal in the use of recursive declarations and block constructs. In this program they can only help, not hurt.

!---------------------------------------------------------------------

module continued_fractions
  !
  ! Continued fractions with memoization.
  !

  implicit none
  private

  public :: cf_generator_proc_t
  public :: cf_generator_t
  public :: cf_t

  public :: cf_generator_make
  public :: cf_make
  public :: cf_generator_make_from_cf

  public :: cf_get_at

  public :: cf2string_max_terms
  public :: cf2string_default_max_terms
  public :: cf2string
  public :: default_max_terms

  integer :: default_max_terms = 20

  interface
     subroutine cf_generator_proc_t (env, term_exists, term)
       class(*), intent(inout) :: env
       logical, intent(out) :: term_exists
       integer, intent(out) :: term
     end subroutine cf_generator_proc_t
  end interface

  type :: cf_generator_t
     procedure(cf_generator_proc_t), pointer, nopass :: proc
     class(*), pointer :: env
     integer :: refcount = 0
   contains
     final :: cf_generator_t_finalize
     procedure :: cf_generator_t_refcount_incr
     procedure :: cf_generator_t_refcount_decr
  end type cf_generator_t

  type :: cf_memo_t
     integer, pointer :: storage(:)
     integer :: refcount = 0
   contains
     final :: cf_memo_t_finalize
     procedure :: cf_memo_t_refcount_incr
     procedure :: cf_memo_t_refcount_decr
  end type cf_memo_t

  type :: cf_t
     logical :: terminated
     integer :: m
     integer :: n
     class(cf_memo_t), pointer :: memo
     class(cf_generator_t), pointer :: gen
   contains
     final :: cf_t_finalize
  end type cf_t

  interface cf2string
     !
     ! Overload the name "cf2string".
     !
     module procedure cf2string_max_terms
     module procedure cf2string_default_max_terms
  end interface

  type :: cf_generator_from_cf_env_t
     class(cf_t), pointer :: cf
     integer :: i
  end type cf_generator_from_cf_env_t

contains

  recursive subroutine cf_generator_make (gen, proc, env)
    type(cf_generator_t), intent(out), pointer :: gen
    interface
       subroutine proc (env, term_exists, term)
         class(*), intent(inout) :: env
         logical, intent(out) :: term_exists
         integer, intent(out) :: term
       end subroutine proc
    end interface
    class(*), pointer, intent(inout) :: env

    allocate (gen)
    gen%proc => proc
    gen%env => env
  end subroutine cf_generator_make

  subroutine cf_generator_t_refcount_incr (gen)
    class(cf_generator_t), intent(inout) :: gen
    gen%refcount = gen%refcount + 1
  end subroutine cf_generator_t_refcount_incr

  subroutine cf_generator_t_refcount_decr (gen)
    class(cf_generator_t), intent(inout) :: gen
    gen%refcount = gen%refcount - 1
  end subroutine cf_generator_t_refcount_decr

  recursive subroutine cf_generator_t_finalize (gen)
    type(cf_generator_t), intent(inout) :: gen
    deallocate (gen%env)
  end subroutine cf_generator_t_finalize

  subroutine cf_memo_t_refcount_incr (memo)
    class(cf_memo_t), intent(inout) :: memo
    memo%refcount = memo%refcount + 1
  end subroutine cf_memo_t_refcount_incr

  subroutine cf_memo_t_refcount_decr (memo)
    class(cf_memo_t), intent(inout) :: memo
    memo%refcount = memo%refcount - 1
  end subroutine cf_memo_t_refcount_decr

  recursive subroutine cf_memo_t_finalize (memo)
    type(cf_memo_t), intent(inout) :: memo
    deallocate (memo%storage)
  end subroutine cf_memo_t_finalize

  recursive subroutine cf_make (cf, gen)
    type(cf_t), pointer, intent(out) :: cf
    type(cf_generator_t), pointer, intent(inout) :: gen

    integer, parameter :: start_size = 8

    allocate (cf)
    allocate (cf%memo)
    allocate (cf%memo%storage(0 : start_size - 1))
    cf%terminated = .false.
    cf%m = 0
    cf%n = start_size
    cf%gen => gen

    call cf%memo%cf_memo_t_refcount_incr
    call cf%gen%cf_generator_t_refcount_incr
  end subroutine cf_make

  recursive subroutine cf_t_finalize (cf)
    type(cf_t), intent(inout) :: cf

    call cf%memo%cf_memo_t_refcount_decr
    if (cf%memo%refcount == 0) deallocate (cf%memo)

    call cf%gen%cf_generator_t_refcount_decr
    if (cf%gen%refcount == 0) deallocate (cf%gen)
  end subroutine cf_t_finalize

  recursive subroutine cf_generator_make_from_cf (gen, cf)
    !
    ! TAKE NOTE: deallocating gen DOES NOT deallocate cf. (Most likely
    ! you would not want it to do so.)
    !
    type(cf_generator_t), intent(out), pointer :: gen
    type(cf_t), pointer, intent(inout) :: cf

    type(cf_generator_from_cf_env_t), pointer :: env
    class(*), pointer :: p

    allocate (env)
    env%cf => cf
    env%i = 0

    p => env
    call cf_generator_make (gen, cf_generator_from_cf_proc, p)
  end subroutine cf_generator_make_from_cf

  recursive subroutine cf_generator_from_cf_proc (env, term_exists, term)
    class(*), intent(inout) :: env
    logical, intent(out) :: term_exists
    integer, intent(out) :: term

    select type (env)
    class is (cf_generator_from_cf_env_t)
       call cf_get_at (env%cf, env%i, term_exists, term)
       env%i = env%i + 1
    end select
  end subroutine cf_generator_from_cf_proc

  recursive subroutine cf_get_more_terms (cf, needed)
    class(cf_t), intent(inout) :: cf
    integer, intent(in) :: needed

    integer :: term_count
    logical :: done

    logical :: term_exists
    integer :: term

    term_count = cf%m
    done = .false.
    do while (.not. done)
       if (term_count == needed) then
          cf%m = needed
          done = .true.
       else
          call cf%gen%proc (cf%gen%env, term_exists, term)
          if (term_exists) then
             cf%memo%storage(term_count) = term
             term_count = term_count + 1
          else
             cf%terminated = .true.
             cf%m = term_count
             done = .true.
          end if
       end if
    end do
  end subroutine cf_get_more_terms

  recursive subroutine cf_update (cf, needed)
    class(cf_t), intent(inout) :: cf
    integer, intent(in) :: needed

    integer, pointer :: storage1(:)

    if (cf%terminated .or. needed <= cf%m) then
       continue
    else if (needed <= cf%n) then
       call cf_get_more_terms (cf, needed)
    else
       ! Provide twice the needed storage.
       cf%n = 2 * needed
       allocate (storage1(0:cf%n - 1))
       storage1(0:cf%m - 1) = cf%memo%storage(0:cf%m - 1)
       deallocate (cf%memo%storage)
       cf%memo%storage => storage1
       call cf_get_more_terms (cf, needed)
    end if
  end subroutine cf_update

  recursive subroutine cf_get_at (cf, i, term_exists, term)
    class(cf_t), intent(inout) :: cf
    integer, intent(in) :: i
    logical, intent(out) :: term_exists
    integer, intent(out) :: term

    call cf_update (cf, i + 1)
    term_exists = (i < cf%m)
    if (term_exists) term = cf%memo%storage(i)
  end subroutine cf_get_at

  recursive function cf2string_max_terms (cf, max_terms) result (s)
    class(cf_t), intent(inout) :: cf
    integer, intent(in) :: max_terms
    character(len = :), allocatable :: s

    integer :: sep
    integer :: i, j
    logical :: done

    logical :: term_exists
    integer :: term

    character(len = 100) :: buf

    s = "["
    sep = 0
    i = 0
    done = .false.
    do while (.not. done)
       if (i == max_terms) then
          s = s // ",...]"
          done = .true.
       else
          call cf_get_at (cf, i, term_exists, term)
          if (term_exists) then
             select case (sep)
             case(0)
                sep = 1
             case(1)
                s = s // ";"
                sep = 2
             case(2)
                s = s // ","
             end select

             write (buf, '(I100)') term
             j = 1
             do while (buf(j:j) == ' ')
                j = j + 1
             end do
             s = s // buf(j:100)

             i = i + 1
          else
             s = s // "]"
             done = .true.
          end if
       end if
    end do
  end function cf2string_max_terms

  recursive function cf2string_default_max_terms (cf) result (s)
    class(cf_t), intent(inout) :: cf
    character(len = :), allocatable :: s
    s = cf2string_max_terms (cf, default_max_terms)
  end function cf2string_default_max_terms

end module continued_fractions

!---------------------------------------------------------------------

module continued_fractions_r2cf
  !
  ! Rational numbers.
  !

  use, non_intrinsic :: continued_fractions

  implicit none

  public :: r2cf_generator_make
  public :: r2cf_make

  type :: r2cf_generator_env_t
     integer :: n, d
  end type r2cf_generator_env_t

contains

  recursive subroutine r2cf_generator_make (gen, n, d)
    type(cf_generator_t), pointer, intent(out) :: gen
    integer, intent(in) :: n, d

    type(r2cf_generator_env_t), pointer :: env
    class(*), pointer :: p

    allocate (env)
    env%n = n
    env%d = d

    p => env
    call cf_generator_make (gen, r2cf_generator_proc, p)
  end subroutine r2cf_generator_make

  recursive subroutine r2cf_generator_proc (env, term_exists, term)
    class(*), intent(inout) :: env
    logical, intent(out) :: term_exists
    integer, intent(out) :: term

    integer :: q, r

    select type (env)
    class is (r2cf_generator_env_t)
       term_exists = (env%d /= 0)
       if (term_exists) then

          ! The direction in which to round the quotient is
          ! arbitrary. We will round it towards negative infinity.
          r = modulo (env%n, env%d)
          q = (env%n - r) / env%d

          env%n = env%d
          env%d = r

          term = q
       end if
    end select
  end subroutine r2cf_generator_proc

  recursive subroutine r2cf_make (cf, n, d)
    type(cf_t), pointer, intent(out) :: cf
    integer, intent(in) :: n, d

    type(cf_generator_t), pointer :: gen

    allocate (gen)
    call r2cf_generator_make (gen, n, d)
    call cf_make (cf, gen)
  end subroutine r2cf_make

end module continued_fractions_r2cf

!---------------------------------------------------------------------

module continued_fractions_sqrt2
  !
  ! The square root of two.
  !

  use, non_intrinsic :: continued_fractions

  implicit none

  public :: sqrt2_generator_make
  public :: sqrt2_make

  type :: sqrt2_generator_env_t
     integer :: term
  end type sqrt2_generator_env_t

contains

  recursive subroutine sqrt2_generator_make (gen)
    type(cf_generator_t), pointer, intent(out) :: gen

    type(sqrt2_generator_env_t), pointer :: env
    class(*), pointer :: p

    allocate (env)
    env%term = 1

    p => env
    call cf_generator_make (gen, sqrt2_generator_proc, p)
  end subroutine sqrt2_generator_make

  recursive subroutine sqrt2_generator_proc (env, term_exists, term)
    class(*), intent(inout) :: env
    logical, intent(out) :: term_exists
    integer, intent(out) :: term

    select type (env)
    class is (sqrt2_generator_env_t)
       term_exists = .true.
       term = env%term
       env%term = 2
    end select
  end subroutine sqrt2_generator_proc

  recursive subroutine sqrt2_make (cf)
    type(cf_t), pointer, intent(out) :: cf

    type(cf_generator_t), pointer :: gen

    allocate (gen)
    call sqrt2_generator_make (gen)
    call cf_make (cf, gen)
  end subroutine sqrt2_make

end module continued_fractions_sqrt2

!---------------------------------------------------------------------

module continued_fractions_hfunc
  !
  ! Homographic functions of cf_t objects.
  !

  use, non_intrinsic :: continued_fractions

  implicit none

  public :: hfunc_make

  type :: hfunc_generator_env_t
     integer :: a1, a, b1, b
     class(cf_generator_t), allocatable :: source_gen
  end type hfunc_generator_env_t

contains

  recursive subroutine hfunc_generator_make (gen, a1, a, b1, b, source_gen)
    type(cf_generator_t), pointer, intent(out) :: gen
    integer, intent(in) :: a1, a, b1, b
    type(cf_generator_t), pointer, intent(inout) :: source_gen

    type(hfunc_generator_env_t), pointer :: env
    class(*), pointer :: p

    allocate (env)
    env%a1 = a1
    env%a = a
    env%b1 = b1
    env%b = b
    env%source_gen = source_gen

    p => env
    call cf_generator_make (gen, hfunc_generator_proc, p)
  end subroutine hfunc_generator_make

  recursive subroutine hfunc_generator_proc (env, term_exists, term)
    class(*), intent(inout) :: env
    logical, intent(out) :: term_exists
    integer, intent(out) :: term

    integer :: q1, q
    logical :: done

    select type (env)
    class is (hfunc_generator_env_t)
       done = .false. 
       do while (.not. done)
          if (env%b1 == 0 .and. env%b == 0) then
             term_exists = .false.
             done = .true.
          else if (env%b1 /= 0 .and. env%b /= 0) then

             ! Because I feel like it, let us round quotients
             ! towards negative infinity.
             q1 = (env%a1 - modulo (env%a1, env%b1)) / env%b1
             q = (env%a - modulo (env%a, env%b)) / env%b

             if (q1 == q) then
                block
                  integer :: a1, a, b1, b
                  a1 = env%a1
                  a = env%a
                  b1 = env%b1
                  b = env%b
                  env%a1 = b1
                  env%a = b
                  env%b1 = a1 - (b1 * q)
                  env%b = a - (b * q)
                  term_exists = .true.
                  term = q
                  done = .true.
                end block
             end if
          end if

          if (.not. done) then
             call env%source_gen%proc (env%source_gen%env, term_exists, term)
             if (term_exists) then
                block
                  integer :: a1, a, b1, b
                  a1 = env%a1
                  a = env%a
                  b1 = env%b1
                  b = env%b
                  env%a1 = a + (a1 * term)
                  env%a = a1
                  env%b1 = b + (b1 * term)
                  env%b = b1
                end block
             else
                env%a = env%a1
                env%b = env%b1
             end if
          end if
       end do

    end select

  end subroutine hfunc_generator_proc

  recursive subroutine hfunc_make (cf, a1, a, b1, b, source_cf)
    type(cf_t), pointer, intent(out) :: cf
    integer, intent(in) :: a1, a, b1, b
    type(cf_t), pointer, intent(inout) :: source_cf

    type(cf_generator_t), pointer :: gen
    class(cf_generator_t), pointer :: source_gen

    call cf_generator_make_from_cf (source_gen, source_cf)
    call hfunc_generator_make (gen, a1, a, b1, b, source_gen)
    call cf_make (cf, gen)
  end subroutine hfunc_make

end module continued_fractions_hfunc

!---------------------------------------------------------------------

program univariate_continued_fraction_task

  use, non_intrinsic :: continued_fractions
  use, non_intrinsic :: continued_fractions_r2cf
  use, non_intrinsic :: continued_fractions_sqrt2
  use, non_intrinsic :: continued_fractions_hfunc

  implicit none

  type(cf_t), pointer :: cf_13_11
  type(cf_t), pointer :: cf_22_7
  type(cf_t), pointer :: cf_sqrt2

  type(cf_t), pointer :: cf_13_11_add_1_2
  type(cf_t), pointer :: cf_22_7_add_1_2
  type(cf_t), pointer :: cf_22_7_div_4
  type(cf_t), pointer :: cf_sqrt2_div_2
  type(cf_t), pointer :: cf_1_div_sqrt2
  type(cf_t), pointer :: cf_one_way
  type(cf_t), pointer :: cf_another_way

  type(cf_t), pointer :: cf_half_of_1_div_sqrt2
  type(cf_t), pointer :: cf_a_third_way

  call r2cf_make (cf_13_11, 13, 11)
  call r2cf_make (cf_22_7, 22, 7)
  call sqrt2_make (cf_sqrt2)

  call hfunc_make (cf_13_11_add_1_2, 2, 1, 0, 2, cf_13_11)
  call hfunc_make (cf_22_7_add_1_2, 2, 1, 0, 2, cf_22_7)
  call hfunc_make (cf_22_7_div_4, 1, 0, 0, 4, cf_22_7)
  call hfunc_make (cf_sqrt2_div_2, 1, 0, 0, 2, cf_sqrt2)
  call hfunc_make (cf_1_div_sqrt2, 0, 1, 1, 0, cf_sqrt2)
  call hfunc_make (cf_one_way, 1, 2, 0, 4, cf_sqrt2)
  call hfunc_make (cf_another_way, 1, 1, 0, 2, cf_1_div_sqrt2)

  call hfunc_make (cf_half_of_1_div_sqrt2, 1, 0, 0, 2, cf_1_div_sqrt2)
  call hfunc_make (cf_a_third_way, 2, 1, 0, 2, cf_half_of_1_div_sqrt2)

  write (*, '("13/11 => ", A)') cf2string (cf_13_11)
  write (*, '("22/7 => ", A)') cf2string (cf_22_7)
  write (*, '("sqrt(2) => ", A)') cf2string (cf_sqrt2)

  write (*, '("13/11 + 1/2 => ", A)') cf2string (cf_13_11_add_1_2)
  write (*, '("22/7 + 1/2 => ", A)') cf2string (cf_22_7_add_1_2)
  write (*, '("(22/7)/4 => ", A)') cf2string (cf_22_7_div_4)
  write (*, '("sqrt(2)/2 => ", A)') cf2string (cf_sqrt2_div_2)
  write (*, '("1/sqrt(2) => ", A)') cf2string (cf_1_div_sqrt2)
  write (*, '("(2 + sqrt(2))/4 => ", A)') cf2string (cf_one_way)
  write (*, '("(1 + 1/sqrt(2))/2 => ", A)') cf2string (cf_another_way)
  write (*, '("(1/sqrt(2))/2 + 1/2 => ", A)') cf2string (cf_a_third_way)

  deallocate (cf_13_11)
  deallocate (cf_22_7)
  deallocate (cf_sqrt2)
  deallocate (cf_13_11_add_1_2)
  deallocate (cf_22_7_add_1_2)
  deallocate (cf_22_7_div_4)
  deallocate (cf_sqrt2_div_2)
  deallocate (cf_1_div_sqrt2)
  deallocate (cf_one_way)
  deallocate (cf_another_way)
  deallocate (cf_half_of_1_div_sqrt2)
  deallocate (cf_a_third_way)

end program univariate_continued_fraction_task

!---------------------------------------------------------------------
Output:
$ gfortran -fbounds-check -Wall -Wextra -g -std=f2018 univariate_continued_fraction_task.f90 && ./a.out
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]
(1/sqrt(2))/2 + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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 ng4.go:

package cf

// A 2×2 matix:
//     [ a₁   a ]
//     [ b₁   b ]
//
// which when "applied" to a continued fraction representing x
// gives a new continued fraction z such that:
//
//         a₁ * x + a
//     z = ----------
//         b₁ * x + b
//
// Examples:
//      NG4{0, 1, 0, 4}.ApplyTo(x) gives 0*x + 1/4 -> 1/4 = [0; 4]
//      NG4{0, 1, 1, 0}.ApplyTo(x) gives 1/x
//      NG4{1, 1, 0, 2}.ApplyTo(x) gives (x+1)/2
//
// Note that several operations (e.g. addition and division)
// can be efficiently done with a single matrix application.
// However, each matrix application may require
// several calculations for each outputed term.
type NG4 struct {
	A1, A int64
	B1, B int64
}

func (ng NG4) needsIngest() bool {
	if ng.isDone() {
		panic("b₁==b==0")
	}
	return ng.B1 == 0 || ng.B == 0 || ng.A1/ng.B1 != ng.A/ng.B
}

func (ng NG4) isDone() bool {
	return ng.B1 == 0 && ng.B == 0
}

func (ng *NG4) ingest(t int64) {
	// [ a₁   a ] becomes [ a + a₁×t   a₁ ]
	// [ b₁   b ]         [ b + b₁×t   b₁ ]
	ng.A1, ng.A, ng.B1, ng.B =
		ng.A+ng.A1*t, ng.A1,
		ng.B+ng.B1*t, ng.B1
}

func (ng *NG4) ingestInfinite() {
	// [ a₁   a ] becomes [ a₁   a₁ ]
	// [ b₁   b ]         [ b₁   b₁ ]
	ng.A, ng.B = ng.A1, ng.B1
}

func (ng *NG4) egest(t int64) {
	// [ a₁   a ] becomes [      b₁         b   ]
	// [ b₁   b ]         [ a₁ - b₁×t   a - b×t ]
	ng.A1, ng.A, ng.B1, ng.B =
		ng.B1, ng.B,
		ng.A1-ng.B1*t, ng.A-ng.B*t
}

// ApplyTo "applies" the matrix `ng` to the continued fraction `cf`,
// returning the resulting continued fraction.
func (ng NG4) ApplyTo(cf ContinuedFraction) ContinuedFraction {
	return func() NextFn {
		next := cf()
		done := false
		return func() (int64, bool) {
			if done {
				return 0, false
			}
			for ng.needsIngest() {
				if t, ok := next(); ok {
					ng.ingest(t)
				} else {
					ng.ingestInfinite()
				}
			}
			t := ng.A1 / ng.B1
			ng.egest(t)
			done = ng.isDone()
			return t, true
		}
	}
}

File ng4_test.go:

package cf

import (
	"fmt"
	"reflect"
	"testing"
)

func Example_NG4() {
	cases := [...]struct {
		r  Rat
		ng NG4
	}{
		{Rat{13, 11}, NG4{2, 1, 0, 2}},
		{Rat{22, 7}, NG4{2, 1, 0, 2}},
		{Rat{22, 7}, NG4{1, 0, 0, 4}},
	}
	for _, tc := range cases {
		cf := tc.r.AsContinuedFraction()
		fmt.Printf("%5s = %-9s with %v gives %v\n", tc.r, cf.String(), tc.ng,
			tc.ng.ApplyTo(cf),
		)
	}

	invSqrt2 := NG4{0, 1, 1, 0}.ApplyTo(Sqrt2)
	fmt.Println("    1/√2 =", invSqrt2)
	result := NG4{1, 1, 0, 2}.ApplyTo(Sqrt2)
	fmt.Println("(1+√2)/2 =", result)

	// Output:
	// 13/11 = [1; 5, 2] with {2 1 0 2} gives [1; 1, 2, 7]
	//  22/7 = [3; 7]    with {2 1 0 2} gives [3; 1, 1, 1, 4]
	//  22/7 = [3; 7]    with {1 0 0 4} gives [0; 1, 3, 1, 2]
	//     1/√2 = [0; 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...]
	// (1+√2)/2 = [1; 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, ...]
}
Output:
13/11 = [1; 5, 2] with {2 1 0 2} gives [1; 1, 2, 7]
 22/7 = [3; 7]    with {2 1 0 2} gives [3; 1, 1, 1, 4]
 22/7 = [3; 7]    with {1 0 0 4} gives [0; 1, 3, 1, 2]
    1/√2 = [0; 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...]
(1+√2)/2 = [1; 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, ...]

Haskell

Works with: GHC version 9.0.2

One might note that a lazy list automatically memoizes terms, but not with O(1) access times.

The continued fraction generated here for sqrt(2) is actually the continued fraction for a close rational approximation to sqrt(2). I borrowed the definition along with that of real2cf. The approximation is probably not what you would want in a practical application, but I thought the implementation was cool, and I did not feel like being pedantic (until writing this commentary). :)

-- A continued fraction is represented as a lazy list of Int.

-- We borrow real2cf from
-- https://rosettacode.org/wiki/Continued_fraction/Arithmetic/Construct_from_rational_number#Haskell
-- though here some names in it are changed.

import Data.Ratio ((%))

real2cf frac =
  let (quotient, remainder) = properFraction frac
  in (quotient : (if remainder == 0
                  then []
                  else real2cf (1 / remainder)))

apply_hfunc (a1, a, b1, b) cf =
  recurs (a1, a, b1, b, cf)
  where recurs (a1, a, b1, b, cf) =
          if b1 == 0 && b == 0 then
            []
          else if b1 /= 0 && b /= 0 then
            let q1 = div a1 b1
                q = div a b
            in
              if q1 == q then
                q : recurs (b1, b, a1 - (b1 * q), a - (b * q), cf)
              else
                recurs (take_term (a1, a, b1, b, cf))
          else recurs (take_term (a1, a, b1, b, cf))
          where take_term (a1, a, b1, b, cf) =
                  case cf of
                    [] -> (a1, a1, b1, b1, cf)
                    (term : cf1) -> (a + (a1 * term), a1,
                                     b + (b1 * term), b1,
                                     cf1)

cf_13_11 = real2cf (13 % 11)
cf_22_7 = real2cf (22 % 7)
cf_sqrt2 = real2cf (sqrt 2)

cfToString cf =
  loop 0 0 "[" cf
  where loop i sep s lst =
          case lst of
            [] -> s ++ "]"
            (term : tail) ->
              if i == 20
              then s ++ ",...]"
              else
                do loop (i + 1) sep1 s1 tail
                     where sepStr = case sep of
                                      0 -> ""
                                      1 -> ";"
                                      _ -> ","
                           sep1 = min (sep + 1) 2
                           termStr = show term
                           s1 = s ++ sepStr ++ termStr

show_cf expr cf =
  do putStr expr;
     putStr " => ";
     putStrLn (cfToString cf)

main =
  do show_cf "13/11" cf_13_11;
     show_cf "22/7" cf_22_7;
     show_cf "sqrt(2)" cf_sqrt2;
     show_cf "13/11 + 1/2" (apply_hfunc (2, 1, 0, 2) cf_13_11);
     show_cf "22/7 + 1/2" (apply_hfunc (2, 1, 0, 2) cf_22_7);
     show_cf "(22/7)/4" (apply_hfunc (1, 0, 0, 4) cf_22_7);
     show_cf "1/sqrt(2)" (apply_hfunc (0, 1, 1, 0) cf_sqrt2);
     show_cf "(2 + sqrt(2))/4" (apply_hfunc (1, 2, 0, 4) cf_sqrt2);
     show_cf "(1 + 1/sqrt(2))/2" (apply_hfunc (2, 1, 0, 2) -- cf + 1/2
                                  (apply_hfunc (1, 0, 0, 2)  -- cf/2
                                   (apply_hfunc (0, 1, 1, 0) -- 1/cf
                                     cf_sqrt2)))
Output:
$ ghc univariate_continued_fraction_task.hs && ./univariate_continued_fraction_task
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

Icon

Works with: Icon version 9.5.22e
Translation of: ATS

This implementation memoizes terms of a continued fraction.

# An implementation in Icon, using co-expressions as generators.

$define YES 1
$define NO  &null

# terminated = are there no more terms to memoize?
# memo       = memoized terms.
# generate   = a co-expression to generate more terms.
record continued_fraction (terminated, memo, generate)

procedure main ()
  local cf_13_11, cf_22_7, cf_sqrt2, cf_1_div_sqrt2

  cf_13_11 := make_cf_rational (13, 11)
  cf_22_7 := make_cf_rational (22, 7)
  cf_sqrt2 := make_cf_sqrt2()
  cf_1_div_sqrt2 := make_cf_hfunc (0, 1, 1, 0, cf_sqrt2)

  show ("13/11", cf_13_11)
  show ("22/7", cf_22_7)
  show ("sqrt(2)", cf_sqrt2)
  show ("13/11 + 1/2", make_cf_hfunc (2, 1, 0, 2, cf_13_11))
  show ("22/7 + 1/2", make_cf_hfunc (2, 1, 0, 2, cf_22_7))
  show ("(22/7)/4", make_cf_hfunc (1, 0, 0, 4, cf_22_7))
  show ("1/sqrt(2)", cf_1_div_sqrt2)
  show ("(2 + sqrt(2))/4", make_cf_hfunc (1, 2, 0, 4, cf_sqrt2))
  show ("(1 + 1/sqrt(2))/2", make_cf_hfunc (1, 1, 0, 2,
                                            cf_1_div_sqrt2))
end

procedure show (expr, cf)
  write (expr, " => ", cf2string (cf))
end

procedure make_cf_sqrt2 ()
  return make_continued_fraction (create gen_sqrt2 ())
end

procedure make_cf_rational (n, d)
  return make_continued_fraction (create gen_rational (n, d))
end

procedure make_cf_hfunc (a1, a, b1, b, other_cf)
  return make_continued_fraction (create gen_hfunc (a1, a, b1, b,
                                                    other_cf))
end

procedure gen_sqrt2 ()
  suspend 1
  repeat suspend 2
end

procedure gen_rational (n, d)
  local q, r

  repeat {
    if d = 0 then fail
    q := n / d
    r := n % d
    n := d
    d := r
    suspend q
  }
end

procedure gen_hfunc (a1, a, b1, b, other_cf)
  local a1_tmp, a_tmp, b1_tmp, b_tmp
  local i, term, skip_getting_a_term
  local q1, q

  i := 0
  repeat {
    skip_getting_a_term := NO
    if b1 = b = 0 then {
      fail
    } else if b1 ~= 0 & b ~= 0 then {
      q1 := a1 / b1
      q := a / b
      if q1 = q then {
        a1_tmp := a1
        a_tmp := a
        b1_tmp := b1
        b_tmp := b
        a1 := b1_tmp
        a := b_tmp
        b1 := a1_tmp - (b1_tmp * q)
        b := a_tmp - (b_tmp * q)
        suspend q
        skip_getting_a_term := YES
      }
    }
    if /skip_getting_a_term then {
      if term := get_term (other_cf, i) then {
        i +:= 1
        a1_tmp := a1
        a_tmp := a
        b1_tmp := b1
        b_tmp := b
        a1 := a_tmp + (a1_tmp * term)
        a := a1_tmp
        b1 := b_tmp + (b1_tmp * term)
        b := b1_tmp
      } else {
        a := a1
        b := b1
      }
    }
  }
end

procedure make_continued_fraction (gen)
  return continued_fraction (NO, [], gen)
end

procedure get_term (cf, i)
  local j, term

  if *cf.memo <= i then {
    if \cf.terminated then {
      fail
    } else {
      every j := *cf.memo to i do {
        if term := @(cf.generate) then {
          put (cf.memo, term)
        } else {
          cf.terminated := YES
          fail
        }
      }
    }
  }
  return cf.memo[i + 1]
end

procedure cf2string (cf, max_terms)
  local s, sep, i, done, term

  /max_terms := 20

  s := "["
  sep := 0
  i := 0
  done := NO
  while /done do {
    if i = max_terms then {
      # We have reached the maximum of terms to print. Stick an
      # ellipsis in the notation.
      s ||:= ",...]"
      done := YES
    } else if term := get_term (cf, i) then {
      # Getting a term succeeded. Include the term.
      s ||:= sep_str (sep) || term
      sep := sep + 1
      if 2 < sep then sep := 2
      i +:= 1
    } else {
      # Getting a term failed. We are done.
      s ||:= "]"
      done := YES
    }
  }
  return s
end

procedure sep_str (sep)
  return (if sep = 0 then "" else if sep = 1 then ";" else ",")
end
Output:
icon  univariate-continued-fraction-task.icn
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

J

Note that the continued fraction representation used here differs from those implemented in the Continued_fraction task. In that task, we alternated a and b values. Here, we only work with a values -- b is implicitly always 1.

Implementation:

ng4cf=: 4 : 0
  cf=. 1000{.!._ y
  ng=. x
  r=.i. ndx=.0
  while. +./0~:{:ng do.
    if.=/<.%/ng do.
      r=.r, t=.{.<.%/ng
      ng=. t (|.@] - ]*0,[) ng 
    else.
      if. _=t=.ndx{cf do.
        ng=. ng+/ .*2 2$1 1 0 0
      else.
        ng=. ng+/ .*2 2$t,1 1 0
      end.
      if. (#cf)=ndx=. ndx+1 do. r return. end.
    end.
  end.
  r
)

Notes:

  • we arbitrarily stop processing continued fractions after 1000 elements. That's more than enough precision for most purposes.
  • we can convert a continued fraction to a rational number using (+%)/ though if we want the full represented precision we should instead use (+%)/@x: (which is slower).
  • we can convert a rational number to a continued fraction using 1 1 {."1@}. ({: , (0 , {:) #: {.)^:(*@{:)^:a: but also this expects a numerator,denominator pair so if you have only a single number use ,&1 to give it a denominator. This works equally well with floating point and arbitrary precision numbers.

Some arbitrary continued fractions and their floating point representations

   arbs=:(,1);(,3);?~&.>3+i.10
   ":@>arbs
1                        
3                        
1 2 0                    
0 2 3 1                  
1 0 3 2 4                
0 2 3 5 1 4              
2 5 0 1 6 3 4            
7 5 6 3 0 4 1 2          
7 0 1 2 6 3 8 4 5        
8 0 5 6 3 7 4 9 1 2      
0 9 8 1 3 10 2 5 6 7 4   
1 7 3 4 5 8 9 10 6 11 0 2
   (+%)/@>arbs
1 3 1 0.444444 4.44444 0.431925 2.16238 7.19368 8.46335 13.1583 0.109719 1.13682

Some NG based cf functions, verifying their behavior against our test set:

   plus1r2=: (2 1,:0 2)&ng4cf
   (plus1r2 each  -&((+%)/@>) ]) arbs 
0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

For every one of our arbitrary continued fractions, the 2 1,:0 2 NG matrix gives us a new continued fraction whose rational value is the original rational value + 1r2.

   times7r22=: (7 0,:0 22)&ng4cf 
   (times7r22 each %&((+%)/@>) ]) arbs 
0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182 0.318182
   (times7r22 each %&((+%)/@x:@>) ]) arbs 
7r22 7r22 7r22 7r22 7r22 7r22 7r22 7r22 7r22 7r22 7r22 7r22

For every one of our arbitrary continued fractions, the 7 0,:0 22 NG matrix gives us a new continued fraction whose rational value is 7r22 times the original rational value.

   times1r4=:(1 0,:0 4)&ng4cf
   (times1r4 each %&((+%)/@>) ]) arbs 
0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25

It seems like a diagonal matrix has the effect of multiplying the numerator by the upper left element and the denominator by the lower right element. And our first experiment suggests that an upper right element in NG with a 0 for the bottom left will add the top right divided by bottom right to our continued fraction.

   reciprocal=:(0 1,:1 0)&ng4cf
   (reciprocal each *&((+%)/@>) ]) arbs 
1 1 1 1 1 1 1 1 1 1 1 1

Looks like we can also divide by our continued fraction...

   plus1r2times1r2=: (1 1,:0 2)&ng4cf
   (plus1r2times1r2 each (= 0.5+0.5*])&((+%)/@>) ]) arbs 
1 1 1 1 1 1 1 1 1 1 1 1

We can add and multiply using a single "ng4" operation.

Task examples:

1r2 + 13r11

   (+%)/1 5 2
1.18182
   plus1r2 1 5 2
1 1 2 7
   (+%)/plus1r2 1 5 2
1.68182

7r22 * 22r7

   (+%)/3 7x
22r7
   times7r22 3 7x
1

1r2 + 22r7

   plus1r2 3 7x
3 1 1 1 4
   (+%)/plus1r2 3 7x
3.64286
   (+%)/x:plus1r2 3 7x
51r14

1r4 * 22r7

   times1r4 3 7x
0 1 3 1 2
   (+%)/x:times1r4 3 7x
11r14

   reciprocal 1,999$2
0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
   (+%)/1,999$2
1.41421
   (+%)/reciprocal 1,999$2
0.707107

   plus1r2times1r2 1,999$2
1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 ...
   (+%)/plus1r2times1r2 1,999$2
1.20711

   plus1r2times1r2 0 1,999$2
0 1 5 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 ...
   (+%)/plus1r2times1r2 0 1,999$2
0.853553

Java

import java.util.List;

public final class ContinuedFractionArithmeticG1 {

	public static void main(String[] aArgs) {
		List<CFData> cfData = List.of(
	        new CFData("[1; 5, 2] + 1 / 2    ", new int[] { 2, 1, 0, 2 }, (CFIterator) new R2cfIterator(13, 11) ),
	        new CFData("[3; 7] + 1 / 2       ", new int[] { 2, 1, 0, 2 }, (CFIterator) new R2cfIterator(22, 7) ),
	        new CFData("[3; 7] divided by 4  ", new int[] { 1, 0, 0, 4 }, (CFIterator) new R2cfIterator(22, 7) ),
	        new CFData("sqrt(2)              ", new int[] { 0, 1, 1, 0 }, (CFIterator) new ReciprocalRoot2() ),
	        new CFData("1 / sqrt(2)          ", new int[] { 0, 1, 1, 0 }, (CFIterator) new Root2() ),
	        new CFData("(1 + sqrt(2)) / 2    ", new int[] { 1, 1, 0, 2 }, (CFIterator) new Root2() ),
	        new CFData("(1 + 1 / sqrt(2)) / 2", new int[] { 1, 1, 0, 2 }, (CFIterator) new ReciprocalRoot2() ) );       
		
		for ( CFData data : cfData ) {
            System.out.print(data.text + " -> ");
			NG ng = new NG(data.arguments);
			CFIterator iterator = data.iterator;			
			int nextTerm = 0;
			
			for ( int i = 1; i <= 20 && iterator.hasNext(); i++ ) {			
				nextTerm = iterator.next();					
		        if ( ! ng.needsTerm() ) {
		        	System.out.print(ng.egress() + " ");
		        } 
		        ng.ingress(nextTerm); 
	        } 
			
			while ( ! ng.done() ) {
	            System.out.print(ng.egressDone() + " ");
	        }
	        System.out.println();
		}

	}
	
	private static class NG {
		
		public NG(int[] aArgs) {
			a1 = aArgs[0]; a = aArgs[1]; b1 = aArgs[2]; b = aArgs[3];
		}

	    public void ingress(int aN) {
	        int temp = a; a = a1; a1 = temp + a1 * aN; 
	            temp = b; b = b1; b1 = temp + b1 * aN;
	    }

	    public int egress() {
	        final int n = a / b;
	        int temp = a;  a  = b;  b  = temp - b * n;
	            temp = a1; a1 = b1; b1 = temp - b1 * n;
	        return n;
	    }

	    public boolean needsTerm() {
	    	return ( b == 0 || b1 == 0 ) || ( a * b1 != a1 * b );
	    }
	    
	    public int egressDone() {	        
            if ( needsTerm() ) {
                a = a1; 
                b = b1;
            }
            return egress();
	    }
	     
	    public boolean done() {
	    	return ( b == 0 || b1 == 0 );
	    }
	    		
	    private int a1, a, b1, b; 
	}

	private static abstract class CFIterator {
		
		public abstract boolean hasNext();
		public abstract int next();
		
	}	
	
	private static class R2cfIterator extends CFIterator {
		
		public R2cfIterator(int aNumerator, int aDenominator) {
			numerator = aNumerator; denominator = aDenominator;
		}
		
		public boolean hasNext() {
			return denominator != 0;
		}
		
		public int next() {
			int div = numerator / denominator;
            int rem = numerator % denominator;
            numerator = denominator;
            denominator = rem;
            return div;
		}
		
		private int numerator, denominator;
		
	}

	private static class Root2 extends CFIterator {
		
		public Root2() {
			firstReturn = true;
		}
		
		public boolean hasNext() {
			return true;
		}
		
		public int next() {
			if ( firstReturn ) {
				firstReturn = false;
				return 1;
			}
			return 2;
		}
		
		private boolean firstReturn;
		
	}
	
	private static class ReciprocalRoot2 extends CFIterator {
		
		public ReciprocalRoot2() {
			firstReturn = true;
			secondReturn = true;
		}
		
		public boolean hasNext() {
			return true;
		}
		
		public int next() {
			if ( firstReturn ) {
				firstReturn = false;
				return 0;
			}
			if ( secondReturn ) {
				secondReturn = false;
				return 1;
			}
			return 2;
		}
		
		private boolean firstReturn, secondReturn;
		
	}
	
	private static record CFData(String text, int[] arguments, CFIterator iterator) {} 
	
}
Output:
[1; 5, 2] + 1 / 2     -> 1 1 2 7 
[3; 7] + 1 / 2        -> 3 1 1 1 4 
[3; 7] divided by 4   -> 0 1 3 1 2 
sqrt(2)               -> 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
1 / sqrt(2)           -> 0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
(1 + sqrt(2)) / 2     -> 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 
(1 + 1 / sqrt(2)) / 2 -> 0 1 5 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 5 

Julia

Translation of: Ruby
function r2cf(n1::Integer, n2::Integer)
    ret = Int[]
    while n2 != 0
        n1, (t1, n2) = n2, divrem(n1, n2)
        push!(ret, t1)
    end
    ret
end
r2cf(r::Rational) = r2cf(numerator(r), denominator(r))

function r2cf(n1, n2, maxiter=20)
    ret = Int[]
    while n2 != 0 && maxiter > 0
        n1, (t1, n2) = n2, divrem(n1, n2)
        push!(ret, t1)
        maxiter -= 1
    end
    ret
end

mutable struct NG
    a1::Int
    a::Int
    b1::Int
    b::Int
end

function ingress(ng, n)
    ng.a, ng.a1= ng.a1, ng.a + ng.a1 * n
    ng.b, ng.b1 = ng.b1, ng.b + ng.b1 * n
end

needterm(ng) = ng.b == 0 || ng.b1 == 0 || !(ng.a // ng.b == ng.a1 // ng.b1)

function egress(ng)
    n = ng.a // ng.b
    ng.a, ng.b = ng.b, ng.a - ng.b * n
    ng.a1, ng.b1 = ng.b1, ng.a1 - ng.b1 * n
    r2cf(n)
end

egress_done(ng) = (if needterm(ng) ng.a, ng.b = ng.a1, ng.b1 end; egress(ng))

done(ng) = ng.b == 0 && ng.b1 == 0

function testng()
    data = [["[1;5,2] + 1/2",      [2,1,0,2], [13,11]],
        ["[3;7] + 1/2",        [2,1,0,2], [22, 7]],
        ["[3;7] divided by 4", [1,0,0,4], [22, 7]],
        ["[1;1] divided by sqrt(2)", [0,1,1,0], [1,sqrt(2)]]]

    for d in data
        str, ng, r = d[1], NG(d[2]...), d[3]
        print(rpad(str, 25), "->")
        for n in r2cf(r...)
            if !needterm(ng)
                print(" $(egress(ng))")
            end
            ingress(ng, n)
        end
        while true
            print(" $(egress_done(ng))")
            if done(ng)
                println()
                break
            end
        end
    end
end

testng()
Output:
[1;5,2] + 1/2            -> [1, 1, 2, 7]
[3;7] + 1/2              -> [3, 1, 1, 1, 4]
[3;7] divided by 4       -> [0, 1, 3, 1, 2]
[1;1] divided by sqrt(2) -> [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

Kotlin

This is based on the Python entry but has been expanded to deal with the '√2' calculations:

// version 1.1.3
// compile with -Xcoroutines=enable flag from command line
 
import kotlin.coroutines.experimental.*

typealias CFGenerator = (Pair<Int, Int>) -> Sequence<Int>

data class CFData( 
    val str: String, 
    val ng: IntArray,
    val r: Pair<Int,Int>,
    val gen: CFGenerator
)
 
fun r2cf(frac: Pair<Int, Int>) = 
    buildSequence {
        var num = frac.first
        var den = frac.second
        while (Math.abs(den) != 0) {
            val div = num / den
            val rem = num % den
            num = den
            den = rem
            yield(div)
        }
    }

fun d2cf(d: Double) = 
    buildSequence {
        var dd  = d
        while (true) {
            val div = Math.floor(dd)
            val rem = dd - div
            yield(div.toInt())
            if (rem == 0.0) break
            dd = 1.0 / rem
        }
    }

@Suppress("UNUSED_PARAMETER")
fun root2(dummy: Pair<Int, Int>) =
    buildSequence {
        yield(1)
        while (true) yield(2)
    }

@Suppress("UNUSED_PARAMETER")
fun recipRoot2(dummy: Pair<Int, Int>) =
    buildSequence {
       yield(0)
       yield(1)
       while (true) yield(2)
    }
 
class NG(var a1: Int, var a: Int, var b1: Int, var b: Int) {

    fun ingress(n: Int) {
        var t = a
        a = a1
        a1 = t + a1 * n
        t = b
        b = b1
        b1 = t + b1 * n
    }

    fun egress(): Int {
        val n = a / b
        var t = a
        a = b
        b = t - b * n
        t = a1
        a1 = b1
        b1 = t - b1 * n
        return n
    }

    val needTerm get() = (b == 0 || b1 == 0) || ((a / b) != (a1 / b1))
    
    val egressDone: Int
        get() {
            if (needTerm) {
                a = a1
                b = b1
            }
            return egress()
        }
        
    val done get() = b == 0 &&  b1 == 0
}

fun main(args: Array<String>) {
    val data = listOf(
        CFData("[1;5,2] + 1/2        ", intArrayOf(2, 1, 0, 2), 13 to 11, ::r2cf),
        CFData("[3;7] + 1/2          ", intArrayOf(2, 1, 0, 2), 22 to 7,  ::r2cf),
        CFData("[3;7] divided by 4   ", intArrayOf(1, 0, 0, 4), 22 to 7,  ::r2cf),
        CFData("sqrt(2)              ", intArrayOf(0, 1, 1, 0),  0 to 0,  ::recipRoot2),
        CFData("1 / sqrt(2)          ", intArrayOf(0, 1, 1, 0),  0 to 0,  ::root2),
        CFData("(1 + sqrt(2)) / 2    ", intArrayOf(1, 1, 0, 2),  0 to 0,  ::root2),
        CFData("(1 + 1 / sqrt(2)) / 2", intArrayOf(1, 1, 0, 2),  0 to 0,  ::recipRoot2)       
    )
    println("Produced by NG class:")
    for ((str, ng, r, gen) in data) {
        print("$str -> ")
        val (a1, a, b1, b) = ng
        val op = NG(a1, a, b1, b)        
        for (n in gen(r).take(20)) {
            if (!op.needTerm) print(" ${op.egress()} ") 
            op.ingress(n) 
        } 
        while (true) {
            print(" ${op.egressDone} ")
            if (op.done) break
        }
        println()
    }
    println("\nProduced by direct calculation:")
    val data2 = listOf(
        Pair("(1 + sqrt(2)) / 2    ", (1 + Math.sqrt(2.0)) / 2), 
        Pair("(1 + 1 / sqrt(2)) / 2", (1 + 1 / Math.sqrt(2.0)) / 2)
    )
    for ((str, d) in data2) {
        println("$str ->  ${d2cf(d).take(20).joinToString("  ")}")
    }
}
Output:
Produced by NG class:
[1;5,2] + 1/2         ->  1  1  2  7 
[3;7] + 1/2           ->  3  1  1  1  4 
[3;7] divided by 4    ->  0  1  3  1  2 
sqrt(2)               ->  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 
1 / sqrt(2)           ->  0  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 
(1 + sqrt(2)) / 2     ->  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4 
(1 + 1 / sqrt(2)) / 2 ->  0  1  5  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  5 

Produced by direct calculation:
(1 + sqrt(2)) / 2     ->  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  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

Mercury

Translation of: Haskell
%%%-------------------------------------------------------------------

:- module univariate_continued_fraction_task_lazy.

:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.

:- implementation.
:- import_module integer.       % Arbitrary-precision integers.
:- import_module lazy.          % Lazy evaluation.
:- import_module list.
:- import_module rational.      % Arbitrary-precision fractions.
:- import_module string.

%%%-------------------------------------------------------------------
%%%
%%% The following lazy list implementation is suggested in the Mercury
%%% Library Reference, although (for convenience) I have changed the
%%% names.
%%%

:- type lzlist(T)
   ---> lzlist(lazy(lzcell(T))).

:- type lzcell(T)
   ---> lzcons(T, lzlist(T))
   ;    lznil.

%%%-------------------------------------------------------------------
%%%
%%% Types of interest.
%%%

:- type cf == lzlist(integer).  % A continued fraction.
:- type hf == {integer, integer,
               integer, integer}. % A homographic function.
:- type ng4 == hf.                % A synonym for hf.


%%%-------------------------------------------------------------------
%%%
%%% Make a "continued fraction" that has no terms.
%%%

:- func cfnil = cf.
cfnil = lzlist(delay((func) = lznil)).

%%%-------------------------------------------------------------------
%%%
%%% Make a continued fraction that repeats the same term forever.
%%%

:- func repeat_forever(integer) = cf.

repeat_forever(N) = CF :-
  CF = lzlist(delay(Cons)),
  Cons = ((func) = lzcons(N, repeat_forever(N))).

%%%-------------------------------------------------------------------
%%%
%%% sqrt2 is a continued fraction for the square root of two.
%%%

:- func sqrt2 = cf.

sqrt2 = lzlist(delay((func) = lzcons(one, repeat_forever(two)))).

%%%-------------------------------------------------------------------
%%%
%%% r2cf takes a fraction, and returns a continued fraction as a lazy
%%% list of terms.
%%%

:- func r2cf(rational) = cf.
:- func r2cf(integer, integer) = cf.

r2cf(Ratnum) = CF :-
  r2cf(numer(Ratnum), denom(Ratnum)) = CF.

r2cf(Numerator, Denominator) = CF :-
  (if (Denominator = zero)
   then (CF = cfnil)
   else (CF = lzlist(delay(Cons)),
         ((func) = X :-
            (X = lzcons(Quotient, r2cf(Denominator, Remainder)),
             %% What follows is division with truncation towards zero.
             divide_with_rem(Numerator, Denominator,
                             Quotient, Remainder))) = Cons)).

%%%-------------------------------------------------------------------
%%%
%%% Homographic functions of continued fractions.
%%%

:- func apply_ng4(ng4, cf) = cf.

:- func add_integer(cf, integer) = cf.
:- func add_rational(cf, rational) = cf.
:- func mul_integer(cf, integer) = cf.
:- func mul_rational(cf, rational) = cf.
:- func div_integer(cf, integer) = cf.
:- func reciprocal(cf) = cf.

add_integer(CF, I) = apply_ng4({one, I, zero, one}, CF).
add_rational(CF, R) = CF1 :-
  N = (rational.numer(R)),
  D = (rational.denom(R)),
  CF1 = apply_ng4({D, N, zero, D}, CF).
mul_integer(CF, I) = apply_ng4({I, zero, zero, one}, CF).
mul_rational(CF, R) = apply_ng4({numer(R), zero, zero, denom(R)}, CF).
div_integer(CF, I) = apply_ng4({one, zero, zero, I}, CF).
reciprocal(CF) = apply_ng4({zero, one, one, zero}, CF).

apply_ng4({ A1, A, B1, B }, Other_CF) = CF :-
  (if (B1 = zero, B = zero)
   then (CF = cfnil)
   else if (B1 \= zero, B \= zero)
   then (
     % The integer divisions here truncate towards zero. Say "div"
     % instead of "//" to truncate towards negative infinity.
     Q1 = A1 // B1,
     Q = A // B,
     (if (Q1 = Q)
      then (CF = lzlist(delay(Cons)),
            Cons = ((func) = lzcons(Q, ng4_eject_term(A1, A, B1, B,
                                                      Other_CF, Q))))
      else (CF = ng4_absorb_term(A1, A, B1, B, Other_CF)))
   )
   else (CF = ng4_absorb_term(A1, A, B1, B, Other_CF))).

:- func ng4_eject_term(integer, integer, integer, integer, cf,
                       integer) = cf.
ng4_eject_term(A1, A, B1, B, Other_CF, Term) = CF :-
  CF = apply_ng4({ B1, B, A1 - (B1 * Term), A - (B * Term) },
                 Other_CF).

:- func ng4_absorb_term(integer, integer, integer, integer, cf) = cf.
ng4_absorb_term(A1, A, B1, B, Other_CF) = CF :-
  (Other_CF = lzlist(Cell),
   CF = (if (force(Cell) = lzcons(Term, Rest))
         then apply_ng4({ A + (A1 * Term), A1,
                          B + (B1 * Term), B1 },
                        Rest)
         else apply_ng4({ A1, A1, B1, B1 }, cfnil))).


%%%-------------------------------------------------------------------
%%%
%%% cf2string and cf2string_with_max_terms convert a continued
%%% fraction to a printable string.
%%%

:- func cf2string(cf) = string.
:- func cf2string_with_max_terms(cf, integer) = string.

cf2string(CF) = cf2string_with_max_terms(CF, integer(20)).

cf2string_with_max_terms(CF, MaxTerms) = S :-
  S = cf2string_loop(CF, MaxTerms, zero, "[").

:- func cf2string_loop(cf, integer, integer, string) = string.
cf2string_loop(CF, MaxTerms, I, Accum) = S :-
  (CF = lzlist(ValCell),
   force(ValCell) = Cell,
   (if (Cell = lzcons(Term, Tail))
    then (if (I = MaxTerms) then (S = Accum ++ ",...]")
          else ((Separator = (if (I = zero) then ""
                              else if (I = one) then ";"
                              else ",")),
                TermStr = to_string(Term),
                S = cf2string_loop(Tail, MaxTerms, I + one,
                                   Accum ++ Separator ++ TermStr)))
    else (S = Accum ++ "]"))).

%%%-------------------------------------------------------------------

:- pred show(string::in, cf::in, io::di, io::uo) is det.
show(Expression, CF, !IO) :-
  print(Expression, !IO),
  print(" => ", !IO),
  print(cf2string(CF), !IO),
  nl(!IO).

main(!IO) :-
  CF_13_11 = r2cf(rational(13, 11)),
  CF_22_7 = r2cf(rational(22, 7)),

  show("13/11", CF_13_11, !IO),
  show("22/7", CF_22_7, !IO),
  show("sqrt(2)", sqrt2, !IO),

  show("13/11 + 1/2", add_rational(CF_13_11, rational(1, 2)), !IO),
  show("22/7 + 1/2", add_rational(CF_22_7, rational(1, 2)), !IO),
  show("(22/7)/4", div_integer(CF_22_7, integer(4)), !IO),
  show("(22/7)*(1/4)", mul_rational(CF_22_7, rational(1, 4)), !IO),
  show("1/sqrt(2)", reciprocal(sqrt2), !IO),
  show("sqrt(2)/2", div_integer(sqrt2, two), !IO),
  show("sqrt(2)*(1/2)", mul_rational(sqrt2, rational(1, 2)), !IO),

  %% Getting (1 + 1/sqrt(2))/2 in a single step.
  show("(2 + sqrt(2))/4",
       apply_ng4({one, two, zero, integer(4)}, sqrt2),
       !IO),

  %% Different ways to compute the same thing.
  show("(1/sqrt(2) + 1)/2",
       div_integer(add_integer(reciprocal(sqrt2), one),
                   two),
       !IO),
  show("(1/sqrt(2))*(1/2) + 1/2",
       add_rational(mul_rational(reciprocal(sqrt2),
                                 rational(1, 2)),
                    rational(1, 2)),
       !IO),
  show("((sqrt(2)/2 + 1)/4)*2", % Contrived, to get in mul_integer.
       mul_integer(div_integer(add_integer(div_integer(sqrt2, two),
                                           one),
                               integer(4)),
                   two),
       !IO),

  true.

%%%-------------------------------------------------------------------
%%% local variables:
%%% mode: mercury
%%% prolog-indent-width: 2
%%% end:
Output:
$ mmc -m univariate_continued_fraction_task_lazy && ./univariate_continued_fraction_task_lazy
Making Mercury/int3s/univariate_continued_fraction_task_lazy.int3
Making Mercury/ints/univariate_continued_fraction_task_lazy.int
Making Mercury/cs/univariate_continued_fraction_task_lazy.c
Making Mercury/os/univariate_continued_fraction_task_lazy.o
Making univariate_continued_fraction_task_lazy
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
(22/7)*(1/4) => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2)*(1/2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(1/sqrt(2) + 1)/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
(1/sqrt(2))*(1/2) + 1/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
((sqrt(2)/2 + 1)/4)*2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]

Nim

Translation of: Kotlin
import math, rationals, strformat

type
  Rat = Rational[int]
  Ng = tuple[a1, a, b1, b: int]

const NS = 1 // 1   # Non significant value.


iterator r2cf(rat: Rat): int {.closure.} =
  var
    num = rat.num
    den = rat.den
  for count in 1..20:
    let d = num div den
    num = num mod den
    swap num, den
    yield d
    if den == 0: break


iterator d2cf(f: float): int {.closure.} =
  var f = f
  for count in 1..20:
    let d = floor(f)
    let r = f - d
    yield int(d)
    if r == 0: break
    f = 1 / r


iterator root2(dummy: Rat): int {.closure.} =
  yield 1
  for count in 1..20: yield 2


iterator recipRoot2(rat: Rat): int {.closure.} =
  yield 0
  yield 1
  for count in 1..20: yield 2


func ingress(ng: var Ng; n: int) =
  ng.a += ng.a1 * n
  swap ng.a, ng.a1
  ng.b += ng.b1 * n
  swap ng.b, ng.b1


func egress(ng: var Ng): int =
  let n = ng.a div ng.b
  ng.a -= ng.b * n
  swap ng.a, ng.b
  ng.a1 -= ng.b1 * n
  swap ng.a1, ng.b1
  result = n


func needTerm(ng: Ng): bool = ng.b == 0 or ng.b1 == 0 or (ng.a // ng.b != ng.a1 // ng.b1)


func egressDone(ng: var Ng): int =
  if ng.needTerm:
    ng.a = ng.a1
    ng.b = ng.b1
  result = ng.egress()


func done(ng: Ng): bool = ng.b == 0 or ng.b1 == 0


when isMainModule:

  let data = [("[1;5,2] + 1/2        ", (2, 1, 0, 2), 13 // 11, r2cf),
              ("[3;7] + 1/2          ", (2, 1, 0, 2), 22 // 7,  r2cf),
              ("[3;7] divided by 4   ", (1, 0, 0, 4), 22 // 7,  r2cf),
              ("sqrt(2)              ", (0, 1, 1, 0), NS,  recipRoot2),
              ("1 / sqrt(2)          ", (0, 1, 1, 0), NS,  root2),
              ("(1 + sqrt(2)) / 2    ", (1, 1, 0, 2), NS,  root2),
              ("(1 + 1 / sqrt(2)) / 2", (1, 1, 0, 2), NS,  recipRoot2)]

  echo "Produced by NG object:"
  for (str, ng, r, gen) in data:
    stdout.write &"{str} → "
    var op = ng
    for n in gen(r):
      if not op.needTerm: stdout.write &" {op.egress()} "
      op.ingress(n)
    while true:
      stdout.write &" {op.egressDone} "
      if op.done: break
    echo()

  echo "\nProduced by direct calculation:"
  let data2 = [("(1 + sqrt(2)) / 2    ", (1 + sqrt(2.0)) / 2),
               ("(1 + 1 / sqrt(2)) / 2", (1 + 1 / sqrt(2.0)) / 2)]
  for (str, d) in data2:
    stdout.write &"{str} →"
    for n in d2cf(d): stdout.write "  ", n
    echo()
Output:
Produced by NG object:
[1;5,2] + 1/2         →  1  1  2  7 
[3;7] + 1/2           →  3  1  1  1  4 
[3;7] divided by 4    →  0  1  3  1  2 
sqrt(2)               →  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 
1 / sqrt(2)           →  0  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 
(1 + sqrt(2)) / 2     →  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4 
(1 + 1 / sqrt(2)) / 2 →  0  1  5  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  5 

Produced by direct calculation:
(1 + sqrt(2)) / 2     →  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  4  1  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

ObjectIcon

Translation of: ATS
Translation of: Icon

This is essentially the Icon implementation, but with the data and procedures encapsulated in classes.

(The generators are likely to run faster than in recent versions of Arizona Icon, due to a faster implementation of co-expressions. Unicon might run the conventional Icon implementation quickly, however.)

# -*- ObjectIcon -*-

import io

procedure main ()
  local cf_13_11, cf_22_7, cf_sqrt2, cf_1_div_sqrt2

  cf_13_11 := CF_rational (13, 11)
  cf_22_7 := CF_rational (22, 7)
  cf_sqrt2 := CF_sqrt2()
  cf_1_div_sqrt2 := CF_hfunc (0, 1, 1, 0, cf_sqrt2)

  show ("13/11", cf_13_11)
  show ("22/7", cf_22_7)
  show ("sqrt(2)", cf_sqrt2)
  show ("13/11 + 1/2", CF_hfunc (2, 1, 0, 2, cf_13_11))
  show ("22/7 + 1/2", CF_hfunc (2, 1, 0, 2, cf_22_7))
  show ("(22/7)/4", CF_hfunc (1, 0, 0, 4, cf_22_7))
  show ("1/sqrt(2)", cf_1_div_sqrt2)
  show ("(2 + sqrt(2))/4", CF_hfunc (1, 2, 0, 4, cf_sqrt2))
  show ("(1 + 1/sqrt(2))/2", CF_hfunc (1, 1, 0, 2,
                                       cf_1_div_sqrt2))
end

procedure show (expr, cf)
  io.write (expr, " => ", cf.to_string())
end

class CF ()                  # A continued fraction.

  private terminated         # Are there no more terms to memoize?
  private memo               # Memoized terms.
  private generate           # A co-expression to generate more terms.

  public new (gen)
    terminated := &no
    memo := []
    generate := gen
    return
  end

  public get_term (i)
    local j, term

    if *memo <= i then {
      if \terminated then {
        fail
      } else {
        every j := *memo to i do {
          if term := @generate then {
            put (memo, term)
          } else {
            terminated := &yes
            fail
          }
        }
      }
    }
    return memo[i + 1]
  end

  public to_string (max_terms)
    local s, sep, i, done, term

    /max_terms := 20

    s := "["
    sep := 0
    i := 0
    done := &no
    while /done do {
      if i = max_terms then {
        # We have reached the maximum of terms to print. Stick an
        # ellipsis in the notation.
        s ||:= ",...]"
        done := &yes
      } else if term := get_term (i) then {
        # Getting a term succeeded. Include the term.
        s ||:= sep_str (sep) || term
        sep := min (sep + 1, 2)
        i +:= 1
      } else {
        # Getting a term failed. We are done.
        s ||:= "]"
        done := &yes
      }
    }
    return s
  end

  private sep_str (sep)
    return (if sep = 0 then "" else if sep = 1 then ";" else ",")
  end

end                             # class CF

class CF_sqrt2 (CF)             # A continued fraction for sqrt(2).
  public override new ()
    CF.new (create gen ())
    return
  end

  private gen ()
    suspend 1
    repeat suspend 2
  end
end                             # class CF_sqrt2

class CF_rational (CF)   # A continued fraction for a rational number.
  public override new (numerator, denominator)
    CF.new (create gen (numerator, denominator))
    return
  end

  private gen (n, d)
    local q, r

    repeat {
      if d = 0 then fail
      q := n / d
      r := n % d
      n := d
      d := r
      suspend q
    }
  end
end                             # class CF_rational

class CF_hfunc (CF) # A continued fraction for a homographic function
                    # of some other continued fraction.

  public override new (a1, a, b1, b, other_cf)
    CF.new (create gen (a1, a, b1, b, other_cf))
    return
  end

  private gen (a1, a, b1, b, other_cf)
    local a1_tmp, a_tmp, b1_tmp, b_tmp
    local i, term, skip_getting_a_term
    local q1, q

    i := 0
    repeat {
      skip_getting_a_term := &no
      if b1 = b = 0 then {
        fail
      } else if b1 ~= 0 & b ~= 0 then {
        q1 := a1 / b1
        q := a / b
        if q1 = q then {
          a1_tmp := a1
          a_tmp := a
          b1_tmp := b1
          b_tmp := b
          a1 := b1_tmp
          a := b_tmp
          b1 := a1_tmp - (b1_tmp * q)
          b := a_tmp - (b_tmp * q)
          suspend q
          skip_getting_a_term := &yes
        }
      }
      if /skip_getting_a_term then {
        if term := other_cf.get_term (i) then {
          i +:= 1
          a1_tmp := a1
          a_tmp := a
          b1_tmp := b1
          b_tmp := b
          a1 := a_tmp + (a1_tmp * term)
          a := a1_tmp
          b1 := b_tmp + (b1_tmp * term)
          b := b1_tmp
        } else {
          a := a1
          b := b1
        }
      }
    }
  end

end                             # class CF_hfunc
Output:
$ oiscript univariate-continued-fraction-task-OI.icn
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]

OCaml

Translation of: ATS

This implementation memoizes terms of a continued fraction.

module CF =                     (* A continued fraction. *)
  struct
    type record_t =
      {
        terminated : bool;   (* Are there no more terms to memoize? *)
        m : int;             (* The number of memoized terms. *)
        memo : int Array.t;  (* Storage for the memoized terms. *)
        gen : unit -> int option; (* A generator of new terms. *)
      }

    type t = record_t ref

    let make gen =
      ref { terminated = false;
            m = 0;
            memo = Array.make (8) 0;
            gen = gen }

    let get cf i =
      let get_more_terms record needed =
        let rec loop j =
          if j = needed then
            { record with terminated = false; m = needed }
          else
            match record.gen () with
            | None -> { record with terminated = true; m = i }
            | Some term ->
               begin
                 record.memo.(i) <- term;
                 loop (j + 1)
               end
        in
        loop record.m
      in
      let update record needed =
        if record.terminated then
          record
        else if needed <= record.m then
          record
        else if needed <= Array.length record.memo then
          get_more_terms record needed
        else
          (* Provide twice the room that might be needed. *)
          let n1 = needed + needed in
          let memo1 = Array.make (n1) 0 in
          let record =
            begin
              for j = 0 to record.m - 1 do
                memo1.(j) <- record.memo.(j)
              done;
              { record with memo = memo1 }
            end
          in
          get_more_terms record needed
      in
      let record = update !cf (i + 1) in
      begin
        cf := record;
        if i < record.m then
          Some record.memo.(i)
        else
          None
      end

    let to_string ?max_terms:(max_terms = 20) cf =
      let rec loop i sep accum =
        if i = max_terms then
          accum ^ ",...]"
        else
          match get cf i with
          | None -> accum ^ "]"
          | Some term ->
             let sep_str =
               match sep with
               | 0 -> ""
               | 1 -> ";"
               | _ -> ","
             in
             let sep = min (sep + 1) 2 in
             let term_str = string_of_int term in
             let accum = accum ^ sep_str ^ term_str in
             loop (i + 1) sep accum
      in
      loop 0 0 "["

    let to_thunk cf =     (* To use a CF.t as a generator of terms. *)
      let index = ref 0 in
      fun () -> let i = !index in
                begin
                  index := i + 1;
                  get cf i
                end
  end

let cf_sqrt2 =                 (* A continued fraction for sqrt(2). *)
  CF.make (let next_term = ref 1 in
           fun () -> let term = !next_term in
                     begin
                       next_term := 2;
                       Some term
                     end)

let cf_rational n d = (* Make a continued fraction for a rational
                         number. *)
  CF.make (let ratnum = ref (n, d) in
           fun () -> let (n, d) = !ratnum in
                     if d = 0 then
                       None
                     else
                       let q = n / d and r = n mod d in
                       begin
                         ratnum := (d, r);
                         Some q
                       end)

let cf_hfunc (a1, a, b1, b) other_cf =
  let gen = CF.to_thunk other_cf in
  let state = ref (a1, a, b1, b, gen) in
  let hgen () =
    let rec loop () =
      let (a1, a, b1, b, gen) = !state in
      let absorb_term () =
        match gen () with
        | None -> state := (a1, a1, b1, b1, gen)
        | Some term -> state := (a + (a1 * term), a1,
                                 b + (b1 * term), b1, 
                                 gen)
      in
      if b1 = 0 && b = 0 then
        None
      else if b1 <> 0 && b <> 0 then
        let q1 = a1 / b1 and q = a / b in
        if q1 = q then
          begin
            state := (b1, b, a1 - (b1 * q), a - (b * q), gen);
            Some q
          end
        else
          begin
            absorb_term ();
            loop ()
          end
      else
        begin
          absorb_term ();
          loop ()
        end
    in
    loop ()
  in
  CF.make hgen

;;

let show expr cf =
  begin
    print_string expr;
    print_string " => ";
    print_string (CF.to_string cf);
    print_newline ()
  end ;;

let hf_cf_add_1_2 = (2, 1, 0, 2) ;;
let hf_cf_add_1 = (1, 1, 0, 1) ;;
let hf_cf_div_2 = (1, 0, 0, 2) ;;
let hf_cf_div_4 = (1, 0, 0, 4) ;;
let hf_1_div_cf = (0, 1, 1, 0) ;;

let cf_13_11 = cf_rational 13 11 ;;
let cf_22_7 = cf_rational 22 7 ;;
let cf_1_div_sqrt2 = cf_hfunc hf_1_div_cf cf_sqrt2 ;;

show "13/11" cf_13_11 ;;
show "22/7" cf_22_7 ;;
show "sqrt(2)" cf_sqrt2 ;;
show "13/11 + 1/2" (cf_hfunc hf_cf_add_1_2 cf_13_11) ;;
show "22/7 + 1/2" (cf_hfunc hf_cf_add_1_2 cf_22_7) ;;
show "(22/7)/4" (cf_hfunc hf_cf_div_4 cf_22_7) ;;
show "1/sqrt(2)" cf_1_div_sqrt2 ;;
show "(2 + sqrt(2))/4" (cf_hfunc (1, 2, 0, 4) cf_sqrt2) ;;

(* Demonstrate a chain of operations. *)
show "(1 + 1/sqrt(2))/2" (cf_1_div_sqrt2
                          |> cf_hfunc hf_cf_add_1
                          |> cf_hfunc hf_cf_div_2) ;;

(* Demonstrate a slightly longer chain of operations. *)
show "((sqrt(2)/2) + 1)/2" (cf_sqrt2
                            |> cf_hfunc hf_cf_div_2
                            |> cf_hfunc hf_cf_add_1
                            |> cf_hfunc hf_cf_div_2) ;;
Output:
$ ocaml univariate_continued_fraction_task.ml
13/11 => [1;5,2]
22/7 => [3;7]
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/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7)/4 => [0;1,3,1,2]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(2 + sqrt(2))/4 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,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,...]
((sqrt(2)/2) + 1)/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]

Phix

Library: Phix/Class
Library: Phix/mpfr

Self-contained. The supporting cast of r2cf(), cf2s(), cf2r() and d2cf() ended up being more code than the task itself.

requires("0.8.2")

class baby_matrix

  integer a1, a, b1, b

  --
  -- used by apply_baby_matrix to yield (a1*cf+a)/(b1*cf+b)
  --
  -- examples: (a1 a  b1 b) => above, simplified:
  -- ========   =  =  =  =
  --           {2, 1, 0, 2} => (2*cf+1)/2, aka cf+1/2
  --           {1, 0, 0, 4} => cf/4
  --           {1, 0, 0, 1} => cf/1, aka cf
  --           {0, 1, 1, 0} => 1/cf
  --           {1, 1, 0, 2} => (cf+1)/2
  --

  function need_term()
    return b==0 or b1==0 or ((a/b)!=(a1/b1))
  end function

  function next_term()
    integer n = floor(a/b)
    {a1,a,b1,b} = {b1,b,a1-b1*n,a-b*n}
    return n
  end function

  procedure in_term(object n={})
    if integer(n) then
        {a1,a,b,b1} = {a+a1*n,a1,b1,b+b1*n}
    else
        {a,b} = {a1,b1}
    end if
  end procedure

  function done()
    return b=0 and b1=0
  end function

end class

function apply_baby_matrix(sequence m, cf)
--
--  for m of integer {a1,a,b1,b}, return (a1*cf+a)/(b1*cf+b):
--
    baby_matrix bm = new(m)
    sequence res = {}
    for i=1 to length(cf) do
        if not bm.need_term() then
            res &= bm.next_term()
        end if  
        bm.in_term(cf[i])
    end for
    while true do   
        if bm.need_term() then
            bm.in_term()
        end if
        res &= bm.next_term()
        if bm.done() then exit end if
    end while
    return res
end function

function r2cf(sequence rat, integer count=20)
    sequence s = {}
    atom {num,den} = rat
    while den!=0 and length(s)<count do
        s &= trunc(num/den)
        {num,den} = {den,num-s[$]*den}
    end while
    return s
end function

function root2(integer count=20)
    return {1}&repeat(2,count-1)
end function
 
function recip_root2(integer count=20)
    return {0,1}&repeat(2,count-2)
end function

function cf2s(sequence cf)
    sequence s = join(apply(cf,sprint),",") -- eg "1,5,2"
    return "["&substitute(s,",",";",1)&"]"  -- => "[1;5,2]"
end function

include mpfr.e

function cf2r(sequence cf)
    mpq res = mpq_init(), -- 0/1
        cfn = mpq_init()
    for n=length(cf) to 1 by -1 do
        mpq_set_si(cfn,cf[n])
        mpq_add(res,res,cfn)
        if n=1 then exit end if
        mpq_inv(res,res)
    end for
    mpz num = mpz_init(),
        den = mpz_init()
    mpq_get_num(num,res)
    mpq_get_den(den,res)
    mpfr x = mpfr_init()
    mpfr_set_q(x,res)
    string xs = mpfr_sprintf("%.15Rf",x),
           ns = mpz_get_str(num),
           ds = mpz_get_str(den),
            s = sprintf("%s (%s/%s)",{xs,ns,ds})
    return s
end function

function d2cf(atom d, integer count=20)
    string res = "["
    integer sep = ';'
    while count do
        atom div = floor(d),
             rem = d - div
        res &= sprintf("%d%c",{div,sep})
        if rem==0 then exit end if
        d = 1/rem
        count -= 1
        sep = ','
    end while
    res[$] = ']'
    return res
end function

constant tests = {
    {"[1;5,2] + 1/2  ", {2, 1, 0, 2}, r2cf({13,11}), 37/22},
    {"[3;7] + 1/2    ", {2, 1, 0, 2}, r2cf({22, 7}), 3+1/7+1/2},
    {"[3;7] / 4      ", {1, 0, 0, 4}, r2cf({22, 7}), (3+1/7)/4},
    {"sqrt(2)        ", {1, 0, 0, 1}, root2(),       sqrt(2)},
    {"sqrt(2) (inv)  ", {0, 1, 1, 0}, recip_root2(), 1/(1/sqrt(2))},
    {"1/sqrt(2)      ", {1, 0, 0, 1}, recip_root2(), 1/sqrt(2)},
    {"1/sqrt(2) (inv)", {0, 1, 1, 0}, root2(),       1/sqrt(2)},
    {"(1+sqrt(2))/2  ", {1, 1, 0, 2}, root2(),       (1+sqrt(2))/2},
    {"(1+1/sqrt(2))/2", {1, 1, 0, 2}, recip_root2(), (1+1/sqrt(2))/2}}

for i=1 to length(tests) do
    {string str, sequence bm, sequence cf, atom eres} = tests[i]
    sequence res = apply_baby_matrix(bm, cf)
    printf(1,"%s ->  %s --> %s\n",{str,cf2s(res),cf2r(res)})
    printf(1,"           direct:  %s ==> %.15f\n",{d2cf(eres,length(res)),eres})
end for
Output:
[1;5,2] + 1/2   ->  [1;1,2,7] --> 1.681818181818182 (37/22)
           direct:  [1;1,2,6] ==> 1.681818181818182
[3;7] + 1/2     ->  [3;1,1,1,4] --> 3.642857142857143 (51/14)
           direct:  [3;1,1,1,3] ==> 3.642857142857143
[3;7] / 4       ->  [0;1,3,1,2] --> 0.785714285714286 (11/14)
           direct:  [0;1,3,1,2] ==> 0.785714285714286
sqrt(2)         ->  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] --> 1.414213562373096 (22619537/15994428)
           direct:  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] ==> 1.414213562373095
sqrt(2) (inv)   ->  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] --> 1.414213562373087 (9369319/6625109)
           direct:  [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] ==> 1.414213562373095
1/sqrt(2)       ->  [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] --> 0.707106781186552 (6625109/9369319)
           direct:  [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] ==> 0.707106781186547
1/sqrt(2) (inv) ->  [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] --> 0.707106781186547 (15994428/22619537)
           direct:  [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] ==> 0.707106781186547
(1+sqrt(2))/2   ->  [1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4] --> 1.207106781186548 (38613965/31988856)
           direct:  [1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4] ==> 1.207106781186547
(1+1/sqrt(2))/2 ->  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,5] --> 0.853553390593276 (7997214/9369319)
           direct:  [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4] ==> 0.853553390593274

The last digits of direct in the first two tests match on 64-bit, ie ,7] and ,4], plus 6/7/8 end in 548.

Python

Translation of: Ruby

Python: NG

class NG:
  def __init__(self, a1, a, b1, b):
    self.a1, self.a, self.b1, self.b = a1, a, b1, b

  def ingress(self, n):
    self.a, self.a1 = self.a1, self.a + self.a1 * n
    self.b, self.b1 = self.b1, self.b + self.b1 * n

  @property
  def needterm(self):
    return (self.b == 0 or self.b1 == 0) or not self.a//self.b == self.a1//self.b1

  @property
  def egress(self):
    n = self.a // self.b
    self.a,  self.b  = self.b,  self.a  - self.b  * n
    self.a1, self.b1 = self.b1, self.a1 - self.b1 * n
    return n

  @property
  def egress_done(self):
    if self.needterm: self.a, self.b = self.a1, self.b1
    return self.egress

  @property
  def done(self):
    return self.b == 0 and self.b1 == 0

Python: Testing

Uses r2cf method from here.

data = [["[1;5,2] + 1/2",      [2,1,0,2], [13,11]],
        ["[3;7] + 1/2",        [2,1,0,2], [22, 7]],
        ["[3;7] divided by 4", [1,0,0,4], [22, 7]]]

for string, ng, r in data:
  print( "%-20s->" % string, end='' )
  op = NG(*ng)
  for n in r2cf(*r):
    if not op.needterm: print( " %r" % op.egress, end='' )
    op.ingress(n)
  while True:
    print( " %r" % op.egress_done, end='' )
    if op.done: break
  print()
Output:
[1;5,2] + 1/2       -> 1 1 2 7
[3;7] + 1/2         -> 3 1 1 1 4
[3;7] divided by 4  -> 0 1 3 1 2

Racket

Translation of: Python
Translation of: C++

Main part of the NG-baby matrices. They are implemented as mutable structs.

#lang racket/base

(struct ng (a1 a b1 b) #:transparent #:mutable)
 
(define (ng-ingress! v t)
  (define a (ng-a v))
  (define a1 (ng-a1 v))
  (define b (ng-b v))
  (define b1 (ng-b1 v))
  (set-ng-a! v a1)
  (set-ng-a1! v (+ a (* a1 t)))
  (set-ng-b! v b1)
  (set-ng-b1! v (+ b (* b1 t))))
 
(define (ng-needterm? v)
  (or (zero? (ng-b v)) 
      (zero? (ng-b1 v)) 
      (not (= (quotient (ng-a v) (ng-b v)) (quotient (ng-a1 v) (ng-b1 v))))))
 
(define (ng-egress! v)
  (define t (quotient (ng-a v) (ng-b v)))
  (define a (ng-a v))
  (define a1 (ng-a1 v))
  (define b (ng-b v))
  (define b1 (ng-b1 v))
  (set-ng-a! v b)
  (set-ng-a1! v b1)
  (set-ng-b! v (- a (* b t)))
  (set-ng-b1! v (- a1 (* b1 t)))
  t)
 
(define (ng-infty! v)
  (when (ng-needterm? v)
    (set-ng-a! v (ng-a1 v))
    (set-ng-b! v (ng-b1 v))))
 
(define (ng-done? v)
  (and (zero? (ng-b v)) (zero? (ng-b1 v))))

Auxiliary functions to create producers of well known continued fractions. The function rational->cf is copied from r2cf task.

(define ((rational->cf n d))
  (and (not (zero? d))
       (let-values ([(q r) (quotient/remainder n d)])
         (set! n d)
         (set! d r)
         q)))

(define (sqrt2->cf)
  (define first? #t)
  (lambda ()
    (if first?
        (begin 
          (set! first? #f)
          1)
        2)))

The function combine-ng-cf->cf combines a ng-matrix and a cf- producer and creates a cf-producer. The cf-producers can represent infinitely long continued fractions. The function cf-showln shows the first coefficients of a continued fraction represented in a cf-producer.

(define (combine-ng-cf->cf ng cf)
  (define empty-producer? #f)
  (lambda ()
    (let loop ()
      (cond 
        [(not empty-producer?) (define t (cf))
                               (cond 
                                   [t (ng-ingress! ng t)
                                      (if (ng-needterm? ng)
                                          (loop)
                                          (ng-egress! ng))]
                                   [else (set! empty-producer? #t)
                                         (loop)])]
        [(ng-done? ng) #f]
        [(ng-needterm? ng) (ng-infty! ng) 
                           (loop)]
        [else (ng-egress! ng)]))))

(define (cf-showln cf n)
  (for ([i (in-range n)])
    (define val (cf))
    (when val
      (printf " ~a" val)))
  (when (cf)
    (printf " ..."))
  (printf "~n"))

Some test

(display "[1;5,2] + 1/2 ->")
(cf-showln (combine-ng-cf->cf (ng 2 1 0 2) (rational->cf 13 11)) 20)

(display "[3;7] + 1/2 ->")
(cf-showln (combine-ng-cf->cf (ng 2 1 0 2) (rational->cf 22 7)) 20)

(display "[3;7] / 4 ->")
(cf-showln (combine-ng-cf->cf (ng 1 0 0 4) (rational->cf 22 7)) 20)

(display "sqrt(2)/2 ->")
(cf-showln (combine-ng-cf->cf (ng 1 0 0 2) (sqrt2->cf)) 20)

(display "1/sqrt(2) ->")
(cf-showln (combine-ng-cf->cf (ng 0 1 1 0) (sqrt2->cf)) 20)

(display "(1+sqrt(2))/2 ->")
(cf-showln (combine-ng-cf->cf (ng 1 1 0 2) (sqrt2->cf)) 20)

Sample output:

[1;5,2] + 1/2 -> 1 1 2 7
[3;7] + 1/2 -> 3 1 1 1 4
[3;7] / 4 -> 0 1 3 1 2
sqrt(2)/2 -> 0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
1/sqrt(2) -> 0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
(1+sqrt(2))/2 -> 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 ...

Raku

(formerly Perl 6)

Works with: Rakudo version 2020.08.1

All the important stuff takes place in the NG object. Everything else is helper subs for testing and display. The NG object is capable of working with infinitely long continued fractions, but displaying them can be problematic. You can pass in a limit to the apply method to get a fixed maximum number of terms though. See the last example: 100 terms from the infinite cf (1+√2)/2 and its Rational representation.

class NG {
    has ( $!a1, $!a, $!b1, $!b );
    submethod BUILD ( :$!a1, :$!a, :$!b1, :$!b ) { }

    # Public methods
    method new( $a1, $a, $b1, $b ) { self.bless( :$a1, :$a, :$b1, :$b ) }
    method apply(@cf, :$limit = Inf) {
        (gather {
            map { take self!extract unless self!needterm; self!inject($_) }, @cf;
            take self!drain until self!done;
        })[ ^ $limit ]
    }

    # Private methods
    method !inject ($n) {
        sub xform($n, $x, $y) { $x, $n * $x + $y }
        ( $!a, $!a1 ) = xform( $n, $!a1, $!a );
        ( $!b, $!b1 ) = xform( $n, $!b1, $!b );
    }
    method !extract {
        sub xform($n, $x, $y) { $y, $x - $y * $n }
        my $n = $!a div $!b;
        ($!a,  $!b ) = xform( $n, $!a,  $!b  );
        ($!a1, $!b1) = xform( $n, $!a1, $!b1 );
        $n
    }
    method !drain { $!a = $!a1, $!b = $!b1 if self!needterm; self!extract }
    method !needterm { so [||] !$!b, !$!b1, $!a/$!b != $!a1/$!b1 }
    method !done { not [||] $!b, $!b1 }
}

sub r2cf(Rat $x is copy) { # Rational to continued fraction
    gather loop {
	$x -= take $x.floor;
	last if !$x;
	$x = 1 / $x;
    }
}

sub cf2r(@a) { # continued fraction to Rational
    my $x = @a[* - 1]; # Use FatRats for arbitrary precision
    $x = ( @a[$_- 1] + 1 / $x ).FatRat for reverse 1 ..^ @a;
    $x
}

sub ppcf(@cf) { # format continued fraction for pretty printing 
    "[{ @cf.join(',').subst(',',';') }]"
}

sub pprat($a) { # format Rational for pretty printing
   # Use FatRats for arbitrary precision
   $a.FatRat.denominator == 1 ?? $a !! $a.FatRat.nude.join('/')
}

sub test_NG ($rat, @ng, $op) { 
    my @cf = $rat.Rat(1e-18).&r2cf;
    my @op = NG.new( |@ng ).apply( @cf );
    say $rat.raku, ' as a cf: ', @cf.&ppcf, " $op = ",
        @op.&ppcf, "\tor ", @op.&cf2r.&pprat, "\n";
}

# Testing
test_NG(|$_) for (
    [ 13/11, [<2 1 0 2>], '+ 1/2 '    ],
    [ 22/7,  [<2 1 0 2>], '+ 1/2    ' ],
    [ 22/7,  [<1 0 0 4>], '/ 4      ' ],
    [ 22/7,  [<7 0 0 22>], '* 7/22   ' ],
    [ 2**.5, [<1 1 0 2>], "\n(1+√2)/2 (approximately)" ]
);

say '100 terms of (1+√2)/2 as a continued fraction and as a rational value:';

my @continued-fraction = NG.new( 1,1,0,2 ).apply( (lazy flat 1, 2 xx * ), limit => 100 );
say @continued-fraction.&ppcf.comb(/ . ** 1..80/).join("\n");
say @continued-fraction.&cf2r.&pprat;
Output:
<13/11> as a cf: [1;5,2] + 1/2  = [1;1,2,7]	or 37/22

<22/7> as a cf: [3;7] + 1/2     = [3;1,1,1,4]	or 51/14

<22/7> as a cf: [3;7] / 4       = [0;1,3,1,2]	or 11/14

<22/7> as a cf: [3;7] * 7/22    = [1]	or 1

1.4142135623731e0 as a cf: [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] 
(1+√2)/2 (approximately) = [1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4]	or 225058681/186444716

100 terms of (1+√2)/2 and its rational value
[1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4
,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4
,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4]
161733217200188571081311986634082331709/133984184101103275326877813426364627544

Ruby

NG

# I define a class to implement baby NG
class NG
  def initialize(a1, a, b1, b)
    @a1, @a, @b1, @b = a1, a, b1, b
  end
  def ingress(n)
    @a, @a1 = @a1, @a + @a1 * n
    @b, @b1 = @b1, @b + @b1 * n
  end
  def needterm?
    return true if @b == 0 or @b1 == 0
    return true unless @a/@b == @a1/@b1
    false
  end
  def egress
    n = @a / @b
    @a,  @b  = @b,  @a  - @b  * n
    @a1, @b1 = @b1, @a1 - @b1 * n
    n
  end
  def egress_done
    @a, @b = @a1, @b1 if needterm?
    egress
  end
  def done?
    @b == 0 and @b1 == 0
  end
end

Testing

Uses r2cf method from here.

data = [["[1;5,2] + 1/2",      [2,1,0,2], [13,11]],
        ["[3;7] + 1/2",        [2,1,0,2], [22, 7]],
        ["[3;7] divided by 4", [1,0,0,4], [22, 7]]]

data.each do |str, ng, r|
  printf "%-20s->", str
  op = NG.new(*ng)
  r2cf(*r) do |n|
    print " #{op.egress}" unless op.needterm?
    op.ingress(n)
  end
  print " #{op.egress_done}" until op.done?
  puts
end
Output:
[1;5,2] + 1/2       -> 1 1 2 7
[3;7] + 1/2         -> 3 1 1 1 4
[3;7] divided by 4  -> 0 1 3 1 2

Scheme

Translated from Racket

Translation of: Racket
Works with: Gauche Scheme version 0.9.12
Works with: CHICKEN Scheme version 5.3.0

For CHICKEN Scheme you need the r7rs egg.

I generate continued fractions differently from how the Racket does. I use an incomprehensible jumble of "calls with the current continuation" that lets you return values but keep going. Thus you can write your generators as loops or recursions. The method would work in Racket as well.

(Use of some procedure less general than call/cc might improve performance in Racket. In CHICKEN, call/cc itself should be efficient.)

;;;-------------------------------------------------------------------
;;;
;;; For R7RS Scheme, translated from the Racket.
;;;

(cond-expand
  (r7rs)
  (chicken (import r7rs)))             ; CHICKEN is not natively R7RS.

;;;-------------------------------------------------------------------
;;;
;;; A partial implementation of Icon-style co-expressions.
;;;
;;; This limited form does not implement co-expressions that receive
;;; inputs.
;;;

(define-library (suspendable-procedures)

  (export suspend)
  (export make-generator-procedure)

  (import (scheme base))

  (begin

    (define *suspend* (make-parameter (lambda (x) x)))
    (define (suspend v) ((*suspend*) v))

    (define (make-generator-procedure thunk)
      ;; This is for making a suspendable procedure that takes no
      ;; arguments when resumed. The result is a simple generator of
      ;; values.
      (define (next-run return)
        (define (my-suspend v)
          (set! return (call/cc (lambda (resumption-point)
                                  (set! next-run resumption-point)
                                  (return v)))))
        (parameterize ((*suspend* my-suspend))
          (suspend (thunk))))
      (lambda () (call/cc next-run)))

    )) ;; end library

;;;-------------------------------------------------------------------
;;;
;;; Let us decide how we wish to do integer division.
;;;

(define-library (division-procedures)
  (export div rem divrem)
  (import (scheme base))
  (begin
    (define div floor-quotient)
    (define rem floor-remainder)
    (define divrem floor/)))

;;;-------------------------------------------------------------------
;;;
;;; The Main part of the baby-NG matrices. They are implemented as
;;; R7RS (SRFI-9) records, made to look like the Racket structs.
;;;

(define-library (baby-ng-matrices)

  (export ng ng?
          ng-a1 set-ng-a1!
          ng-a set-ng-a!
          ng-b1 set-ng-b1!
          ng-b set-ng-b!)
  (export ng-ingress!
          ng-needterm?
          ng-egress!
          ng-infty!
          ng-done?)

  (import (scheme base))
  (import (division-procedures))

  (begin

    (define-record-type <ng>
      (ng a1 a b1 b)
      ng?
      (a1 ng-a1 set-ng-a1!)
      (a ng-a set-ng-a!)
      (b1 ng-b1 set-ng-b1!)
      (b ng-b set-ng-b!))

    (define (ng-ingress! v t)
      (define a (ng-a v))
      (define a1 (ng-a1 v))
      (define b (ng-b v))
      (define b1 (ng-b1 v))
      (set-ng-a! v a1)
      (set-ng-a1! v (+ a (* a1 t)))
      (set-ng-b! v b1)
      (set-ng-b1! v (+ b (* b1 t))))
    
    (define (ng-needterm? v)
      (or (zero? (ng-b v)) 
          (zero? (ng-b1 v)) 
          (not (= (div (ng-a v) (ng-b v))
                  (div (ng-a1 v) (ng-b1 v))))))
    
    (define (ng-egress! v)
      (define t (div (ng-a v) (ng-b v)))
      (define a (ng-a v))
      (define a1 (ng-a1 v))
      (define b (ng-b v))
      (define b1 (ng-b1 v))
      (set-ng-a! v b)
      (set-ng-a1! v b1)
      (set-ng-b! v (- a (* b t)))
      (set-ng-b1! v (- a1 (* b1 t)))
      t)
    
    (define (ng-infty! v)
      (when (ng-needterm? v)
        (set-ng-a! v (ng-a1 v))
        (set-ng-b! v (ng-b1 v))))
    
    (define (ng-done? v)
      (and (zero? (ng-b v))
           (zero? (ng-b1 v))))

    )) ;; end library

;;;-------------------------------------------------------------------
;;;
;;; Procedures to create generators of continued fractions. (The
;;; Racket implementations could have been adapted, but I like to use
;;; my suspendable-procedures library.)
;;;

(define-library (cf-generators)

  (export make-generator:rational->cf
          make-generator:sqrt2->cf
          make-generator:apply-baby-ng)

  (import (scheme base))
  (import (baby-ng-matrices))
  (import (suspendable-procedures))
  (import (division-procedures))

  (begin

    ;; Generate n/d.
    (define (make-generator:rational->cf n d)
      (make-generator-procedure
       (lambda ()
         (let loop ((n n)
                    (d d))
           (if (zero? d)
               (begin
                 ;; One might reasonably (suspend +inf.0) instead of
                 ;; (suspend #f)
                 (suspend #f)
                 (loop n d))
               (let-values (((q r) (divrem n d)))
                 (suspend q)
                 (loop d r)))))))

    ;; Generate sqrt(2).
    (define (make-generator:sqrt2->cf)
      (make-generator-procedure
       (lambda ()
         (suspend 1)
         (let loop ()
           (suspend 2)
           (loop)))))

    ;; Apply a baby NG to a generator, resulting in a new generator.
    (define (make-generator:apply-baby-ng ng gen)
      (make-generator-procedure
       (lambda ()
         (let loop ()
           (let ((t (gen)))
             (when t
               (ng-ingress! ng t)
               (unless (ng-needterm? ng)
                 (suspend (ng-egress! ng)))
               (loop))))
         (let loop ()
           (cond ((ng-done? ng)
                  (suspend #f)
                  (loop))
                 ((ng-needterm? ng)
                  (ng-infty! ng)
                  (loop))
                 (else
                  (suspend (ng-egress! ng))
                  (loop)))))))

    )) ;; end library

;;;-------------------------------------------------------------------
;;;
;;; Demo.
;;;

(define-library (demonstration)
  (export demonstration)

  (import (scheme base))
  (import (scheme cxr))
  (import (scheme write))
  (import (baby-ng-matrices))
  (import (cf-generators))

  (begin

    (define (display-cf max-digits)
      (lambda (gen)
        (let loop ((i 0)
                   (sep "["))
          (if (= i max-digits)
              (display ",...")
              (let ((digit (gen)))
                (when digit
                  (display sep)
                  (display digit)
                  (loop (+ i 1) (if (string=? sep "[") ";" ","))))))
        (display "]")))

    (define demonstration-instances
      (let ((rat make-generator:rational->cf)
            (sr2 make-generator:sqrt2->cf))
        `(("[1;5,2] + 1/2" ,(ng 2 1 0 2) ,(rat 13 11))
          ("[3;7] + 1/2" ,(ng 2 1 0 2) ,(rat 22 7))
          ("[3;7] / 4" ,(ng 1 0 0 4) ,(rat 22 7))
          ("sqrt(2)/2", (ng 1 0 0 2) ,(sr2))
          ("1/sqrt(2)" ,(ng 0 1 1 0) ,(sr2))
          ("(1+sqrt(2))/2" ,(ng 1 1 0 2) ,(sr2))
          ("(2+sqrt(2))/4 = (1+1/sqrt(2))/2" ,(ng 1 2 0 4) ,(sr2)))))

    (define (demonstration max-digits)
      (define dsply (display-cf max-digits))
      (do ((p demonstration-instances (cdr p)))
          ((null? p))
        (let ((expr-string (caar p))
              (baby-ng (cadar p))
              (gen (caddar p)))
          (display expr-string)
          (display " => ")
          (dsply (make-generator:apply-baby-ng baby-ng gen))
          (newline))))

    )) ;; end library

(import (demonstration))
(demonstration 20)

;;;-------------------------------------------------------------------
Output:
$ gosh single-continued-fraction-task.scm
[1;5,2] + 1/2 => [1;1,2,7]
[3;7] + 1/2 => [3;1,1,1,4]
[3;7] / 4 => [0;1,3,1,2]
sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
(1+sqrt(2))/2 => [1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,...]
(2+sqrt(2))/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,...]

Translated from ATS

Translation of: ATS
Works with: Gauche Scheme version 0.9.12
Works with: CHICKEN Scheme version 5.3.0

For CHICKEN Scheme you need the r7rs egg.

This implementation memoizes terms of a continued fraction.

(cond-expand
  (r7rs)
  (chicken (import (r7rs))))

(define-library (continued-fraction)

  (export make-continued-fraction
          continued-fraction?
          continued-fraction-ref
          continued-fraction->thunk)
  (export continued-fraction->string
          continued-fraction-max-terms)

  (import (scheme base)
          (scheme case-lambda))

  (begin

    (define-record-type <cf-record>
      ;; terminated? -- are these all the terms there are?
      ;; m           -- how many terms are memoized so far?
      ;; memo        -- where terms are memoized.
      ;; gen         -- a thunk that generates terms.
      (cf-record terminated? m memo gen)
      cf-record?
      (terminated? cf-record-terminated?
                   set-cf-record-terminated?!)
      (m cf-record-m set-cf-record-m!)
      (memo cf-record-memo set-cf-record-memo!)
      (gen cf-record-gen set-cf-record-gen!))

    (define cf-record-memo-start-size 8)

    (define (make-continued-fraction gen)
      (cf-record #f 0 (make-vector cf-record-memo-start-size) gen))

    (define continued-fraction? cf-record?)

    ;; The following is an updating operation, but nevertheless I
    ;; leave out the "!" from the name.
    (define (continued-fraction-ref cf i)
      (cf-update! cf (+ i 1))
      (and (< i (cf-record-m cf))
           (vector-ref (cf-record-memo cf) i)))

    (define (cf-get-more-terms! cf needed)
      (define (loop i)
        (if (= i needed)
            (begin
              (set-cf-record-terminated?! cf #f)
              (set-cf-record-m! cf needed))
            (let ((term ((cf-record-gen cf))))
              (if term
                  (begin
                    (vector-set! (cf-record-memo cf) i term)
                    (loop (+ i 1)))
                  (begin
                    (set-cf-record-terminated?! cf #t)
                    (set-cf-record-m! cf i))))))
      (loop (cf-record-m cf)))

    (define (cf-update! cf needed)
      (cond ((cf-record-terminated? cf) (begin))
            ((<= needed (cf-record-m cf)) (begin))
            ((<= needed (vector-length (cf-record-memo cf)))
             (cf-get-more-terms! cf needed))
            (else
             ;; Provide twice the room that might be needed.
             (let* ((n1 (+ needed needed))
                    (memo1 (make-vector n1)))
               (vector-copy! memo1 0 (cf-record-memo cf))
               (set-cf-record-memo! cf memo1)
               (cf-get-more-terms! cf needed)))))

    (define (continued-fraction->thunk cf)
      ;; Make a generator from a continued fraction.
      (define i 0)
      (lambda ()
        (let ((term (continued-fraction-ref cf i)))
          (set! i (+ i 1))
          term)))

    (define continued-fraction-max-terms (make-parameter 20))

    ;; The following is an updating operation, but nevertheless I
    ;; leave out the "!" from the name.
    (define continued-fraction->string
      (case-lambda
        ((cf) (continued-fraction->string
               cf (continued-fraction-max-terms)))
        ((cf max-terms)
         (let loop ((i 0)
                    (sep 0)
                    (accum "["))
           (if (= i max-terms)
               (string-append accum ",...]")
               (let ((term (continued-fraction-ref cf i)))
                 (if (not term)
                     (string-append accum "]")
                     (let* ((term-str (number->string term))
                            (sep-str (case sep
                                       ((0) "")
                                       ((1) ";")
                                       ((2) ",")))
                            (accum (string-append accum sep-str
                                                  term-str))
                            (sep (min (+ sep 1) 2)))
                       (loop (+ i 1) sep accum)))))))))

    )) ;; end library (continued-fraction)

(define-library (number->continued-fraction)

  (export number->continued-fraction)

  (import (scheme base))
  (import (continued-fraction))

  (begin

    (define (number->continued-fraction x)
      ;; This algorithm works directly with exact rationals, rather
      ;; than numerator and denominator separately.
      (unless (real? x)
        (error "number->continued-fraction: argument must be real" x))
      (let ((ratnum (exact x))
            (terminated? #f))
        (make-continued-fraction
         (lambda ()
           (and (not terminated?)
                (let* ((q (floor ratnum))
                       (diff (- ratnum q)))
                  (if (zero? diff)
                      (set! terminated? #t)
                      (set! ratnum (/ diff)))
                  q))))))

    )) ;; end library (number->continued-fraction)

(define-library (homographic-function)

  (export make-homographic-function
          homographic-function?
          homographic-function-ref
          homographic-function-set!
          homographic-function-copy
          apply-homographic-function
          make-homographic-function-operator)

  (import (scheme base)
          (scheme case-lambda))
  (import (continued-fraction))

  (begin

    (define-record-type <homographic-function>
      (make-homographic-function a1 a b1 b)
      homographic-function?
      (a1 homographic-function-a1 set-homographic-function-a1!)
      (a homographic-function-a set-homographic-function-a!)
      (b1 homographic-function-b1 set-homographic-function-b1!)
      (b homographic-function-b set-homographic-function-b!))

    (define (homographic-function-ref hfunc i)
      (case i
        ((0) (homographic-function-a1 hfunc))
        ((1) (homographic-function-a hfunc))
        ((2) (homographic-function-b1 hfunc))
        ((3) (homographic-function-b hfunc))
        (else
         (error "homographic-function-ref: index out of range" i))))

    (define (homographic-function-set! hfunc i x)
      (case i
        ((0) (set-homographic-function-a1! hfunc x))
        ((1) (set-homographic-function-a! hfunc x))
        ((2) (set-homographic-function-b1! hfunc x))
        ((3) (set-homographic-function-b! hfunc x))
        (else
         (error "homographic-function-set!: index out of range" i))))

    (define (homographic-function-copy hfunc)
      (make-homographic-function (homographic-function-ref hfunc 0)
                                 (homographic-function-ref hfunc 1)
                                 (homographic-function-ref hfunc 2)
                                 (homographic-function-ref hfunc 3)))

    (define (apply-homographic-function hfunc cf)
      (define gen (continued-fraction->thunk cf))
      (define state (homographic-function-copy hfunc))
      (make-continued-fraction
       (lambda ()
         (let loop ()
           (let ((a1 (homographic-function-ref state 0))
                 (a (homographic-function-ref state 1))
                 (b1 (homographic-function-ref state 2))
                 (b (homographic-function-ref state 3)))
             (define (take-term)
               (let ((term (gen)))
                 (if term
                     (set! state
                       (make-homographic-function
                        (+ a (* a1 term)) a1 (+ b (* b1 term)) b1))
                     (begin
                       (homographic-function-set! state 1 a1)
                       (homographic-function-set! state 3 b1)))))
             (cond
              ((and (zero? b1) (zero? b)) #f)
              ((and (not (zero? b1)) (not (zero? b)))
               (let ((q1 (floor-quotient a1 b1))
                     (q (floor-quotient a b)))
                 (if (= q1 q)
                     (begin
                       (set! state
                         (make-homographic-function
                          b1 b (- a1 (* b1 q)) (- a (* b q))))
                       q)
                     (begin
                       (take-term)
                       (loop)))))
              (else
               (take-term)
               (loop))))))))

    (define make-homographic-function-operator
      (case-lambda
        ((hfunc) (lambda (cf)
                   (apply-homographic-function hfunc cf)))
        ((a1 a b1 b) (make-homographic-function-operator
                      (make-homographic-function a1 a b1 b)))))

    )) ;; end library (number->continued-fraction)

(define-library (demonstration)

  (export run-demonstration)

  (import (scheme base)
          (scheme write))
  (import (continued-fraction)
          (number->continued-fraction)
          (homographic-function))

  (begin

    (define (run-demonstration)

      (define cf+1/2 (make-homographic-function-operator 2 1 0 2))
      (define cf/2 (make-homographic-function-operator 1 0 0 2))
      (define cf/4 (make-homographic-function-operator 1 0 0 4))
      (define 1/cf (make-homographic-function-operator 0