Continued fraction/Arithmetic/G(matrix ng, continued fraction n): Difference between revisions
Content added Content deleted
| Line 2,971: | Line 2,971: | ||
{{trans|ATS}} |
{{trans|ATS}} |
||
{{trans|C}} |
{{trans|C}} |
||
'''WARNING: THIS IS BROKEN. I am fixing it right now.''' |
|||
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. |
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. |
||
Revision as of 01:49, 24 February 2023
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
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
----------------------------------------------------------------------
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 = "+"
val dot_operator = "⋅"
val centered_ellipsis = "⋯"
val right_arrow = "→"
(*------------------------------------------------------------------*)
(* 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:
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:
C
With a garbage collector
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
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
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
Fortran
WARNING: THIS IS BROKEN. I am fixing it right now.
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.)
!---------------------------------------------------------------------
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
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
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
subroutine cf_memo_t_finalize (memo)
type(cf_memo_t), intent(inout) :: memo
deallocate (memo%storage)
end subroutine cf_memo_t_finalize
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
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
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
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
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
subroutine cf_update (cf, needed)
class(cf_t), intent(inout) :: cf
integer, intent(in) :: needed
integer, pointer :: storage1(:)
integer :: i
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) = cf%memo%storage(0:cf%m)
deallocate (cf%memo%storage)
cf%memo%storage => storage1
call cf_get_more_terms (cf, needed)
end if
end subroutine cf_update
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
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
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
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
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
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
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
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
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
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
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)
end block
term_exists = .true.
term = q
done = .true.
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
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
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)
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)
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)
end program univariate_continued_fraction_task
!---------------------------------------------------------------------
- Output:
$ gfortran -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,...]
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, ...]
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,&1to 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
Julia
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
Nim
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
Phix
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
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
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)
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
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,...]
Tcl
This uses the Generator class, R2CF class and printcf procedure from the r2cf task.
# The single-operand version of the NG operator, using our little generator framework
oo::class create NG1 {
superclass Generator
variable a1 a b1 b cf
constructor args {
next
lassign $args a1 a b1 b
}
method Ingress n {
lassign [list [expr {$a + $a1*$n}] $a1 [expr {$b + $b1*$n}] $b1] \
a1 a b1 b
}
method NeedTerm? {} {
expr {$b1 == 0 || $b == 0 || $a/$b != $a1/$b1}
}
method Egress {} {
set n [expr {$a/$b}]
lassign [list $b1 $b [expr {$a1 - $b1*$n}] [expr {$a - $b*$n}]] \
a1 a b1 b
return $n
}
method EgressDone {} {
if {[my NeedTerm?]} {
set a $a1
set b $b1
}
tailcall my Egress
}
method Done? {} {
expr {$b1 == 0 && $b == 0}
}
method operand {N} {
set cf $N
return [self]
}
method Produce {} {
while 1 {
set n [$cf]
if {![my NeedTerm?]} {
yield [my Egress]
}
my Ingress $n
}
while {![my Done?]} {
yield [my EgressDone]
}
}
}
Demonstrating:
# The square root of 2 as a continued fraction in the framework
oo::class create Root2 {
superclass Generator
method apply {} {
yield 1
while {[self] ne ""} {
yield 2
}
}
}
set op [[NG1 new 2 1 0 2] operand [R2CF new 13/11]]
printcf "\[1;5,2\] + 1/2" $op
set op [[NG1 new 7 0 0 22] operand [R2CF new 22/7]]
printcf "\[3;7\] * 7/22" $op
set op [[NG1 new 2 1 0 2] operand [R2CF new 22/7]]
printcf "\[3;7\] + 1/2" $op
set op [[NG1 new 1 0 0 4] operand [R2CF new 22/7]]
printcf "\[3;7\] / 4" $op
set op [[NG1 new 0 1 1 0] operand [Root2 new]]
printcf "1/\u221a2" $op 20
set op [[NG1 new 1 1 0 2] operand [Root2 new]]
printcf "(1+\u221a2)/2" $op 20
printcf "approx val" [R2CF new 24142136 20000000]
- Output:
[1;5,2] + 1/2 -> 1,1,2,7 [3;7] * 7/22 -> 1 [3;7] + 1/2 -> 3,1,1,1,4 [3;7] / 4 -> 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,… approx val -> 1,4,1,4,1,4,1,4,1,4,3,2,1,9,5
Wren
import "/dynamic" for Tuple
var CFData = Tuple.create("Tuple", ["str", "ng", "r", "gen"])
var r2cf = Fn.new { |frac|
var num = frac[0]
var den = frac[1]
while (den.abs != 0) {
var div = (num/den).truncate
var rem = num % den
num = den
den = rem
Fiber.yield(div)
}
}
var d2cf = Fn.new { |d|
while (true) {
var div = d.floor
var rem = d - div
Fiber.yield(div)
if (rem == 0) break
d = 1 / rem
}
}
var root2 = Fn.new {
Fiber.yield(1)
while (true) Fiber.yield(2)
}
var recipRoot2 = Fn.new {
Fiber.yield(0)
Fiber.yield(1)
while (true) Fiber.yield(2)
}
class NG {
construct new(a1, a, b1, b) {
_a1 = a1
_a = a
_b1 = b1
_b = b
}
ingress(n) {
var t = _a
_a = _a1
_a1 = t + _a1 * n
t = _b
_b = _b1
_b1 = t + _b1 * n
}
egress() {
var n = (_a/_b).truncate
var t = _a
_a = _b
_b = t - _b * n
t = _a1
_a1 = _b1
_b1 = t - _b1 * n
return n
}
needTerm { (_b == 0 || _b1 == 0) || ((_a / _b) != (_a1 / _b1)) }
egressDone {
if (needTerm) {
_a = _a1
_b = _b1
}
return egress()
}
done { _b == 0 && _b1 == 0 }
}
var data = [
CFData.new("[1;5,2] + 1/2 ", [2, 1, 0, 2], [13, 11], r2cf),
CFData.new("[3;7] + 1/2 ", [2, 1, 0, 2], [22, 7], r2cf),
CFData.new("[3;7] divided by 4 ", [1, 0, 0, 4], [22, 7], r2cf),
CFData.new("sqrt(2) ", [0, 1, 1, 0], [ 0, 0], recipRoot2),
CFData.new("1 / sqrt(2) ", [0, 1, 1, 0], [ 0, 0], root2),
CFData.new("(1 + sqrt(2)) / 2 ", [1, 1, 0, 2], [ 0, 0], root2),
CFData.new("(1 + 1 / sqrt(2)) / 2", [1, 1, 0, 2], [ 0, 0], recipRoot2)
]
System.print("Produced by NG class:")
for (cfd in data) {
System.write("%(cfd.str) -> ")
var a1 = cfd.ng[0]
var a = cfd.ng[1]
var b1 = cfd.ng[2]
var b = cfd.ng[3]
var op = NG.new(a1, a, b1, b)
var seq = []
var i = 0
var fib = Fiber.new(cfd.gen)
while (i < 20) {
var j = fib.call(cfd.r)
if (j) seq.add(j) else break
i = i + 1
}
for (n in seq) {
if (!op.needTerm) System.write(" %(op.egress()) ")
op.ingress(n)
}
while (true) {
System.write(" %(op.egressDone) ")
if (op.done) break
}
System.print()
}
System.print("\nProduced by direct calculation:")
var data2 = [
["(1 + sqrt(2)) / 2 ", (1 + 2.sqrt) / 2],
["(1 + 1 / sqrt(2)) / 2", (1 + 1 / 2.sqrt) / 2]
]
for (p in data2) {
var seq = []
var fib = Fiber.new(d2cf)
var i = 0
while (i < 20) {
var j = fib.call(p[1])
if (j) seq.add(j) else break
i = i + 1
}
System.print("%(p[0]) -> %(seq.join(" "))")
}- 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
zkl
class NG{
fcn init(_a1,_a, _b1,_b){ var a1=_a1,a=_a, b1=_b1,b=_b; }
var [proxy] done =fcn{ b==0 and b1==0 };
var [proxy] needterm=fcn{ (b==0 or b1==0) or (a/b!=a1/b1) };
fcn ingress(n){
t:=self.copy(True); // tmp copy of vars for eager vs late evaluation
a,a1=t.a1, t.a + n*t.a1;
b,b1=t.b1, t.b + n*t.b1;
}
fcn egress{
n,t:=a/b,self.copy(True);
a,b =t.b, t.a - n*t.b;
a1,b1=t.b1,t.a1 - n*t.b1;
n
}
fcn egress_done{
if(needterm) a,b=a1,b1;
egress()
}
} // from task: Continued fraction/Arithmetic/Construct from rational number
fcn r2cf(nom,dnom){ // -->Walker (iterator)
Walker.tweak(fcn(_,state){
nom,dnom:=state;
if(dnom==0) return(Void.Stop);
n,d:=nom.divr(dnom);
state.clear(dnom,d);
n
}.fp1(List(nom,dnom)))
}data:=T(T("[1;5,2] + 1/2", T(2,1,0,2), T(13,11)),
T("[3;7] + 1/2", T(2,1,0,2), T(22, 7)),
T("[3;7] divided by 4", T(1,0,0,4), T(22, 7)));
foreach string,ng,r in (data){
print("%-20s-->".fmt(string));
op:=NG(ng.xplode());
foreach n in (r2cf(r.xplode())){
if(not op.needterm) print(" %s".fmt(op.egress()));
op.ingress(n);
}
do{ print(" ",op.egress_done()) }while(not op.done);
println();
}- 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