How to resolve the algorithm Continued fraction/Arithmetic/G(matrix ng, continued fraction n) step by step in the ATS programming language
How to resolve the algorithm Continued fraction/Arithmetic/G(matrix ng, continued fraction n) step by step in the ATS programming language
Table of Contents
Problem Statement
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: I output a term if the integer parts of
a b
{\displaystyle {\frac {a}{b}}}
and
a
1
b
1
{\displaystyle {\frac {a_{1}}{b_{1}}}}
are equal. Otherwise I input a term from N. If I need a term from N but N has no more terms I inject
∞
{\displaystyle \infty }
. When I input a term t my internal state:
[
a
1
a
b
1
b
]
{\displaystyle {\begin{bmatrix}a_{1}&a\b_{1}&b\end{bmatrix}}}
is transposed thus
[
a +
a
1
∗ t
a
1
b +
b
1
∗ t
b
1
]
{\displaystyle {\begin{bmatrix}a+a_{1}*t&a_{1}\b+b_{1}*t&b_{1}\end{bmatrix}}}
When I output a term t my internal state:
[
a
1
a
b
1
b
]
{\displaystyle {\begin{bmatrix}a_{1}&a\b_{1}&b\end{bmatrix}}}
is transposed thus
[
b
1
b
a
1
−
b
1
∗ t
a − b ∗ t
]
{\displaystyle {\begin{bmatrix}b_{1}&b\a_{1}-b_{1}t&a-bt\end{bmatrix}}}
When I need a term t but there are no more my internal state:
[
a
1
a
b
1
b
]
{\displaystyle {\begin{bmatrix}a_{1}&a\b_{1}&b\end{bmatrix}}}
is transposed thus
[
a
1
a
1
b
1
b
1
]
{\displaystyle {\begin{bmatrix}a_{1}&a_{1}\b_{1}&b_{1}\end{bmatrix}}}
I am done when b1 and b are zero. Demonstrate your solution by calculating: Using a generator for
2
{\displaystyle {\sqrt {2}}}
(e.g., from Continued fraction) calculate
1
2
{\displaystyle {\frac {1}{\sqrt {2}}}}
. 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
1 +
1
2
2
{\displaystyle {\frac {1+{\frac {1}{\sqrt {2}}}}{2}}}
do this now to cross the starting line and begin the race.
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Continued fraction/Arithmetic/G(matrix ng, continued fraction n) step by step in the ATS programming language
Source code in the ats programming language
(*------------------------------------------------------------------*)
#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) =
() - 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}
..
(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 (cf, needed)
else
let (* Provides twice the room needed. *)
val n1 = needed + needed
val memo1 = arrayref_make_elt (n1, g0i2i 0)
val () =
let
var i : [i : nat] size_t i
in
for (i := i2sz 0; i < (cf.m); i := succ i)
memo1[i] := memo[i]
end
val cf1 =
'{
terminated = false,
m = (cf.m),
n = n1,
memo = memo1,
gen = (cf.gen)
}
in
_cf_get_more_terms (cf1, needed)
end
end
in (* local *)
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 (cf, succ i)
in
cf := cf1;
if i < (cf1.m) then
Some (arrayref_get_at (cf1.memo, i))
else
None ()
end
implement {tk} {tki}
cf_get_at_gint (cf, i) =
cf_get_at_guint (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, ", ",
centered_ellipsis, "] "))
else
begin
case+ cf[i] of
| None () =>
strptr2string (string_append (accum, "] "))
| Some term =>
let
val term_str = tostring_val term
val accum =
strptr2string (string_append (accum, sep, "",
term_str, ""))
val sep =
if sep = "[ " then
"; "
else
", "
in
loop (succ i, cf, sep, accum)
end
end
in
loop (i2sz 0, cf, "[ ", "")
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 () =
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 (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 cfref
val gen2 = hfunc_apply_generator2generator (f, gen1)
in
cfref_make 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 n
val d = ref d
fn
gen () : 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 () =
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
macdef cfref2ml = cfref2mathml
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
macdef apply = hfunc_apply
macdef text (s) =
strptr2string (string_append ("", ,(s), "
"))
macdef becomes =
strptr2string (string_append ("", right_arrow, " "))
macdef start_math =
"
macdef stop_math = ""
macdef start_table = ""
macdef stop_table = ""
macdef left_side (s) =
strptr2string
(string_append
("", ,(s), " "))
macdef middle (s) =
strptr2string
(string_append
("", ,(s), " "))
macdef right_side (s) =
strptr2string
(string_append
("", ,(s), " "))
macdef entry (left, right) =
strptr2string
(string_append
("",
left_side (,(left)),
middle becomes,
right_side (,(right)),
""))
macdef num s =
strptr2string (string_append ("", ,(s), " "))
macdef id s =
strptr2string (string_append ("", ,(s), " "))
macdef oper s =
strptr2string (string_append ("", ,(s), " "))
macdef frac (n, d) =
strptr2string (string_append ("", ,(n), ,(d),
""))
macdef numfrac (n, d) = frac (num ,(n), num ,(d))
macdef sqrt_of (s) =
strptr2string (string_append ("", ,(s), " "))
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
(*------------------------------------------------------------------*)
(*------------------------------------------------------------------*)
(* 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 = 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 (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 (workspace, 0);
{prval () = opt_some term};
arrayptr_set_at (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 (workspace, 2)
and b = arrayptr_get_at (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
(*------------------------------------------------------------------*)
(*------------------------------------------------------------------*)
#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"
(* For convenience (because the prelude provides it), we will use
integer division with truncation towards zero. *)
infixl ( / ) div
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod
(*------------------------------------------------------------------*)
(* The definition of a continued fraction, and a few simple ones. *)
typedef cf (tk : tkind) = stream (g0int tk)
(* A "continued fraction" with no terms. *)
fn {tk : tkind}
cfnil ()
: cf tk =
stream_make_nil ()
(* A continued fraction of one term followed by more terms. *)
fn {tk : tkind}
cfcons (term : g0int tk,
more : cf tk)
: cf tk =
stream_make_cons (term, more)
(* A continued fraction with all terms equal. *)
fn {tk : tkind}
repeat_forever (term : g0int tk)
: cf tk =
let
fun recurs () : stream_con (g0int tk) =
stream_cons (term, $delay recurs ())
in
$delay recurs ()
end
(* The square root of two. *)
fn {tk : tkind}
sqrt2 ()
: cf tk =
cfcons (g0i2i 1, repeat_forever (g0i2i 2))
(*------------------------------------------------------------------*)
(* A continued fraction for a rational number. *)
typedef ratnum (tk : tkind) = @(g0int tk, g0int tk)
fn {tk : tkind}
r2cf_integers (n : g0int tk,
d : g0int tk)
: cf tk =
let
fun recurs (n : g0int tk,
d : g0int tk)
: cf tk =
if iseqz d then
cfnil ()
else
cfcons (n div d, recurs (d, n rem d))
in
recurs (n, d)
end
fn {tk : tkind}
r2cf_ratnum (r : ratnum tk)
: cf tk =
r2cf_integers (r.0, r.1)
overload r2cf with r2cf_integers
overload r2cf with r2cf_ratnum
(*------------------------------------------------------------------*)
(* Application of a homographic function to a continued fraction. *)
typedef ng4 (tk : tkind) = @(g0int tk, g0int tk,
g0int tk, g0int tk)
fn {tk : tkind}
apply_ng4 (ng4 : ng4 tk,
other_cf : cf tk)
: cf tk =
let
typedef t = g0int tk
fun
recurs (a1 : t,
a : t,
b1 : t,
b : t,
other_cf : cf tk)
: stream_con t =
let
fn {}
eject_term (a1 : t,
a : t,
b1 : t,
b : t,
other_cf : cf tk,
term : t)
: stream_con t =
stream_cons (term, $delay recurs (b1, b, a1 - (b1 * term),
a - (b * term), other_cf))
fn {}
absorb_term (a1 : t,
a : t,
b1 : t,
b : t,
other_cf : cf tk)
: stream_con t =
case+ !other_cf of
| stream_nil () =>
recurs (a1, a1, b1, b1, other_cf)
| stream_cons (term, rest) =>
recurs (a + (a1 * term), a1, b + (b1 * term), b1, rest)
in
if iseqz b1 && iseqz b then
stream_nil ()
else if iseqz b1 || iseqz b then
absorb_term (a1, a, b1, b, other_cf)
else
let
val q1 = a1 div b1
and q = a div b
in
if q1 = q then
eject_term (a1, a, b1, b, other_cf, q)
else
absorb_term (a1, a, b1, b, other_cf)
end
end
val @(a1, a, b1, b) = ng4
in
$delay recurs (a1, a, b1, b, other_cf)
end
(*------------------------------------------------------------------*)
(* Some special cases of homographic functions. *)
fn {tk : tkind}
integer_add_cf (n : g0int tk,
cf : cf tk)
: cf tk =
apply_ng4 (@(g0i2i 1, n, g0i2i 0, g0i2i 1), cf)
fn {tk : tkind}
cf_add_ratnum (cf : cf tk,
r : ratnum tk)
: cf tk =
let
val @(n, d) = r
in
apply_ng4 (@(d, n, g0i2i 0, d), cf)
end
fn {tk : tkind}
cf_mul_ratnum (cf : cf tk,
r : ratnum tk)
: cf tk =
let
val @(n, d) = r
in
apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), cf)
end
fn {tk : tkind}
cf_div_integer (cf : cf tk,
n : g0int tk)
: cf tk =
apply_ng4 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i n), cf)
fn {tk : tkind}
integer_div_cf (n : g0int tk,
cf : cf tk)
: cf tk =
apply_ng4 (@(g0i2i 0, g0i2i n, g0i2i 1, g0i2i 0), cf)
overload + with integer_add_cf
overload + with cf_add_ratnum
overload * with cf_mul_ratnum
overload / with cf_div_integer
overload / with integer_div_cf
(*------------------------------------------------------------------*)
(* cf2string: convert a continued fraction to a string. *)
fn {tk : tkind}
cf2string_max_terms_given
(cf : cf tk,
max_terms : intGte 1)
: string =
let
fun
loop (i : intGte 0,
cf : cf tk,
accum : List0 string)
: List0 string =
case+ !cf of
| stream_nil () => list_cons ("]", accum)
| stream_cons (term, rest) =>
if i = max_terms then
list_cons (",...]", accum)
else
let
val accum =
list_cons
(tostring_val term,
(case+ i of
| 0 => accum
| 1 => list_cons (";", accum)
| _ => list_cons (",", accum)) : List0 string)
in
loop (succ i, rest, accum)
end
val string_lst = list_vt2t (reverse (loop (0, cf, list_sing "[")))
in
strptr2string (stringlst_concat string_lst)
end
extern fn {tk : tkind}
cf2string$max_terms :
() -> intGte 1
implement {tk} cf2string$max_terms () = 20
fn {tk : tkind}
cf2string_max_terms_default
(cf : cf tk)
: string =
cf2string_max_terms_given (cf, cf2string$max_terms ())
overload cf2string with cf2string_max_terms_given
overload cf2string with cf2string_max_terms_default
(*------------------------------------------------------------------*)
fn {tk : tkind}
show (expression : string,
cf : cf tk)
: void =
begin
print! expression;
print! " => ";
println! (cf2string cf);
end
implement
main () =
let
val cf_13_11 = r2cf (13, 11)
val cf_22_7 = r2cf (22, 7)
val cf_sqrt2 = sqrt2 ()
val cf_1_sqrt2 = 1 / cf_sqrt2
in
show ("13/11", cf_13_11);
show ("22/7", cf_22_7);
show ("sqrt(2)", cf_sqrt2);
show ("13/11 + 1/2", cf_13_11 + @(1, 2));
show ("22/7 + 1/2", cf_22_7 + @(1, 2));
show ("(22/7)/4", cf_22_7 * @(1, 4));
show ("1/sqrt(2)", cf_1_sqrt2);
show ("(2 + sqrt(2))/4", apply_ng4 (@(1, 2, 0, 4), cf_sqrt2));
(* To show it can be done, write the following without using
results already obtained: *)
show ("(1 + 1/sqrt(2))/2", (1 + 1/sqrt2())/2);
0
end
(*------------------------------------------------------------------*)
You may also check:How to resolve the algorithm Program termination step by step in the Simula programming language
You may also check:How to resolve the algorithm Verify distribution uniformity/Naive step by step in the AutoHotkey programming language
You may also check:How to resolve the algorithm Quine step by step in the zkl programming language
You may also check:How to resolve the algorithm Constrained random points on a circle step by step in the Julia programming language
You may also check:How to resolve the algorithm Ethiopian multiplication step by step in the Erlang programming language