Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2)
Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2) is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
This task performs the basic mathematical functions on 2 continued fractions. This requires the full version of matrix NG:
I may perform perform the following operations:
- Input the next term of continued fraction N1
- Input the next term of continued fraction N2
- Output a term of the continued fraction resulting from the operation.
I output a term if the integer parts of and and and are equal. Otherwise I input a term from continued fraction N1 or continued fraction N2. If I need a term from N but N has no more terms I inject .
When I input a term t from continued fraction N1 I change my internal state:
- is transposed thus
When I need a term from exhausted continued fraction N1 I change my internal state:
- is transposed thus
When I input a term t from continued fraction N2 I change my internal state:
- is transposed thus
When I need a term from exhausted continued fraction N2 I change my internal state:
- is transposed thus
When I output a term t I change my internal state:
- is transposed thus
When I need to choose to input from N1 or N2 I act:
- if b and b2 are zero I choose N1
- if b is zero I choose N2
- if b2 is zero I choose N2
- if abs( is greater than abs( I choose N1
- otherwise I choose N2
When performing arithmetic operation on two potentially infinite continued fractions it is possible to generate a rational number. eg * should produce 2. This will require either that I determine that my internal state is approaching infinity, or limiting the number of terms I am willing to input without producing any output.
ATS
Using long long int as the integer type
The following code uses GNU C extensions: integer operations that detect overflow. All the calculations (except when converting a result to floating point) are done using ordinary machine integers.
When overflow occurs—for instance when one is calculating sqrt(2)*sqrt(2)—a large term appears where probably there should have been an infinity. Heuristics might be used to elide the erroneous term, but I have not implemented such heuristics here.
(*------------------------------------------------------------------*)
#include "share/atspre_staload.hats"
%{#
#include <limits.h>
#include <float.h>
#include <math.h>
%}
(* Fill in something that is missing from the prelude. (See
https://sourceforge.net/p/chemoelectric/ats2-xprelude for a way to
fill in many gaps in the prelude.) *)
implement g0int2float<llintknd,dblknd> x = $UNSAFE.cast x
#define NIL list_nil ()
#define :: list_cons
(* Truncate quotients towards zero, simply because this is all that
the prelude has functions for. *)
infixl ( / ) div divrem
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod
macdef divrem (a, b) =
let
val x = ,(a)
and y = ,(b)
in
(* Optimizing C compilers will compute both quotient and remainder
at the same time. *)
@(x \g0int_div y, x \g0int_mod y)
end
(*------------------------------------------------------------------*)
(* We have to catch integer overflows! For this, we will use GNU
extensions to the C language. *)
%{#
#define add_overflow_llint(x, y, pz) \
(__builtin_saddll_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define sub_overflow_llint(x, y, pz) \
(__builtin_ssubll_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define mul_overflow_llint(x, y, pz) \
(__builtin_smulll_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
%}
exception gint_overflow of ()
local
extern fn
add_overflow_llint :
(llint, llint, &llint? >> llint) -< !wrt > bool
= "mac#add_overflow_llint"
extern fn
sub_overflow_llint :
(llint, llint, &llint? >> llint) -< !wrt > bool
= "mac#sub_overflow_llint"
extern fn
mul_overflow_llint :
(llint, llint, &llint? >> llint) -< !wrt > bool
= "mac#mul_overflow_llint"
in (* local *)
fn
g0int_add_overflow_exn_llint (x : llint, y : llint) :<!exn> llint =
let
var z : llint?
val overflow = $effmask_wrt add_overflow_llint (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
fn
g0int_sub_overflow_exn_llint (x : llint, y : llint) :<!exn> llint =
let
var z : llint?
val overflow = $effmask_wrt sub_overflow_llint (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
fn
g0int_mul_overflow_exn_llint (x : llint, y : llint) :<!exn> llint =
let
var z : llint?
val overflow = $effmask_wrt mul_overflow_llint (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
end (* local *)
extern fn {tk : tkind}
g0int_add_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_sub_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_mul_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
implement
g0int_add_overflow_exn<llintknd> (x, y) =
g0int_add_overflow_exn_llint (x, y)
implement
g0int_sub_overflow_exn<llintknd> (x, y) =
g0int_sub_overflow_exn_llint (x, y)
implement
g0int_mul_overflow_exn<llintknd> (x, y) =
g0int_mul_overflow_exn_llint (x, y)
infixl ( + ) +!
infixl ( - ) -!
infixl ( * ) *!
overload +! with g0int_add_overflow_exn
overload -! with g0int_sub_overflow_exn
overload *! with g0int_mul_overflow_exn
(*------------------------------------------------------------------*)
(* Continued fractions.
cf_generator tk -- A closure that produces terms of type g0int tk,
sequentially.
cf tk -- A structure from which one can get the ith
term of a continued fraction. It gets terms
from a cf_generator tk. *)
typedef integer = llint
(* I am going to let the most negative value of an integer represent
"no term", or in other words "negative infinity". Safer would be to
use "Option(integer)", but I want to play with the more efficient
representation. *)
typedef cf_generator =
() -<cloref1> integer
extern fn {} neginf : () -<> integer
implement {} neginf () = $extval (integer, "LLONG_MIN")
local
typedef _cf (terminated : bool,
m : int,
n : int) =
[m <= n]
@{
terminated = bool terminated, (* No more terms? *)
m = size_t m, (* The number of terms computed so far. *)
n = size_t n, (* The size of memo storage.*)
memo = arrayref (integer, n), (* Memoized terms. *)
gen = cf_generator (* A thunk to generate terms. *)
}
typedef _cf (m : int) =
[terminated : bool]
[n : int | m <= n]
_cf (terminated, m, n)
typedef _cf =
[m : int]
_cf m
fn
_cf_get_more_terms
{terminated : bool}
{m : int}
{n : int}
{needed : int | m <= needed; needed <= n}
(cf : _cf (terminated, m, n),
needed : size_t needed)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)
macdef memo = cf.memo
fun
loop {i : int | m <= i; i <= needed}
.<needed - i>.
(i : size_t i)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
if i = needed then
@{
terminated = false,
m = needed,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
else
let
val term = (cf.gen) ()
in
if term <> neginf () then
begin
memo[i] := term;
loop (succ i)
end
else
@{
terminated = true,
m = i,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
end
in
loop (cf.m)
end
fn
_cf_update {terminated : bool}
{m : int}
{n : int | m <= n}
{needed : int}
(cf : _cf (terminated, m, n),
needed : size_t needed)
: [terminated1 : bool]
[m1 : int | m <= m1]
[n1 : int | m1 <= n1]
_cf (terminated1, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)
macdef memo = cf.memo
macdef gen = cf.gen
in
if (cf.terminated) then
cf
else if needed <= (cf.m) then
cf
else if needed <= (cf.n) then
_cf_get_more_terms (cf, needed)
else
let (* Provides twice the room needed. *)
val n1 = needed + needed
val memo1 = arrayref_make_elt (n1, g0i2i 0)
val () =
let
var i : [i : nat] size_t i
in
for (i := i2sz 0; i < (cf.m); i := succ i)
memo1[i] := memo[i]
end
val cf1 =
@{
terminated = false,
m = (cf.m),
n = n1,
memo = memo1,
gen = (cf.gen)
}
in
_cf_get_more_terms (cf1, needed)
end
end
in (* local *)
typedef cf = ref _cf
extern fn
cf_make :
cf_generator -> cf
extern fn
cf_get_at_size :
{i : int}
(cf, size_t i) -> integer
extern fn
cf_get_at_int :
{i : nat}
(cf, int i) -> integer
(* The precedence of the overloads has to exceed that of ref[] *)
overload cf_get_at with cf_get_at_size of 1
overload cf_get_at with cf_get_at_int of 1
overload [] with cf_get_at of 1
implement
cf_make gen =
let
#ifndef CF_START_SIZE #then
#define CF_START_SIZE 32
#endif
in
ref
@{
terminated = false,
m = i2sz 0,
n = i2sz CF_START_SIZE,
memo = arrayref_make_elt (i2sz CF_START_SIZE, g0i2i 0),
gen = gen
}
end
implement
cf_get_at_size (cf, i) =
let
prval () = lemma_g1uint_param i
val cf1 = _cf_update (!cf, succ i)
in
!cf := cf1;
if i < (cf1.m) then
arrayref_get_at<integer> (cf1.memo, i)
else
neginf ()
end
implement
cf_get_at_int (cf, i) =
cf_get_at_size (cf, g1i2u i)
end (* local *)
(*------------------------------------------------------------------*)
(* Converting a continued fraction to a string. *)
extern fn cf2string_with_default_max_terms : cf -> string
extern fn cf2string_given_max_terms : (cf, Size_t) -> string
overload cf2string with cf2string_with_default_max_terms
overload cf2string with cf2string_given_max_terms
val cf2string_default_max_terms : ref Size_t = ref (i2sz 20)
implement
cf2string_with_default_max_terms cf =
cf2string_given_max_terms (cf, !cf2string_default_max_terms)
implement
cf2string_given_max_terms (cf, max_terms) =
let
fun
loop (i : Size_t,
accum : List0 string)
: List0 string =
let
val term = cf[i]
in
if i = max_terms then
begin
if term = neginf () then
"]" :: accum
else
",...]" :: accum
end
else if term = neginf () then
"]" :: accum
else
let
val separator =
(if i = i2sz 0 then
""
else if i = i2sz 1 then
";"
else
",")
and term_str = tostring_val<integer> term
in
loop (succ i, term_str :: separator :: accum)
end
end
val string_lst = loop (i2sz 0, "[" :: NIL)
in
strptr2string (stringlst_concat (list_vt2t (reverse string_lst)))
end
(*------------------------------------------------------------------*)
(* Compute the ith convergent as a floating-point number. Convergents
are indexed starting at 0. *)
extern fn
cf_convergent_double_size :
{i : int}
(cf, size_t i) -> double
extern fn
cf_convergent_double_int :
{i : nat}
(cf, int i) -> double
overload cf_convergent_double with cf_convergent_double_size
overload cf_convergent_double with cf_convergent_double_int
implement
cf_convergent_double_size {i} (cf, i) =
let
prval () = lemma_g1uint_param i
fun
loop {j : nat | j <= i + 1}
.<(i + 1) - j>.
(j : size_t j,
accum : double)
: double =
if j = succ i then
accum
else
let
val term = cf[i - j]
val x =
(if term = neginf () then
$extval (double, "INFINITY")
else
g0i2f term) : double
in
loop (succ j, x + (1.0 / accum))
end
in
loop (i2sz 0, 0.0)
end
implement
cf_convergent_double_int (cf, i) =
cf_convergent_double_size (cf, g1i2u i)
(*------------------------------------------------------------------*)
(* A continued fraction for the square root of two. *)
val sqrt2 : cf =
cf_make
let
val first : ref bool = ref true
in
lam () =<cloref1>
let
val fst = !first
in
!first := false;
(if fst then g0i2i 1 else g0i2i 2) : integer
end
end
(*------------------------------------------------------------------*)
(* The continued fraction for a rational number. *)
typedef ratnum = @(integer, integer)
fn
r2cf (ratnum : ratnum) : cf =
cf_make
let
val ratnum_ref : ref ratnum = ref ratnum
in
lam () =<cloref1>
let
val @(n, d) = !ratnum_ref
in
if iseqz d then
neginf ()
else
let
val @(q, r) = n divrem d
in
!ratnum_ref := @(d, r);
q
end
end
end
(*------------------------------------------------------------------*)
(* Application of a homographic function to a continued fraction. *)
typedef ng4 = @(integer, integer, integer, integer)
fn
apply_ng4 (ng4 : ng4, x : cf) : cf =
cf_make
let
val state : ref @(ng4, Size_t) = ref @(ng4, i2sz 0)
in
lam () =<cloref1>
let
fnx
loop (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
if (iseqz b1) * (iseqz b) then
neginf ()
else if (iseqz b1) + (iseqz b) then
absorb_term (a1, a, b1, b, i)
else
let
val q1 = a1 div b1
and q = a div b
in
if q1 <> q then
absorb_term (a1, a, b1, b, i)
else
begin
!state :=
@(@(b1, b, a1 - (b1 * q), a - (b * q)), i);
q
end
end
and
absorb_term (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
let
val term = x[i]
in
if term <> neginf () then
loop (a + (a1 * term), a1,
b + (b1 * term), b1, succ i)
else
loop (a1, a1, b1, b1, succ i)
end
val @(@(a1, a, b1, b), i) = !state
in
loop (a1, a, b1, b, i)
end
end
(*------------------------------------------------------------------*)
(* Some basic operations involving only one continued fraction. *)
fn
cf_neg (x : cf) : cf =
apply_ng4 (@(g0i2i ~1, g0i2i 0, g0i2i 0, g0i2i 1), x)
fn
cf_add_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(d, n, g0i2i 0, d), x)
end
fn
ratnum_add_cf (x : ratnum, y : cf) : cf =
cf_add_ratnum (y, x)
fn
cf_add_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(y, g0i2i 1))
fn
integer_add_cf (x : integer, y : cf) : cf =
cf_add_ratnum (y, @(x, g0i2i 1))
fn
cf_sub_ratnum (x : cf, y : ratnum) : cf =
cf_add_ratnum (x, @(~(y.0), (y.1)))
fn
ratnum_sub_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(~d, n, g0i2i 0, d), y)
end
fn
cf_sub_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(~y, g0i2i 1))
fn
integer_sub_cf (x : integer, y : cf) : cf =
ratnum_sub_cf (@(x, g0i2i 1), y)
fn
cf_mul_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), x)
end
fn
ratnum_mul_cf (x : ratnum, y : cf) : cf =
cf_mul_ratnum (y, x)
fn
cf_mul_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(y, g0i2i 1))
fn
integer_mul_cf (x : integer, y : cf) : cf =
cf_mul_ratnum (y, @(x, g0i2i 1))
fn
cf_div_ratnum (x : cf, y : ratnum) : cf =
cf_mul_ratnum (x, @(y.1, y.0))
fn
ratnum_div_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(g0i2i 0, n, d, g0i2i 0), y)
end
fn
cf_div_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(g0i2i 1, y))
fn
integer_div_cf (x : integer, y : cf) : cf =
ratnum_div_cf (@(x, g0i2i 1), y)
overload ~ with cf_neg
overload + with cf_add_ratnum
overload + with ratnum_add_cf
overload + with cf_add_integer
overload + with integer_add_cf
overload - with cf_sub_ratnum
overload - with ratnum_sub_cf
overload - with cf_sub_integer
overload - with integer_sub_cf
overload * with cf_mul_ratnum
overload * with ratnum_mul_cf
overload * with cf_mul_integer
overload * with integer_mul_cf
overload / with cf_div_ratnum
overload / with ratnum_div_cf
overload / with cf_div_integer
overload / with integer_div_cf
(*------------------------------------------------------------------*)
(* Application of a bihomographic function to a continued fraction. *)
typedef ng8 = @(integer, integer, integer, integer,
integer, integer, integer, integer)
fn
apply_ng8 (ng8 : ng8, x : cf, y : cf) : cf =
cf_make
let
typedef state = '(ng8, Size_t, Size_t, bool)
val state_ref : ref state = ref '(ng8, i2sz 0, i2sz 0, false)
in
lam () =<cloref1>
let
fnx
loop (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
in
if iseqz b12 * iseqz b1 * iseqz b2 * iseqz b then
neginf ()
else if iseqz b * iseqz b2 then
absorb_x_term state
else if iseqz b + iseqz b2 then
absorb_y_term state
else if iseqz b1 then
absorb_x_term state
else
let
val @(q, r) = a divrem b // a/b = q + r/b
and @(q1, r1) = a1 divrem b1 // a1/b1 = q1 + r1/b1
and @(q2, r2) = a2 divrem b2 // a2/b2 = q2 + r2/b2
val q1_diff = abs (q1 -! q)
and q2_diff = abs (q2 -! q)
in
if q1_diff > q2_diff then
absorb_x_term state
else if q1_diff < q2_diff then
absorb_y_term state
else if iseqz b12 || q <> a12 div b12 then
let
(* Because q = q1 = q2, the following also are
true:
a1/b1 - a/b = r1/b1 - r/b
a2/b2 - a/b = r2/b2 - r/b
But now the numerators are smaller, and so
overflow is less likely to occur when we put
terms over a common denominator. *)
(* Put r, r1, r2 over a common denominator,
b*b1*b2. *)
val n = r *! b1 *! b2
and n1 = r1 *! b *! b2
and n2 = r2 *! b *! b1
in
if abs (n1 -! n) > abs (n2 -! n) then
absorb_x_term state
else
absorb_y_term state
end
else
begin
!state_ref :=
'(@(b12, b1, b2, b,
a12 -! (b12 *! q), a1 -! (b1 *! q),
a2 -! (b2 *! q), a -! (b *! q)),
i, j, overflow);
q
end
end
end
and
absorb_x_term (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
val term = (if overflow then neginf () else x[i])
and i = succ i
in
if term <> neginf () then
begin
try
loop '(@(a2 +! (a12 *! term),
a +! (a1 *! term), a12, a1,
b2 +! (b12 *! term),
b +! (b1 *! term), b12, b1),
i, j, overflow)
with
| ~ gint_overflow () =>
loop '(@(a12, a1, a12, a1, b12, b1, b12, b1),
pred i, j, true)
end
else
loop '(@(a12, a1, a12, a1, b12, b1, b12, b1),
i, j, overflow)
end
and
absorb_y_term (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
val term = (if overflow then neginf () else y[j])
and j = succ j
in
if term <> neginf () then
begin
try
loop '(@(a1 +! (a12 *! term), a12,
a +! (a2 *! term), a2,
b1 +! (b12 *! term), b12,
b +! (b2 *! term), b2),
i, j, overflow)
with
| ~ gint_overflow () =>
loop '(@(a12, a12, a2, a2, b12, b12, b2, b2),
i, pred j, true)
end
else
loop '(@(a12, a12, a2, a2, b12, b12, b2, b2),
i, j, overflow)
end
in
loop (!state_ref)
end
end
(*------------------------------------------------------------------*)
(* Some basic operations on two continued fractions. *)
fn
cf_add_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_sub_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i ~1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_mul_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_div_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 1, g0i2i 0), x, y)
overload + with cf_add_cf
overload - with cf_sub_cf
overload * with cf_mul_cf
overload / with cf_div_cf
(*------------------------------------------------------------------*)
fn
show (expression : string,
cf : cf)
:void =
begin
print! (expression);
print! (" => ");
print! (cf2string cf);
print! (" => ");
ignoret ($extfcall (int, "printf", "%.10lf",
(cf_convergent_double (cf, 25))));
println! ()
end
implement
main () =
begin
println! ();
print! ("A few examples.\n\n");
show ("13/11 + 1/2", r2cf @(13LL, 11LL) + r2cf @(1LL, 2LL));
show ("22/7 + 1/2", r2cf @(22LL, 7LL) + r2cf @(1LL, 2LL));
show ("13/11 - 22/7", r2cf @(13LL, 11LL) - r2cf @(22LL, 7LL));
show ("(484/49)/(22/7)",
(r2cf @(484LL, 1LL) / r2cf @(49LL, 1LL))
/ (r2cf @(22LL, 1LL) / r2cf @(7LL, 1LL)));
println! ();
print!
("Some of the following, due to truncation of their\n",
"continued fractions, are not exact but \"should\" be.\n",
"The continued fractions were truncated because integer\n",
"overflows occurred. Perhaps one could use heuristics to\n",
"elide the erroneous large terms.\n\n");
show ("sqrt(2)", sqrt2);
show ("sqrt(2)/2", sqrt2 / 2LL);
show ("1/sqrt(2)", 1LL / sqrt2);
show ("sqrt(2)*sqrt(2)", sqrt2 * sqrt2);
show ("(1/sqrt(2))*sqrt(2)", (1LL / sqrt2) * sqrt2);
show ("(1/sqrt(2))*(1/sqrt(2))", (1LL / sqrt2) * (1LL / sqrt2));
show ("sqrt(2)/sqrt(2)", sqrt2 / sqrt2);
println! ();
0
end
(*------------------------------------------------------------------*)- Output:
$ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW bivariate-continued-fraction-task-memoizing.dats -lgc && ./a.out A few examples. 13/11 + 1/2 => [1;1,2,7] => 1.6818181818 22/7 + 1/2 => [3;1,1,1,4] => 3.6428571429 13/11 - 22/7 => [-1;-1,-24,-1,-2] => -1.9610389610 (484/49)/(22/7) => [3;7] => 3.1428571429 Some of the following, due to truncation of their continued fractions, are not exact but "should" be. The continued fractions were truncated because integer overflows occurred. Perhaps one could use heuristics to elide the erroneous large terms. sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 1.4142135624 sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 0.7071067812 1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 0.7071067812 sqrt(2)*sqrt(2) => [1;1,79570259] => 1.9999999874 (1/sqrt(2))*sqrt(2) => [1;5874690117631298042] => 1.0000000000 (1/sqrt(2))*(1/sqrt(2)) => [0;1,1,79570259] => 0.5000000031 sqrt(2)/sqrt(2) => [1;5874690117631298042] => 1.0000000000
Using __int128 as the integer type
Here is essentially the same program as above, but using the GNU extension for 128-bit integers.
(*------------------------------------------------------------------*)
#include "share/atspre_staload.hats"
%{#
#include <stdint.h>
#include <limits.h>
#include <float.h>
#include <math.h>
%}
#define NIL list_nil ()
#define :: list_cons
exception gint_overflow of ()
(*------------------------------------------------------------------*)
extern fn {tk : tkind}
g0int_neginf :
() -<> g0int tk
extern fn {tk : tkind}
g0int_add_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_sub_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_mul_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
infixl ( + ) +!
infixl ( - ) -!
infixl ( * ) *!
overload neginf with g0int_neginf
overload +! with g0int_add_overflow_exn
overload -! with g0int_sub_overflow_exn
overload *! with g0int_mul_overflow_exn
(*------------------------------------------------------------------*)
(* 128-bit integers. *)
%{#
/* The most negative int128 will be treated as "neginf" or "negative
infinity". For our purposes the sign will not matter, though. */
#define neginf_int128() (((__int128) 1) << 127)
#define neg_c(x) (-(x))
#define add_c(x, y) ((x) + (y))
#define sub_c(x, y) ((x) - (y))
#define mul_c(x, y) ((x) * (y))
#define div_c(x, y) ((x) / (y))
#define mod_c(x, y) ((x) % (y))
#define eq_c(x, y) (((x) == (y)) ? atsbool_true : atsbool_false)
#define neq_c(x, y) (((x) != (y)) ? atsbool_true : atsbool_false)
#define lt_c(x, y) (((x) < (y)) ? atsbool_true : atsbool_false)
#define lte_c(x, y) (((x) <= (y)) ? atsbool_true : atsbool_false)
#define gt_c(x, y) (((x) > (y)) ? atsbool_true : atsbool_false)
#define gte_c(x, y) (((x) >= (y)) ? atsbool_true : atsbool_false)
/* GNU extensions for detection of integer overflow. */
#define add_overflow(x, y, pz) \
(__builtin_add_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define sub_overflow(x, y, pz) \
(__builtin_sub_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define mul_overflow(x, y, pz) \
(__builtin_mul_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
%}
tkindef int128_kind = "__int128" (* A GNU extension. *)
stadef int128knd = int128_kind
typedef int128_0 = g0int int128knd
typedef int128_1 (i : int) = g1int (int128knd, i)
stadef int128 = int128_1 // 2nd-select
stadef int128 = int128_0 // 1st-select
stadef Int128 = [i : int] int128_1 i
extern fn g0int_neginf_int128 : () -<> int128 = "mac#neginf_int128"
extern fn g0int_neg_int128 : int128 -<> int128 = "mac#neg_c"
extern fn g0int_add_int128 : (int128, int128) -<> int128 = "mac#add_c"
extern fn g0int_sub_int128 : (int128, int128) -<> int128 = "mac#sub_c"
extern fn g0int_mul_int128 : (int128, int128) -<> int128 = "mac#mul_c"
extern fn g0int_div_int128 : (int128, int128) -<> int128 = "mac#div_c"
extern fn g0int_mod_int128 : (int128, int128) -<> int128 = "mac#mod_c"
extern fn g0int_eq_int128 : (int128, int128) -<> bool = "mac#eq_c"
extern fn g0int_neq_int128 : (int128, int128) -<> bool = "mac#neq_c"
extern fn g0int_lt_int128 : (int128, int128) -<> bool = "mac#lt_c"
extern fn g0int_lte_int128 : (int128, int128) -<> bool = "mac#lte_c"
extern fn g0int_gt_int128 : (int128, int128) -<> bool = "mac#gt_c"
extern fn g0int_gte_int128 : (int128, int128) -<> bool = "mac#gte_c"
implement g0int_neginf<int128knd> () = g0int_neginf_int128 ()
implement g0int_neg<int128knd> x = g0int_neg_int128 x
implement g0int_add<int128knd> (x, y) = g0int_add_int128 (x, y)
implement g0int_sub<int128knd> (x, y) = g0int_sub_int128 (x, y)
implement g0int_mul<int128knd> (x, y) = g0int_mul_int128 (x, y)
implement g0int_div<int128knd> (x, y) = g0int_div_int128 (x, y)
implement g0int_mod<int128knd> (x, y) = g0int_mod_int128 (x, y)
implement g0int_eq<int128knd> (x, y) = g0int_eq_int128 (x, y)
implement g0int_neq<int128knd> (x, y) = g0int_neq_int128 (x, y)
implement g0int_lt<int128knd> (x, y) = g0int_lt_int128 (x, y)
implement g0int_lte<int128knd> (x, y) = g0int_lte_int128 (x, y)
implement g0int_gt<int128knd> (x, y) = g0int_gt_int128 (x, y)
implement g0int_gte<int128knd> (x, y) = g0int_gte_int128 (x, y)
implement g0int2int<intknd,int128knd> i = $UNSAFE.cast i
implement g0int2int<int128knd,intknd> i = $UNSAFE.cast i
implement g0int2float<int128knd,ldblknd> i = $UNSAFE.cast i
implement g0int_iseqz<int128knd> x = (x = g0i2i 0)
implement g0int_isneqz<int128knd> x = (x <> g0i2i 0)
implement g0int_isltz<int128knd> x = (x < g0i2i 0)
implement g0int_isltez<int128knd> x = (x <= g0i2i 0)
implement g0int_isgtz<int128knd> x = (x > g0i2i 0)
implement g0int_isgtez<int128knd> x = (x >= g0i2i 0)
implement g0int_abs<int128knd> x = (if isltz x then ~x else x)
local
extern fn
add_overflow_int128 :
(int128, int128, &int128? >> int128) -< !wrt > bool
= "mac#add_overflow"
extern fn
sub_overflow_int128 :
(int128, int128, &int128? >> int128) -< !wrt > bool
= "mac#sub_overflow"
extern fn
mul_overflow_int128 :
(int128, int128, &int128? >> int128) -< !wrt > bool
= "mac#mul_overflow"
in (* local *)
fn
g0int_add_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt add_overflow_int128 (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
fn
g0int_sub_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt sub_overflow_int128 (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
fn
g0int_mul_overflow_exn_int128 (x : int128, y : int128)
:<!exn> int128 =
let
var z : int128?
val overflow = $effmask_wrt mul_overflow_int128 (x, y, z)
in
if ~overflow then
z
else
$raise gint_overflow ()
end
end (* local *)
implement
g0int_add_overflow_exn<int128knd> (x, y) =
g0int_add_overflow_exn_int128 (x, y)
implement
g0int_sub_overflow_exn<int128knd> (x, y) =
g0int_sub_overflow_exn_int128 (x, y)
implement
g0int_mul_overflow_exn<int128knd> (x, y) =
g0int_mul_overflow_exn_int128 (x, y)
(*------------------------------------------------------------------*)
(* We will truncate quotients towards zero. *)
infixl ( / ) div divrem
infixl ( mod ) rem
macdef div = g0int_div
macdef rem = g0int_mod
macdef divrem (a, b) =
let
val x = ,(a)
and y = ,(b)
in
(* Optimizing C compilers will compute both quotient and remainder
at the same time. *)
@(x \g0int_div y, x \g0int_mod y)
end
(*------------------------------------------------------------------*)
(* Continued fractions.
cf_generator tk -- A closure that produces terms of type g0int tk,
sequentially.
cf tk -- A structure from which one can get the ith
term of a continued fraction. It gets terms
from a cf_generator tk. *)
typedef integer = int128
stadef integerknd = int128knd
typedef cf_generator = () -<cloref1> integer
local
typedef _cf (terminated : bool,
m : int,
n : int) =
[m <= n]
@{
terminated = bool terminated, (* No more terms? *)
m = size_t m, (* The number of terms computed so far. *)
n = size_t n, (* The size of memo storage.*)
memo = arrayref (integer, n), (* Memoized terms. *)
gen = cf_generator (* A thunk to generate terms. *)
}
typedef _cf (m : int) =
[terminated : bool]
[n : int | m <= n]
_cf (terminated, m, n)
typedef _cf =
[m : int]
_cf m
fn
_cf_get_more_terms
{terminated : bool}
{m : int}
{n : int}
{needed : int | m <= needed; needed <= n}
(cf : _cf (terminated, m, n),
needed : size_t needed)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)
macdef memo = cf.memo
fun
loop {i : int | m <= i; i <= needed}
.<needed - i>.
(i : size_t i)
: [m1 : int | m <= m1; m1 <= needed]
[n1 : int | m1 <= n1]
_cf (m1 < needed, m1, n1) =
if i = needed then
@{
terminated = false,
m = needed,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
else
let
val term = (cf.gen) ()
in
if term <> neginf<integerknd> () then
begin
memo[i] := term;
loop (succ i)
end
else
@{
terminated = true,
m = i,
n = (cf.n),
memo = memo,
gen = (cf.gen)
}
end
in
loop (cf.m)
end
fn
_cf_update {terminated : bool}
{m : int}
{n : int | m <= n}
{needed : int}
(cf : _cf (terminated, m, n),
needed : size_t needed)
: [terminated1 : bool]
[m1 : int | m <= m1]
[n1 : int | m1 <= n1]
_cf (terminated1, m1, n1) =
let
prval () = lemma_g1uint_param (cf.m)
macdef memo = cf.memo
macdef gen = cf.gen
in
if (cf.terminated) then
cf
else if needed <= (cf.m) then
cf
else if needed <= (cf.n) then
_cf_get_more_terms (cf, needed)
else
let (* Provides twice the room needed. *)
val n1 = needed + needed
val memo1 = arrayref_make_elt (n1, g0i2i 0)
val () =
let
var i : [i : nat] size_t i
in
for (i := i2sz 0; i < (cf.m); i := succ i)
memo1[i] := memo[i]
end
val cf1 =
@{
terminated = false,
m = (cf.m),
n = n1,
memo = memo1,
gen = (cf.gen)
}
in
_cf_get_more_terms (cf1, needed)
end
end
in (* local *)
typedef cf = ref _cf
extern fn cf_make : cf_generator -> cf
extern fn cf_get_at_size : {i : int} (cf, size_t i) -> integer
extern fn cf_get_at_int : {i : nat} (cf, int i) -> integer
(* The precedence of the overloads has to exceed that of ref[] *)
overload cf_get_at with cf_get_at_size of 1
overload cf_get_at with cf_get_at_int of 1
overload [] with cf_get_at of 1
implement
cf_make gen =
let
#ifndef CF_START_SIZE #then
#define CF_START_SIZE 32
#endif
in
ref
@{
terminated = false,
m = i2sz 0,
n = i2sz CF_START_SIZE,
memo = arrayref_make_elt (i2sz CF_START_SIZE, g0i2i 0),
gen = gen
}
end
implement
cf_get_at_size (cf, i) =
let
prval () = lemma_g1uint_param i
val cf1 = _cf_update (!cf, succ i)
in
!cf := cf1;
if i < (cf1.m) then
arrayref_get_at<integer> (cf1.memo, i)
else
neginf<integerknd> ()
end
implement cf_get_at_int (cf, i) = cf_get_at_size (cf, g1i2u i)
end (* local *)
(*------------------------------------------------------------------*)
(* Make a string from an int128. *)
fn
int128_2string (i : int128) : string =
let
fun
loop (i : int128, accum : List0 charNZ) : List0 charNZ =
if iseqz i then
accum
else
let
val @(i, digit) = i divrem (g0i2i 10)
val digit = (g0i2i digit) : int
val digit = g1ofg0 (int2char0 (digit + (char2int0 '0')))
val () = assertloc (digit <> '\0')
in
loop (i, digit :: accum)
end
val minus_sign : Char = '-'
val () = assertloc (minus_sign <> '\0')
in
if iseqz i then
"0"
else if i = neginf<integerknd> () then
"neginf"
else if isltz i then
strnptr2string (string_implode (minus_sign :: (loop (~i, NIL))))
else
strnptr2string (string_implode (loop (i, NIL)))
end
implement tostring_val<integer> i = $effmask_all int128_2string i
implement fprint_val<integer> (f, i) = fprint! (f, tostring_val<integer> i)
fn fprint_integer (f : FILEref, i : integer) = fprint_val<integer> (f, i)
fn print_integer (i : integer) = fprint_integer (stdout_ref, i)
fn prerr_integer (i : integer) = fprint_integer (stderr_ref, i)
overload fprint with fprint_integer
overload print with print_integer
overload prerr with prerr_integer
(*------------------------------------------------------------------*)
(* Converting a continued fraction to a string. *)
extern fn cf2string_with_default_max_terms : cf -> string
extern fn cf2string_given_max_terms : (cf, Size_t) -> string
overload cf2string with cf2string_with_default_max_terms
overload cf2string with cf2string_given_max_terms
val cf2string_default_max_terms : ref Size_t = ref (i2sz 20)
implement
cf2string_with_default_max_terms cf =
cf2string_given_max_terms (cf, !cf2string_default_max_terms)
implement
cf2string_given_max_terms (cf, max_terms) =
let
fun
loop (i : Size_t,
accum : List0 string)
: List0 string =
let
val term = cf[i]
in
if i = max_terms then
begin
if term = neginf<integerknd> () then
"]" :: accum
else
",...]" :: accum
end
else if term = neginf<integerknd> () then
"]" :: accum
else
let
val separator =
(if i = i2sz 0 then
""
else if i = i2sz 1 then
";"
else
",")
and term_str = tostring_val<integer> term
in
loop (succ i, term_str :: separator :: accum)
end
end
val string_lst = loop (i2sz 0, "[" :: NIL)
in
strptr2string (stringlst_concat (list_vt2t (reverse string_lst)))
end
(*------------------------------------------------------------------*)
(* Compute the ith convergent as a floating-point number. Convergents
are indexed starting at 0. *)
extern fn
cf_convergent_ldouble_size :
{i : int}
(cf, size_t i) -> ldouble
extern fn
cf_convergent_ldouble_int :
{i : nat}
(cf, int i) -> ldouble
overload cf_convergent_ldouble with cf_convergent_ldouble_size
overload cf_convergent_ldouble with cf_convergent_ldouble_int
implement
cf_convergent_ldouble_size {i} (cf, i) =
let
prval () = lemma_g1uint_param i
fun
loop {j : nat | j <= i + 1}
.<(i + 1) - j>.
(j : size_t j,
accum : ldouble)
: ldouble =
if j = succ i then
accum
else
let
val term = cf[i - j]
val x =
(if term = neginf<integerknd> () then
$extval (ldouble, "INFINITY")
else
g0i2f term) : ldouble
in
loop (succ j, x + (1.0L / accum))
end
in
loop (i2sz 0, 0.0L)
end
implement
cf_convergent_ldouble_int (cf, i) =
cf_convergent_ldouble_size (cf, g1i2u i)
(*------------------------------------------------------------------*)
(* A continued fraction for the square root of two. *)
val sqrt2 : cf =
cf_make
let
val first : ref bool = ref true
in
lam () =<cloref1>
let
val fst = !first
in
!first := false;
(if fst then g0i2i 1 else g0i2i 2) : integer
end
end
(*------------------------------------------------------------------*)
(* The continued fraction for a rational number. *)
typedef ratnum = @(integer, integer)
fn
r2cf (ratnum : ratnum) : cf =
cf_make
let
val ratnum_ref : ref ratnum = ref ratnum
in
lam () =<cloref1>
let
val @(n, d) = !ratnum_ref
in
if iseqz d then
neginf<integerknd> ()
else
let
val @(q, r) = n divrem d
in
!ratnum_ref := @(d, r);
q
end
end
end
(*------------------------------------------------------------------*)
(* Application of a homographic function to a continued fraction. *)
typedef ng4 = @(integer, integer, integer, integer)
fn
apply_ng4 (ng4 : ng4, x : cf) : cf =
cf_make
let
val state : ref @(ng4, Size_t) = ref @(ng4, i2sz 0)
in
lam () =<cloref1>
let
fnx
loop (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
if (iseqz b1) * (iseqz b) then
neginf<integerknd> ()
else if (iseqz b1) + (iseqz b) then
absorb_term (a1, a, b1, b, i)
else
let
val q1 = a1 div b1
and q = a div b
in
if q1 <> q then
absorb_term (a1, a, b1, b, i)
else
begin
!state :=
@(@(b1, b, a1 - (b1 * q), a - (b * q)), i);
q
end
end
and
absorb_term (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
let
val term = x[i]
in
if term <> neginf<integerknd> () then
loop (a + (a1 * term), a1,
b + (b1 * term), b1, succ i)
else
loop (a1, a1, b1, b1, succ i)
end
val @(@(a1, a, b1, b), i) = !state
in
loop (a1, a, b1, b, i)
end
end
(*------------------------------------------------------------------*)
(* Some basic operations involving only one continued fraction. *)
fn
cf_neg (x : cf) : cf =
apply_ng4 (@(g0i2i ~1, g0i2i 0, g0i2i 0, g0i2i 1), x)
fn
cf_add_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(d, n, g0i2i 0, d), x)
end
fn
ratnum_add_cf (x : ratnum, y : cf) : cf =
cf_add_ratnum (y, x)
fn
cf_add_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(y, g0i2i 1))
fn
integer_add_cf (x : integer, y : cf) : cf =
cf_add_ratnum (y, @(x, g0i2i 1))
fn
cf_sub_ratnum (x : cf, y : ratnum) : cf =
cf_add_ratnum (x, @(~(y.0), (y.1)))
fn
ratnum_sub_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(~d, n, g0i2i 0, d), y)
end
fn
cf_sub_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(~y, g0i2i 1))
fn
integer_sub_cf (x : integer, y : cf) : cf =
ratnum_sub_cf (@(x, g0i2i 1), y)
fn
cf_mul_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), x)
end
fn
ratnum_mul_cf (x : ratnum, y : cf) : cf =
cf_mul_ratnum (y, x)
fn
cf_mul_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(y, g0i2i 1))
fn
integer_mul_cf (x : integer, y : cf) : cf =
cf_mul_ratnum (y, @(x, g0i2i 1))
fn
cf_div_ratnum (x : cf, y : ratnum) : cf =
cf_mul_ratnum (x, @(y.1, y.0))
fn
ratnum_div_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(g0i2i 0, n, d, g0i2i 0), y)
end
fn
cf_div_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(g0i2i 1, y))
fn
integer_div_cf (x : integer, y : cf) : cf =
ratnum_div_cf (@(x, g0i2i 1), y)
overload ~ with cf_neg
overload + with cf_add_ratnum
overload + with ratnum_add_cf
overload + with cf_add_integer
overload + with integer_add_cf
overload - with cf_sub_ratnum
overload - with ratnum_sub_cf
overload - with cf_sub_integer
overload - with integer_sub_cf
overload * with cf_mul_ratnum
overload * with ratnum_mul_cf
overload * with cf_mul_integer
overload * with integer_mul_cf
overload / with cf_div_ratnum
overload / with ratnum_div_cf
overload / with cf_div_integer
overload / with integer_div_cf
(*------------------------------------------------------------------*)
(* Application of a bihomographic function to a continued fraction. *)
typedef ng8 = @(integer, integer, integer, integer,
integer, integer, integer, integer)
fn
apply_ng8 (ng8 : ng8, x : cf, y : cf) : cf =
cf_make
let
typedef state = '(ng8, Size_t, Size_t, bool)
val state_ref : ref state = ref '(ng8, i2sz 0, i2sz 0, false)
in
lam () =<cloref1>
let
fnx
loop (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
in
if iseqz b12 * iseqz b1 * iseqz b2 * iseqz b then
neginf<integerknd> ()
else if iseqz b * iseqz b2 then
absorb_x_term state
else if iseqz b + iseqz b2 then
absorb_y_term state
else if iseqz b1 then
absorb_x_term state
else
let
val @(q, r) = a divrem b // a/b = q + r/b
and @(q1, r1) = a1 divrem b1 // a1/b1 = q1 + r1/b1
and @(q2, r2) = a2 divrem b2 // a2/b2 = q2 + r2/b2
val q1_diff = abs (q1 -! q)
and q2_diff = abs (q2 -! q)
in
if q1_diff > q2_diff then
absorb_x_term state
else if q1_diff < q2_diff then
absorb_y_term state
else if iseqz b12 || q <> a12 div b12 then
let
(* Because q = q1 = q2, the following also are
true:
a1/b1 - a/b = r1/b1 - r/b
a2/b2 - a/b = r2/b2 - r/b
But now the numerators are smaller, and so
overflow is less likely to occur when we put
terms over a common denominator. *)
(* Put r, r1, r2 over a common denominator,
b*b1*b2. *)
val n = r *! b1 *! b2
and n1 = r1 *! b *! b2
and n2 = r2 *! b *! b1
in
if abs (n1 -! n) > abs (n2 -! n) then
absorb_x_term state
else
absorb_y_term state
end
else
begin
!state_ref :=
'(@(b12, b1, b2, b,
a12 -! (b12 *! q), a1 -! (b1 *! q),
a2 -! (b2 *! q), a -! (b *! q)),
i, j, overflow);
q
end
end
end
and
absorb_x_term (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
val term = (if overflow then neginf<integerknd> () else x[i])
and i = succ i
in
if term <> neginf<integerknd> () then
begin
try
loop '(@(a2 +! (a12 *! term),
a +! (a1 *! term), a12, a1,
b2 +! (b12 *! term),
b +! (b1 *! term), b12, b1),
i, j, overflow)
with
| ~ gint_overflow () =>
loop '(@(a12, a1, a12, a1, b12, b1, b12, b1),
pred i, j, true)
end
else
loop '(@(a12, a1, a12, a1, b12, b1, b12, b1),
i, j, overflow)
end
and
absorb_y_term (state : state) : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
i, j, overflow) = state
val term = (if overflow then neginf<integerknd> () else y[j])
and j = succ j
in
if term <> neginf<integerknd> () then
begin
try
loop '(@(a1 +! (a12 *! term), a12,
a +! (a2 *! term), a2,
b1 +! (b12 *! term), b12,
b +! (b2 *! term), b2),
i, j, overflow)
with
| ~ gint_overflow () =>
loop '(@(a12, a12, a2, a2, b12, b12, b2, b2),
i, pred j, true)
end
else
loop '(@(a12, a12, a2, a2, b12, b12, b2, b2),
i, j, overflow)
end
in
loop (!state_ref)
end
end
(*------------------------------------------------------------------*)
(* Some basic operations on two continued fractions. *)
fn
cf_add_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_sub_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i ~1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_mul_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_div_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 1, g0i2i 0), x, y)
overload + with cf_add_cf
overload - with cf_sub_cf
overload * with cf_mul_cf
overload / with cf_div_cf
(*------------------------------------------------------------------*)
fn
show (expression : string,
cf : cf)
:void =
begin
print! (expression);
print! (" => ");
print! (cf2string cf);
print! (" => ");
ignoret ($extfcall (int, "printf", "%.15Lf",
(cf_convergent_ldouble (cf, 50))));
println! ()
end
implement
main () =
let
macdef i2integer i = (g0i2i ,(i)) : integer
macdef I = i2integer
in
println! ();
print! ("A few examples.\n\n");
show ("13/11 + 1/2", r2cf @(I 13, I 11) + r2cf @(I 1, I 2));
show ("22/7 + 1/2", r2cf @(I 22, I 7) + r2cf @(I 1, I 2));
show ("13/11 - 22/7", r2cf @(I 13, I 11) - r2cf @(I 22, I 7));
show ("(484/49)/(22/7)",
(r2cf @(I 484, I 1) / r2cf @(I 49, I 1))
/ (r2cf @(I 22, I 1) / r2cf @(I 7, I 1)));
println! ();
print!
("Some of the following, due to truncation of their\n",
"continued fractions, are not exact but \"should\" be.\n",
"The continued fractions were truncated because integer\n",
"overflows occurred. Perhaps one could use heuristics to\n",
"elide the erroneous large terms.\n\n");
show ("sqrt(2)", sqrt2);
show ("sqrt(2)/2", sqrt2 / I 2);
show ("1/sqrt(2)", (I 1) / sqrt2);
show ("sqrt(2)*sqrt(2)", sqrt2 * sqrt2);
show ("(1/sqrt(2))*sqrt(2)", ((I 1) / sqrt2) * sqrt2);
show ("(1/sqrt(2))*(1/sqrt(2))", ((I 1) / sqrt2) * ((I 1) / sqrt2));
show ("sqrt(2)/sqrt(2)", sqrt2 / sqrt2);
println! ();
0
end
(*------------------------------------------------------------------*)- Output:
$ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW bivariate-continued-fraction-task-int128.dats -lgc && ./a.out A few examples. 13/11 + 1/2 => [1;1,2,7] => 1.681818181818182 22/7 + 1/2 => [3;1,1,1,4] => 3.642857142857143 13/11 - 22/7 => [-1;-1,-24,-1,-2] => -1.961038961038961 (484/49)/(22/7) => [3;7] => 3.142857142857143 Some of the following, due to truncation of their continued fractions, are not exact but "should" be. The continued fractions were truncated because integer overflows occurred. Perhaps one could use heuristics to elide the erroneous large terms. sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 1.414213562373095 sqrt(2)/2 => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 0.707106781186548 1/sqrt(2) => [0;1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] => 0.707106781186548 sqrt(2)*sqrt(2) => [1;1,105964828892651] => 1.999999999999991 (1/sqrt(2))*sqrt(2) => [1;80866608600094285540655993317041165854] => 1.000000000000000 (1/sqrt(2))*(1/sqrt(2)) => [0;1,1,105964828892651] => 0.500000000000002 sqrt(2)/sqrt(2) => [1;80866608600094285540655993317041165854] => 1.000000000000000
C++
Uses matrixNG, NG_4 and NG from Continued_fraction/Arithmetic/G(matrix_NG,_Contined_Fraction_N)#C++, and r2cf from Continued_fraction/Arithmetic/Construct_from_rational_number#C++
/* Implement matrix NG
Nigel Galloway, February 12., 2013
*/
class NG_8 : public matrixNG {
private: int a12, a1, a2, a, b12, b1, b2, b, t;
double ab, a1b1, a2b2, a12b12;
const int chooseCFN(){return fabs(a1b1-ab) > fabs(a2b2-ab)? 0 : 1;}
const bool needTerm() {
if (b1==0 and b==0 and b2==0 and b12==0) return false;
if (b==0){cfn = b2==0? 0:1; return true;} else ab = ((double)a)/b;
if (b2==0){cfn = 1; return true;} else a2b2 = ((double)a2)/b2;
if (b1==0){cfn = 0; return true;} else a1b1 = ((double)a1)/b1;
if (b12==0){cfn = chooseCFN(); return true;} else a12b12 = ((double)a12)/b12;
thisTerm = (int)ab;
if (thisTerm==(int)a1b1 and thisTerm==(int)a2b2 and thisTerm==(int)a12b12){
t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm; t=a2; a2=b2; b2=t-b2*thisTerm; t=a12; a12=b12; b12=t-b12*thisTerm;
haveTerm = true; return false;
}
cfn = chooseCFN();
return true;
}
void consumeTerm(){if(cfn==0){a=a1; a2=a12; b=b1; b2=b12;} else{a=a2; a1=a12; b=b2; b1=b12;}}
void consumeTerm(int n){
if(cfn==0){t=a; a=a1; a1=t+a1*n; t=a2; a2=a12; a12=t+a12*n; t=b; b=b1; b1=t+b1*n; t=b2; b2=b12; b12=t+b12*n;}
else{t=a; a=a2; a2=t+a2*n; t=a1; a1=a12; a12=t+a12*n; t=b; b=b2; b2=t+b2*n; t=b1; b1=b12; b12=t+b12*n;}
}
public:
NG_8(int a12, int a1, int a2, int a, int b12, int b1, int b2, int b): a12(a12), a1(a1), a2(a2), a(a), b12(b12), b1(b1), b2(b2), b(b){
}};
Testing
[3;7] + [0;2]
int main() {
NG_8 a(0,1,1,0,0,0,0,1);
r2cf n2(22,7);
r2cf n1(1,2);
for(NG n(&a, &n1, &n2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
NG_4 a3(2,1,0,2);
r2cf n3(22,7);
for(NG n(&a3, &n3); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 1 1 1 4 3 1 1 1 4
[1:5,2] * [3;7]
int main() {
NG_8 b(1,0,0,0,0,0,0,1);
r2cf b1(13,11);
r2cf b2(22,7);
for(NG n(&b, &b1, &b2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(NG n(&a, &b2, &b1); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(286,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 1 2 2 3 1 2 2
[1:5,2] - [3;7]
int main() {
NG_8 c(0,1,-1,0,0,0,0,1);
r2cf c1(13,11);
r2cf c2(22,7);
for(NG n(&c, &c1, &c2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(-151,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
-1 -1 -24 -1 -2 -1 -1 -24 -1 -2
Divide [] by [3;7]
int main() {
NG_8 d(0,1,0,0,0,0,1,0);
r2cf d1(22*22,7*7);
r2cf d2(22,7);
for(NG n(&d, &d1, &d2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 7
([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])
int main() {
r2cf a1(2,7);
r2cf a2(13,11);
NG_8 na(0,1,1,0,0,0,0,1);
NG A(&na, &a1, &a2); //[0;3,2] + [1;5,2]
r2cf b1(2,7);
r2cf b2(13,11);
NG_8 nb(0,1,-1,0,0,0,0,1);
NG B(&nb, &b1, &b2); //[0;3,2] - [1;5,2]
NG_8 nc(1,0,0,0,0,0,0,1); //A*B
for(NG n(&nc, &A, &B); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(2,7); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(13,11); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(-7797,5929); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
Go
Adding to the existing package from the
Continued_fraction/Arithmetic/Construct_from_rational_number#Go
task, re-uses cf.go and rat.go as given in that task.
File ng8.go:
package cf
import "math"
// A 2×4 matix:
// [ a₁₂ a₁ a₂ a ]
// [ b₁₂ b₁ b₂ b ]
//
// which when "applied" to two continued fractions N1 and N2
// gives a new continued fraction z such that:
//
// a₁₂ * N1 * N2 + a₁ * N1 + a₂ * N2 + a
// z = -------------------------------------------
// b₁₂ * N1 * N2 + b₁ * N1 + b₂ * N2 + b
//
// Examples:
// NG8{0,1,1,0, 0,0,0,1} gives N1 + N2
// NG8{0,1,-1,0, 0,0,0,1} gives N1 - N2
// NG8{1,0,0,0, 0,0,0,1} gives N1 * N2
// NG8{0,1,0,0, 0,0,1,0} gives N1 / N2
// NG8{21,-15,28,-20, 0,0,0,1} gives 21*N1*N2 -15*N1 +28*N2 -20
// which is (3*N1 + 4) * (7*N2 - 5)
type NG8 struct {
A12, A1, A2, A int64
B12, B1, B2, B int64
}
// Basic identities as NG8 matrices.
var (
NG8Add = NG8{0, 1, 1, 0, 0, 0, 0, 1}
NG8Sub = NG8{0, 1, -1, 0, 0, 0, 0, 1}
NG8Mul = NG8{1, 0, 0, 0, 0, 0, 0, 1}
NG8Div = NG8{0, 1, 0, 0, 0, 0, 1, 0}
)
func (ng *NG8) needsIngest() bool {
if ng.B12 == 0 || ng.B1 == 0 || ng.B2 == 0 || ng.B == 0 {
return true
}
x := ng.A / ng.B
return ng.A1/ng.B1 != x || ng.A2/ng.B2 != x && ng.A12/ng.B12 != x
}
func (ng *NG8) isDone() bool {
return ng.B12 == 0 && ng.B1 == 0 && ng.B2 == 0 && ng.B == 0
}
func (ng *NG8) ingestWhich() bool { // true for N1, false for N2
if ng.B == 0 && ng.B2 == 0 {
return true
}
if ng.B == 0 || ng.B2 == 0 {
return false
}
x1 := float64(ng.A1) / float64(ng.B1)
x2 := float64(ng.A2) / float64(ng.B2)
x := float64(ng.A) / float64(ng.B)
return math.Abs(x1-x) > math.Abs(x2-x)
}
func (ng *NG8) ingest(isN1 bool, t int64) {
if isN1 {
// [ a₁₂ a₁ a₂ a ] becomes [ a₂+a₁₂*t a+a₁*t a₁₂ a₁]
// [ b₁₂ b₁ b₂ b ] [ b₂+b₁₂*t b+b₁*t b₁₂ b₁]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.A2+ng.A12*t, ng.A+ng.A1*t, ng.A12, ng.A1,
ng.B2+ng.B12*t, ng.B+ng.B1*t, ng.B12, ng.B1
} else {
// [ a₁₂ a₁ a₂ a ] becomes [ a₁+a₁₂*t a₁₂ a+a₂*t a₂]
// [ b₁₂ b₁ b₂ b ] [ b₁+b₁₂*t b₁₂ b+b₂*t b₂]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.A1+ng.A12*t, ng.A12, ng.A+ng.A2*t, ng.A2,
ng.B1+ng.B12*t, ng.B12, ng.B+ng.B2*t, ng.B2
}
}
func (ng *NG8) ingestInfinite(isN1 bool) {
if isN1 {
// [ a₁₂ a₁ a₂ a ] becomes [ a₁₂ a₁ a₁₂ a₁ ]
// [ b₁₂ b₁ b₂ b ] [ b₁₂ b₁ b₁₂ b₁ ]
ng.A2, ng.A, ng.B2, ng.B =
ng.A12, ng.A1,
ng.B12, ng.B1
} else {
// [ a₁₂ a₁ a₂ a ] becomes [ a₁₂ a₁₂ a₂ a₂ ]
// [ b₁₂ b₁ b₂ b ] [ b₁₂ b₁₂ b₂ b₂ ]
ng.A1, ng.A, ng.B1, ng.B =
ng.A12, ng.A2,
ng.B12, ng.B2
}
}
func (ng *NG8) egest(t int64) {
// [ a₁₂ a₁ a₂ a ] becomes [ b₁₂ b₁ b₂ b ]
// [ b₁₂ b₁ b₂ b ] [ a₁₂-b₁₂*t a₁-b₁*t a₂-b₂*t a-b*t ]
ng.A12, ng.A1, ng.A2, ng.A,
ng.B12, ng.B1, ng.B2, ng.B =
ng.B12, ng.B1, ng.B2, ng.B,
ng.A12-ng.B12*t, ng.A1-ng.B1*t, ng.A2-ng.B2*t, ng.A-ng.B*t
}
// ApplyTo "applies" the matrix `ng` to the continued fractions
// `N1` and `N2`, returning the resulting continued fraction.
// After ingesting `limit` terms without any output terms the resulting
// continued fraction is terminated.
func (ng NG8) ApplyTo(N1, N2 ContinuedFraction, limit int) ContinuedFraction {
return func() NextFn {
next1, next2 := N1(), N2()
done := false
sinceEgest := 0
return func() (int64, bool) {
if done {
return 0, false
}
for ng.needsIngest() {
sinceEgest++
if sinceEgest > limit {
done = true
return 0, false
}
isN1 := ng.ingestWhich()
next := next2
if isN1 {
next = next1
}
if t, ok := next(); ok {
ng.ingest(isN1, t)
} else {
ng.ingestInfinite(isN1)
}
}
sinceEgest = 0
t := ng.A / ng.B
ng.egest(t)
done = ng.isDone()
return t, true
}
}
}
File ng8_test.go:
package cf
import "fmt"
func ExampleNG8() {
cases := [...]struct {
op string
r1, r2 Rat
ng NG8
}{
{"+", Rat{22, 7}, Rat{1, 2}, NG8Add},
{"*", Rat{13, 11}, Rat{22, 7}, NG8Mul},
{"-", Rat{13, 11}, Rat{22, 7}, NG8Sub},
{"/", Rat{22 * 22, 7 * 7}, Rat{22, 7}, NG8Div},
}
for _, tc := range cases {
n1 := tc.r1.AsContinuedFraction()
n2 := tc.r2.AsContinuedFraction()
z := tc.ng.ApplyTo(n1, n2, 1000)
fmt.Printf("%v %s %v is %v %v %v gives %v\n",
tc.r1, tc.op, tc.r2,
tc.ng, n1, n2, z,
)
}
z := NG8Mul.ApplyTo(Sqrt2, Sqrt2, 1000)
fmt.Println("√2 * √2 =", z)
// Output:
// 22/7 + 1/2 is {0 1 1 0 0 0 0 1} [3; 7] [0; 2] gives [3; 1, 1, 1, 4]
// 13/11 * 22/7 is {1 0 0 0 0 0 0 1} [1; 5, 2] [3; 7] gives [3; 1, 2, 2]
// 13/11 - 22/7 is {0 1 -1 0 0 0 0 1} [1; 5, 2] [3; 7] gives [-1; -1, -24, -1, -2]
// 484/49 / 22/7 is {0 1 0 0 0 0 1 0} [9; 1, 7, 6] [3; 7] gives [3; 7]
// √2 * √2 = [1; 0, 1]
}
- Output:
(Note, [1; 0, 1] = 1 + 1 / (0 + 1/1 ) = 2, so the answer is correct, it should however be normalised to the more reasonable form of [2].)
22/7 + 1/2 is {0 1 1 0 0 0 0 1} [3; 7] [0; 2] gives [3; 1, 1, 1, 4]
13/11 * 22/7 is {1 0 0 0 0 0 0 1} [1; 5, 2] [3; 7] gives [3; 1, 2, 2]
13/11 - 22/7 is {0 1 -1 0 0 0 0 1} [1; 5, 2] [3; 7] gives [-1; -1, -24, -1, -2]
484/49 / 22/7 is {0 1 0 0 0 0 1 0} [9; 1, 7, 6] [3; 7] gives [3; 7]
√2 * √2 = [1; 0, 1]
Julia
abstract type MatrixNG end
mutable struct NG4 <: MatrixNG
cfn::Int
thisterm::Int
haveterm::Bool
a1::Int
a::Int
b1::Int
b::Int
NG4(a1, a, b1, b) = new(0, 0, false, a1, a, b1, b)
end
mutable struct NG8 <: MatrixNG
cfn::Int
thisterm::Int
haveterm::Bool
a12::Int
a1::Int
a2::Int
a::Int
b12::Int
b1::Int
b2::Int
b::Int
NG8(a12, a1, a2, a, b12, b1, b2, b) = new(0, 0, false, a12, a1, a2, a, b12, b1, b2, b)
end
function needterm(m::NG4)::Bool
m.b1 == m.b == 0 && return false
(m.b1 == 0 || m.b == 0) && return true
(m.thisterm = m.a ÷ m.b) != m.a1 ÷ m.b1 && return true
m.a, m.b = m.b, m.a - m.b * m.thisterm
m.a1, m.b1, m.haveterm = m.b1, m.a1 - m.b1 * m.thisterm, true
return false
end
consumeterm(m::NG4) = (m.a = m.a1; m.b = m.b1)
function consumeterm(m::NG4, n)
m.a, m.a1 = m.a1, m.a + m.a1 * n
m.b, m.b1 = m.b1, m.b + m.b1 * n
end
function needterm(m::NG8)::Bool
m.b1 == m.b == m.b2 == m.b12 == 0 && return false
if m.b == 0
m.cfn = m.b2 == 0 ? 0 : 1
return true
elseif m.b2 == 0
m.cfn = 1
return true
elseif m.b1 == 0
m.cfn = 0
return true
end
ab = m.a / m.b
a1b1 = m.a1 / m.b1
a2b2 = m.a2 / m.b2
if m.b12 == 0
m.cfn = abs(a1b1 - ab) > abs(a2b2 - ab) ? 0 : 1
return true
end
m.thisterm = m.a ÷ m.b
if m.thisterm == m.a1 ÷ m.b1 == m.a2 ÷ m.b2 == m.a12 ÷ m.b12
m.a, m.b = m.b, m.a - m.b * m.thisterm
m.a1, m.b1 = m.b1, m.a1 - m.b1 * m.thisterm
m.a2, m.b2 = m.b2, m.a2 - m.b2 * m.thisterm
m.a12, m.b12, m.haveterm = m.b12, m.a12 - m.b12 * m.thisterm, true
return false
end
m.cfn = abs(a1b1 - ab) > abs(a2b2 - ab) ? 0 : 1
return true
end
function consumeterm(m::NG8)
if m.cfn == 0
m.a, m.a2 = m.a1, m.a12
m.b, m.b2 = m.b1, m.b12
else
m.a, m.a1 = m.a2, m.a12
m.b, m.b1 = m.b2, m.b12
end
end
function consumeterm(m::NG8, n)
if m.cfn == 0
m.a, m.a1 = m.a1, m.a + m.a1 * n
m.a2, m.a12 = m.a12, m.a2 + m.a12 * n
m.b, m.b1 = m.b1, m.b + m.b1 * n
m.b2, m.b12 = m.b12, m.b2 + m.b12 * n
else
m.a, m.a2 = m.a2, m.a + m.a2 * n
m.a1, m.a12 = m.a12, m.a1 + m.a12 * n
m.b, m.b2 = m.b2, m.b + m.b2 * n
m.b1, m.b12 = m.b12, m.b1 + m.b12 * n
end
end
abstract type ContinuedFraction end
mutable struct R2cf <: ContinuedFraction
n1::Int
n2::Int
end
function nextterm(x::R2cf)
term = x.n1 ÷ x.n2
x.n1, x.n2 = x.n2, x.n1 - term * x.n2
return term
end
moreterms(x::R2cf) = abs(x.n2) > 0
mutable struct NG <: ContinuedFraction
ng::MatrixNG
n::Vector{ContinuedFraction}
end
NG(ng, n1::ContinuedFraction) = NG(ng, [n1])
NG(ng, n1::ContinuedFraction, n2) = NG(ng, [n1, n2])
nextterm(x::NG) = (x.ng.haveterm = false; x.ng.thisterm)
function moreterms(x::NG)::Bool
while needterm(x.ng)
if moreterms(x.n[x.ng.cfn + 1])
consumeterm(x.ng, nextterm(x.n[x.ng.cfn + 1]))
else
consumeterm(x.ng)
end
end
return x.ng.haveterm
end
function testcfs()
function test(desc, cfs)
println("TESTING -> $desc")
for cf in cfs
while moreterms(cf)
print(nextterm(cf), " ")
end
println()
end
println()
end
a = NG8(0, 1, 1, 0, 0, 0, 0, 1)
n2 = R2cf(22, 7)
n1 = R2cf(1, 2)
a3 = NG4(2, 1, 0, 2)
n3 = R2cf(22, 7)
test("[3;7] + [0;2]", [NG(a, n1, n2), NG(a3, n3)])
b = NG8(1, 0, 0, 0, 0, 0, 0, 1)
b1 = R2cf(13, 11)
b2 = R2cf(22, 7)
test("[1;5,2] * [3;7]", [NG(b, b1, b2), R2cf(286, 77)])
c = NG8(0, 1, -1, 0, 0, 0, 0, 1)
c1 = R2cf(13, 11)
c2 = R2cf(22, 7)
test("[1;5,2] - [3;7]", [NG(c, c1, c2), R2cf(-151, 77)])
d = NG8(0, 1, 0, 0, 0, 0, 1, 0)
d1 = R2cf(22 * 22, 7 * 7)
d2 = R2cf(22, 7)
test("Divide [] by [3;7]", [NG(d, d1, d2)])
na = NG8(0, 1, 1, 0, 0, 0, 0, 1)
a1 = R2cf(2, 7)
a2 = R2cf(13, 11)
aa = NG(na, a1, a2)
nb = NG8(0, 1, -1, 0, 0, 0, 0, 1)
b3 = R2cf(2, 7)
b4 = R2cf(13, 11)
bb = NG(nb, b3, b4)
nc = NG8(1, 0, 0, 0, 0, 0, 0, 1)
desc = "([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])"
test(desc, [NG(nc, aa, bb), R2cf(-7797, 5929)])
end
testcfs()
- Output:
TESTING -> [3;7] + [0;2] 3 1 1 1 4 3 1 1 1 4 TESTING -> [1;5,2] * [3;7] 3 1 2 2 3 1 2 2 TESTING -> [1;5,2] - [3;7] -1 -1 -24 -1 -2 -1 -1 -24 -1 -2 TESTING -> Divide [] by [3;7] 3 7 TESTING -> ([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) -1 -3 -5 -1 -2 -1 -26 -3 -1 -3 -5 -1 -2 -1 -26 -3
Kotlin
The C++ entry uses a number of classes which have been coded in other "Continued Fraction" tasks. I've pulled all these into my Kotlin translation and unified the tests so that the whole thing can, hopefully, be understood and run as a single program.
// version 1.2.10
import kotlin.math.abs
abstract class MatrixNG {
var cfn = 0
var thisTerm = 0
var haveTerm = false
abstract fun consumeTerm()
abstract fun consumeTerm(n: Int)
abstract fun needTerm(): Boolean
}
class NG4(
var a1: Int, var a: Int, var b1: Int, var b: Int
) : MatrixNG() {
private var t = 0
override fun needTerm(): Boolean {
if (b1 == 0 && b == 0) return false
if (b1 == 0 || b == 0) return true
thisTerm = a / b
if (thisTerm == a1 / b1) {
t = a; a = b; b = t - b * thisTerm
t = a1; a1 = b1; b1 = t - b1 * thisTerm
haveTerm = true
return false
}
return true
}
override fun consumeTerm() {
a = a1
b = b1
}
override fun consumeTerm(n: Int) {
t = a; a = a1; a1 = t + a1 * n
t = b; b = b1; b1 = t + b1 * n
}
}
class NG8(
var a12: Int, var a1: Int, var a2: Int, var a: Int,
var b12: Int, var b1: Int, var b2: Int, var b: Int
) : MatrixNG() {
private var t = 0
private var ab = 0.0
private var a1b1 = 0.0
private var a2b2 = 0.0
private var a12b12 = 0.0
fun chooseCFN() = if (abs(a1b1 - ab) > abs(a2b2-ab)) 0 else 1
override fun needTerm(): Boolean {
if (b1 == 0 && b == 0 && b2 == 0 && b12 == 0) return false
if (b == 0) {
cfn = if (b2 == 0) 0 else 1
return true
}
else ab = a.toDouble() / b
if (b2 == 0) {
cfn = 1
return true
}
else a2b2 = a2.toDouble() / b2
if (b1 == 0) {
cfn = 0
return true
}
else a1b1 = a1.toDouble() / b1
if (b12 == 0) {
cfn = chooseCFN()
return true
}
else a12b12 = a12.toDouble() / b12
thisTerm = ab.toInt()
if (thisTerm == a1b1.toInt() && thisTerm == a2b2.toInt() &&
thisTerm == a12b12.toInt()) {
t = a; a = b; b = t - b * thisTerm
t = a1; a1 = b1; b1 = t - b1 * thisTerm
t = a2; a2 = b2; b2 = t - b2 * thisTerm
t = a12; a12 = b12; b12 = t - b12 * thisTerm
haveTerm = true
return false
}
cfn = chooseCFN()
return true
}
override fun consumeTerm() {
if (cfn == 0) {
a = a1; a2 = a12
b = b1; b2 = b12
}
else {
a = a2; a1 = a12
b = b2; b1 = b12
}
}
override fun consumeTerm(n: Int) {
if (cfn == 0) {
t = a; a = a1; a1 = t + a1 * n
t = a2; a2 = a12; a12 = t + a12 * n
t = b; b = b1; b1 = t + b1 * n
t = b2; b2 = b12; b12 = t + b12 * n
}
else {
t = a; a = a2; a2 = t + a2 * n
t = a1; a1 = a12; a12 = t + a12 * n
t = b; b = b2; b2 = t + b2 * n
t = b1; b1 = b12; b12 = t + b12 * n
}
}
}
interface ContinuedFraction {
fun nextTerm(): Int
fun moreTerms(): Boolean
}
class R2cf(var n1: Int, var n2: Int) : ContinuedFraction {
override fun nextTerm(): Int {
val thisTerm = n1 /n2
val t2 = n2
n2 = n1 - thisTerm * n2
n1 = t2
return thisTerm
}
override fun moreTerms() = abs(n2) > 0
}
class NG : ContinuedFraction {
val ng: MatrixNG
val n: List<ContinuedFraction>
constructor(ng: NG4, n1: ContinuedFraction) {
this.ng = ng
n = listOf(n1)
}
constructor(ng: NG8, n1: ContinuedFraction, n2: ContinuedFraction) {
this.ng = ng
n = listOf(n1, n2)
}
override fun nextTerm(): Int {
ng.haveTerm = false
return ng.thisTerm
}
override fun moreTerms(): Boolean {
while (ng.needTerm()) {
if (n[ng.cfn].moreTerms())
ng.consumeTerm(n[ng.cfn].nextTerm())
else
ng.consumeTerm()
}
return ng.haveTerm
}
}
fun test(desc: String, vararg cfs: ContinuedFraction) {
println("TESTING -> $desc")
for (cf in cfs) {
while (cf.moreTerms()) print ("${cf.nextTerm()} ")
println()
}
println()
}
fun main(args: Array<String>) {
val a = NG8(0, 1, 1, 0, 0, 0, 0, 1)
val n2 = R2cf(22, 7)
val n1 = R2cf(1, 2)
val a3 = NG4(2, 1, 0, 2)
val n3 = R2cf(22, 7)
test("[3;7] + [0;2]", NG(a, n1, n2), NG(a3, n3))
val b = NG8(1, 0, 0, 0, 0, 0, 0, 1)
val b1 = R2cf(13, 11)
val b2 = R2cf(22, 7)
test("[1;5,2] * [3;7]", NG(b, b1, b2), R2cf(286, 77))
val c = NG8(0, 1, -1, 0, 0, 0, 0, 1)
val c1 = R2cf(13, 11)
val c2 = R2cf(22, 7)
test("[1;5,2] - [3;7]", NG(c, c1, c2), R2cf(-151, 77))
val d = NG8(0, 1, 0, 0, 0, 0, 1, 0)
val d1 = R2cf(22 * 22, 7 * 7)
val d2 = R2cf(22,7)
test("Divide [] by [3;7]", NG(d, d1, d2))
val na = NG8(0, 1, 1, 0, 0, 0, 0, 1)
val a1 = R2cf(2, 7)
val a2 = R2cf(13, 11)
val aa = NG(na, a1, a2)
val nb = NG8(0, 1, -1, 0, 0, 0, 0, 1)
val b3 = R2cf(2, 7)
val b4 = R2cf(13, 11)
val bb = NG(nb, b3, b4)
val nc = NG8(1, 0, 0, 0, 0, 0, 0, 1)
val desc = "([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])"
test(desc, NG(nc, aa, bb), R2cf(-7797, 5929))
}
- Output:
TESTING -> [3;7] + [0;2] 3 1 1 1 4 3 1 1 1 4 TESTING -> [1;5,2] * [3;7] 3 1 2 2 3 1 2 2 TESTING -> [1;5,2] - [3;7] -1 -1 -24 -1 -2 -1 -1 -24 -1 -2 TESTING -> Divide [] by [3;7] 3 7 TESTING -> ([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) -1 -3 -5 -1 -2 -1 -26 -3 -1 -3 -5 -1 -2 -1 -26 -3
Nim
import strformat
####################################################################################################
type MatrixNG = ref object of RootObj
cfn: int
thisTerm: int
haveTerm: bool
method consumeTerm(m: MatrixNG) {.base.} =
raise newException(CatchableError, "Method without implementation override")
method consumeTerm(m: MatrixNG; n: int) {.base.} =
raise newException(CatchableError, "Method without implementation override")
method needTerm(m: MatrixNG): bool {.base.} =
raise newException(CatchableError, "Method without implementation override")
####################################################################################################
type NG4 = ref object of MatrixNG
a1, a, b1, b: int
proc newNG4(a1, a, b1, b: int): NG4 =
NG4(a1: a1, a: a, b1: b1, b: b)
method needTerm(ng: NG4): bool =
if ng.b1 == 0 and ng.b == 0: return false
if ng.b1 == 0 or ng.b == 0: return true
ng.thisTerm = ng.a div ng.b
if ng.thisTerm == ng.a1 div ng.b1:
ng.a -= ng.b * ng.thisTerm; swap ng.a, ng.b
ng.a1 -= ng.b1 * ng.thisTerm; swap ng.a1, ng.b1
ng.haveTerm = true
return false
return true
method consumeTerm(ng: NG4) =
ng.a = ng.a1
ng.b = ng.b1
method consumeTerm(ng: NG4; n: int) =
ng.a += ng.a1 * n; swap ng.a, ng.a1
ng.b += ng.b1 * n; swap ng.b, ng.b1
####################################################################################################
type NG8 = ref object of MatrixNG
a12, a1, a2, a: int
b12, b1, b2, b: int
proc newNG8(a12, a1, a2, a, b12, b1, b2, b: int): NG8 =
NG8(a12: a12, a1: a1, a2: a2, a: a, b12: b12, b1: b1, b2: b2, b: b)
method needTerm(ng: NG8): bool =
if ng.b1 == 0 and ng.b == 0 and ng.b2 == 0 and ng.b12 == 0: return false
if ng.b == 0:
ng.cfn = ord(ng.b2 != 0)
return true
if ng.b2 == 0:
ng. cfn = 1
return true
if ng.b1 == 0:
ng.cfn = 0
return true
let
ab = ng.a / ng.b
a1b1 = ng.a1 / ng.b1
a2b2 = ng.a2 / ng.b2
if ng.b12 == 0:
ng.cfn = ord(abs(a1b1 - ab) <= abs(a2b2 - ab))
return true
ng.thisTerm = int(ab)
if ng.thisTerm == int(a1b1) and ng.thisTerm == int(a2b2) and ng.thisTerm == ng.a12 div ng.b12:
ng.a -= ng.b * ng.thisTerm; swap ng.a, ng.b
ng.a1 -= ng.b1 * ng.thisTerm; swap ng.a1, ng.b1
ng.a2 -= ng.b2 * ng.thisTerm; swap ng.a2, ng.b2
ng.a12 -= ng.b12 * ng.thisTerm; swap ng.a12, ng.b12
ng.haveTerm = true
return false
ng.cfn = ord(abs(a1b1 - ab) <= abs(a2b2 - ab))
result = true
method consumeTerm(ng: NG8) =
if ng.cfn == 0:
ng.a = ng.a1
ng.a2 = ng.a12
ng.b = ng.b1
ng.b2 = ng.b12
else:
ng.a = ng.a2
ng.a1 = ng.a12
ng.b = ng.b2
ng.b1 = ng.b12
method consumeTerm(ng: NG8; n: int) =
if ng.cfn == 0:
ng.a += ng.a1 * n; swap ng.a, ng.a1
ng.a2 += ng.a12 * n; swap ng.a2, ng.a12
ng.b += ng.b1 * n; swap ng.b, ng.b1
ng.b2 += ng.b12 * n; swap ng.b2, ng.b12
else:
ng.a += ng.a2 * n; swap ng.a, ng.a2
ng.a1 += ng.a12 * n; swap ng.a1, ng.a12
ng.b += ng.b2 * n; swap ng.b, ng.b2
ng.b1 += ng.b12 * n; swap ng.b1, ng.b12
####################################################################################################
type ContinuedFraction = ref object of RootObj
method nextTerm(cf: ContinuedFraction): int {.base.} =
raise newException(CatchableError, "Method without implementation override")
method moreTerms(cf: ContinuedFraction): bool {.base.} =
raise newException(CatchableError, "Method without implementation override")
####################################################################################################
type R2Cf = ref object of ContinuedFraction
n1, n2: int
proc newR2Cf(n1, n2: int): R2Cf =
R2Cf(n1: n1, n2: n2)
method nextTerm(x: R2Cf): int =
result = x.n1 div x.n2
x.n1 -= result * x.n2
swap x.n1, x.n2
method moreTerms(x: R2Cf): bool =
abs(x.n2) > 0
####################################################################################################
type NG = ref object of ContinuedFraction
ng: MatrixNG
n: seq[ContinuedFraction]
proc newNG(ng: NG4; n1: ContinuedFraction): NG =
NG(ng: ng, n: @[n1])
proc newNG(ng: NG8; n1, n2: ContinuedFraction): NG =
NG(ng: ng, n: @[n1, n2])
method nextTerm(x: NG): int =
x.ng.haveTerm = false
result = x.ng.thisTerm
method moreTerms(x: NG): bool =
while x.ng.needTerm():
if x.n[x.ng.cfn].moreTerms():
x.ng.consumeTerm(x.n[x.ng.cfn].nextTerm())
else:
x.ng.consumeTerm()
result = x.ng.haveTerm
#———————————————————————————————————————————————————————————————————————————————————————————————————
when isMainModule:
proc test(desc: string; cfs: varargs[ContinuedFraction]) =
echo &"TESTING → {desc}"
for cf in cfs:
while cf.moreTerms(): stdout.write &"{cf.nextTerm()} "
echo()
echo()
let
a = newNG8(0, 1, 1, 0, 0, 0, 0, 1)
n2 = newR2Cf(22, 7)
n1 = newR2Cf(1, 2)
a3 = newNG4(2, 1, 0, 2)
n3 = newR2cf(22, 7)
test("[3;7] + [0;2]", newNG(a, n1, n2), newNG(a3, n3))
let
b = newNG8(1, 0, 0, 0, 0, 0, 0, 1)
b1 = newR2cf(13, 11)
b2 = newR2cf(22, 7)
test("[1;5,2] * [3;7]", newNG(b, b1, b2), newR2cf(286, 77))
let
c = newNG8(0, 1, -1, 0, 0, 0, 0, 1)
c1 = newR2cf(13, 11)
c2 = newR2cf(22, 7)
test("[1;5,2] - [3;7]", newNG(c, c1, c2), newR2cf(-151, 77))
let
d = newNG8(0, 1, 0, 0, 0, 0, 1, 0)
d1 = newR2cf(22 * 22, 7 * 7)
d2 = newR2cf(22,7)
test("Divide [] by [3;7]", newNG(d, d1, d2))
let
na = newNG8(0, 1, 1, 0, 0, 0, 0, 1)
a1 = newR2cf(2, 7)
a2 = newR2cf(13, 11)
aa = newNG(na, a1, a2)
nb = newNG8(0, 1, -1, 0, 0, 0, 0, 1)
b3 = newR2cf(2, 7)
b4 = newR2cf(13, 11)
bb = newNG(nb, b3, b4)
nc = newNG8(1, 0, 0, 0, 0, 0, 0, 1)
desc = "([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])"
test(desc, newNG(nc, aa, bb), newR2cf(-7797, 5929))
- Output:
TESTING → [3;7] + [0;2] 3 1 1 1 4 3 1 1 1 4 TESTING → [1;5,2] * [3;7] 3 1 2 2 3 1 2 2 TESTING → [1;5,2] - [3;7] -1 -1 -24 -1 -2 -1 -1 -24 -1 -2 TESTING → Divide [] by [3;7] 3 7 TESTING → ([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) -1 -3 -5 -1 -2 -1 -26 -3 -1 -3 -5 -1 -2 -1 -26 -3
Phix
(self-contained)
class full_matrix -- -- Used by apply_full_matrix() -- Note that each instance of full_matrix should be discarded after use. -- integer a12, a1, a2, a, b12, b1, b2, b function need_term() if b12==0 or b1==0 or b2==0 or b==0 then return true end if atom ab = a/b return ab!=a1/b1 or ab!=a1/b2 or ab!=a12/b12 end function function which_term() -- returns true for cf1, false for cf2 if b==0 and b2==0 then return true end if if b==0 or b2==0 then return false end if if b1==0 then return true end if atom ab = a/b return abs(a1/b1-ab) > abs(a2/b2-ab) end function function next_term() integer t = floor(a/b) sequence newas = {b12,b1,b2,b}, newbs = {a12-b12*t,a1-b1*t,a2-b2*t,a-b*t} {a12,a1,a2,a} = newas {b12,b1,b2,b} = newbs return t end function procedure in_term(bool is_cf1, object t={}) if integer(t) then sequence newas = iff(is_cf1?{a2+a12*t, a+a1*t, a12, a1} :{a1+a12*t, a12, a+a2*t, a2}), newbs = iff(is_cf1?{b2+b12*t, b+b1*t, b12, b1} :{b1+b12*t, b12, b+b2*t, b2}) {a12, a1, a2, a} = newas {b12, b1, b2, b} = newbs elsif is_cf1 then {a2, a, b2, b} = {a12, a1, b12, b1} else {a1, a, b1, b} = {a12, a2, b12, b2} end if end procedure function done() return b12==0 and b1==0 and b2==0 and b==0 end function end class function apply_full_matrix(sequence ctrl, cf1, cf2) -- -- If ctrl is {a12, a1, a2, a, -- b12, b1, b2, b} -- -- Then the result of apply_full_matrix(ctrl,cf1,cf2) would be -- -- (a12*cf1*cf2 + a1*cf1 + a2*cf2 + a) -- ----------------------------------- -- (b12*cf1*cf2 + b1*cf1 + b2*cf2 + b) -- -- For instance: -- { 0, 1, 1, 0, calculates cf1 + cf2 -- 0, 0, 0, 1} (divided by 1) -- -- { 0, 1,-1, 0, calculates cf1 - cf2 -- 0, 0, 0, 1} (divided by 1) -- -- { 1, 0, 0, 0, calculates cf1 * cf2 -- 0, 0, 0, 1} (divided by 1) -- -- { 0, 1, 0, 0, calculates cf1 -- 0, 0, 1, 0} divided by cf2 -- full_matrix fm = new(ctrl) sequence res = {} integer l1 = length(cf1), dx1=1, l2 = length(cf2), dx2=1 while true do if fm.need_term() then bool is_cf1 = fm.which_term() object t = {} if is_cf1 then if dx1<=l1 then t = cf1[dx1] dx1 += 1 end if else if dx2<=l2 then t = cf2[dx2] dx2 += 1 end if end if fm.in_term(is_cf1,t) else res &= fm.next_term() end if if fm.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 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 constant fmAdd = { 0, 1, 1, 0, 0, 0, 0, 1}, fmSub = { 0, 1,-1, 0, 0, 0, 0, 1}, fmMul = { 1, 0, 0, 0, 0, 0, 0, 1}, fmDiv = { 0, 1, 0, 0, 0, 0, 1, 0}, tests = {{"+",{22, 7},{ 1,2},fmAdd,22/7+1/2}, {"-",{13,11},{22,7},fmSub,13/11-22/7}, {"*",{13,11},{22,7},fmMul,13/11*22/7}, {"/",{22*22,7*7},{22,7},fmDiv,22/7}} for i=1 to length(tests) do {string op, sequence rat1, sequence rat2, sequence m, atom eres2} = tests[i] sequence cf1 = r2cf(rat1), cf2 = r2cf(rat2), cfr = apply_full_matrix(m,cf1,cf2) string bop = sprintf("%s %s %s",{cf2s(cf1),op,cf2s(cf2)}) printf(1,"%s is %s -> %s (est %g)\n",{bop,cf2s(cfr),cf2r(cfr),eres2}) end for
- Output:
[3;7] + [0;2] is [3;1,1,1,4] -> 3.642857142857143 (51/14) (est 3.64286) [1;5,2] - [3;7] is [-2;25,1,2] -> -1.961038961038961 (-151/77) (est -1.96104) [1;5,2] * [3;7] is [3;1,2,2] -> 3.714285714285714 (26/7) (est 3.71429) [9;1,7,6] / [3;7] is [3;7] -> 3.142857142857143 (22/7) (est 3.14286)
Raku
(formerly Perl 6)
The NG2 object can work with infinitely long continued fractions, it does lazy evaluation. By default, it is limited to returning the first 30 terms. Pass in a limit value if you want something other than default.
class NG2 {
has ( $!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b );
# Public methods
method operator($!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b ) { self }
method apply(@cf1, @cf2, :$limit = 30) {
my @cfs = [@cf1], [@cf2];
gather {
while @cfs[0] or @cfs[1] {
my $term;
(take $term if $term = self!extract) unless self!needterm;
my $from = self!from;
$from = @cfs[$from] ?? $from !! $from +^ 1;
self!inject($from, @cfs[$from].shift);
}
take self!drain while $!b;
}[ ^$limit ].grep: *.defined;
}
# Private methods
method !inject ($n, $t) {
multi sub xform(0, $t, $x12, $x1, $x2, $x) { $x2 + $x12 * $t, $x + $x1 * $t, $x12, $x1 }
multi sub xform(1, $t, $x12, $x1, $x2, $x) { $x1 + $x12 * $t, $x12, $x + $x2 * $t, $x2 }
( $!a12, $!a1, $!a2, $!a ) = xform($n, $t, $!a12, $!a1, $!a2, $!a );
( $!b12, $!b1, $!b2, $!b ) = xform($n, $t, $!b12, $!b1, $!b2, $!b );
}
method !extract {
my $t = $!a div $!b;
( $!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b ) =
$!b12, $!b1, $!b2, $!b,
$!a12 - $!b12 * $t,
$!a1 - $!b1 * $t,
$!a2 - $!b2 * $t,
$!a - $!b * $t;
$t;
}
method !from {
return $!b == $!b2 == 0 ?? 0 !!
$!b == 0 || $!b2 == 0 ?? 1 !!
abs($!a1*$!b*$!b2 - $!a*$!b1*$!b2) > abs($!a2*$!b*$!b1 - $!a*$!b1*$!b2) ?? 0 !! 1;
}
method !needterm {
so !([&&] $!b12, $!b1, $!b2, $!b) or $!a/$!b != $!a1/$!b1 != $!a2/$!b2 != $!a12/$!b1;
}
method !noterms($which) {
$which ?? (($!a1, $!a, $!b1, $!b ) = $!a12, $!a2, $!b12, $!b2)
!! (($!a2, $!a, $!b2, $!b ) = $!a12, $!a1, $!b12, $!b1);
}
method !drain {
self!noterms(self!from) if self!needterm;
self!extract;
}
}
sub r2cf(Rat $x is copy) { # Rational to continued fraction
gather loop {
$x -= take $x.floor;
last unless $x;
$x = 1 / $x;
}
}
sub cf2r(@a) { # continued fraction to Rational
my $x = @a[* - 1].FatRat; # Use FatRats for arbitrary precision
$x = @a[$_- 1] + 1 / $x for reverse 1 ..^ @a;
$x
}
# format continued fraction for pretty printing
sub ppcf(@cf) { "[{ @cf.join(',').subst(',',';') }]" }
# format Rational for pretty printing. Use FatRats for arbitrary precision
sub pprat($a) { $a.FatRat.denominator == 1 ?? $a !! $a.FatRat.nude.join('/') }
my %ops = ( # convenience hash of NG matrix operators
'+' => (0,1,1,0,0,0,0,1),
'-' => (0,1,-1,0,0,0,0,1),
'*' => (1,0,0,0,0,0,0,1),
'/' => (0,1,0,0,0,0,1,0)
);
sub test_NG2 ($rat1, $op, $rat2) {
my @cf1 = $rat1.&r2cf;
my @cf2 = $rat2.&r2cf;
my @result = NG2.new.operator(|%ops{$op}).apply( @cf1, @cf2 );
say "{$rat1.&pprat} $op {$rat2.&pprat} => {@cf1.&ppcf} $op ",
"{@cf2.&ppcf} = {@result.&ppcf} => {@result.&cf2r.&pprat}\n";
}
# Testing
test_NG2(|$_) for
[ 22/7, '+', 1/2 ],
[ 23/11, '*', 22/7 ],
[ 13/11, '-', 22/7 ],
[ 484/49, '/', 22/7 ];
# Sometimes you may want to limit the terms in the continued fraction to something other than default.
# Here a lazy infinite continued fraction for √2, then multiply it by itself. We'll limit the result
# to 6 terms for brevity’s' sake. We'll then convert that continued fraction back to an arbitrary precision
# FatRat Rational number. (Raku stores FatRats internally as a ratio of two arbitrarily long integers.
# We need to exercise a little caution because they can eat up all of your memory if allowed to grow unchecked,
# hence the limit of 6 terms in continued fraction.) We'll then convert that number to a normal precision
# Rat, which is accurate to the nearest 1 / 2^64,
say "√2 expressed as a continued fraction, then squared: ";
my @root2 = lazy flat 1, 2 xx *;
my @result = NG2.new.operator(|%ops{'*'}).apply( @root2, @root2, limit => 6 );
say @root2.&ppcf, "² = \n";
say @result.&ppcf;
say "\nConverted back to an arbitrary (ludicrous) precision Rational: ";
say @result.&cf2r.nude.join(" /\n");
say "\nCoerced to a standard precision Rational: ", @result.&cf2r.Num.Rat;
- Output:
22/7 + 1/2 => [3;7] + [0;2] = [3;1,1,1,4] => 51/14 23/11 * 22/7 => [2;11] * [3;7] = [6;1,1,3] => 46/7 13/11 - 22/7 => [1;5,2] - [3;7] = [-2;25,1,2] => -151/77 484/49 / 22/7 => [9;1,7,6] / [3;7] = [3;7] => 22/7 √2 expressed as a continued fraction, then squared: [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,...]² = [1;1,-58451683124983302025,-1927184886226364356176,-65467555105469489418600,-2223969688699736275876224] Converted back to an arbitrary (ludicrous) precision Rational: 32802382178012409621354320392819425499699206367450594986122623570838188983519955166754002 / 16401191089006204810536863200564985394427741343927508600629139291039556821665755787817601 Coerced to a standard precision Rational: 2
Tcl
This uses the Generator class, R2CF class and printcf procedure from the r2cf task.
oo::class create NG2 {
variable a b a1 b1 a2 b2 a12 b12 cf1 cf2
superclass Generator
constructor {args} {
lassign $args a12 a1 a2 a b12 b1 b2 b
next
}
method operands {N1 N2} {
set cf1 $N1
set cf2 $N2
return [self]
}
method Ingress1 t {
lassign [list [expr {$a2+$a12*$t}] [expr {$a+$a1*$t}] $a12 $a1 \
[expr {$b2+$b12*$t}] [expr {$b+$b1*$t}] $b12 $b1] \
a12 a1 a2 a b12 b1 b2 b
}
method Exhaust1 {} {
lassign [list $a12 $a1 $a12 $a1 $b12 $b1 $b12 $b1] \
a12 a1 a2 a b12 b1 b2 b
}
method Ingress2 t {
lassign [list [expr {$a1+$a12*$t}] $a12 [expr {$a+$a2*$t}] $a2 \
[expr {$b1+$b12*$t}] $b12 [expr {$b+$b2*$t}] $b2] \
a12 a1 a2 a b12 b1 b2 b
}
method Exhaust2 {} {
lassign [list $a12 $a12 $a2 $a2 $b12 $b12 $b2 $b2] \
a12 a1 a2 a b12 b1 b2 b
}
method Egress {} {
set t [expr {$a/$b}]
lassign [list $b12 $b1 $b2 $b \
[expr {$a12 - $b12*$t}] [expr {$a1 - $b1*$t}] \
[expr {$a2 - $b2*$t}] [expr {$a - $b*$t}]] \
a12 a1 a2 a b12 b1 b2 b
return $t
}
method DoIngress1 {} {
try {tailcall my Ingress1 [$cf1]} on break {} {}
oo::objdefine [self] forward DoIngress1 my Exhaust1
set cf1 ""
tailcall my Exhaust1
}
method DoIngress2 {} {
try {tailcall my Ingress2 [$cf2]} on break {} {}
oo::objdefine [self] forward DoIngress2 my Exhaust2
set cf2 ""
tailcall my Exhaust2
}
method Ingress {} {
if {$b==0} {
if {$b2 == 0} {
tailcall my DoIngress1
} else {
tailcall my DoIngress2
}
}
if {!$b2} {
tailcall my DoIngress2
}
if {!$b1} {
tailcall my DoIngress1
}
if {[my FirstSource?]} {
tailcall my DoIngress1
} else {
tailcall my DoIngress2
}
}
method FirstSource? {} {
expr {abs($a1*$b*$b2 - $a*$b1*$b2) > abs($a2*$b*$b1 - $a*$b1*$b2)}
}
method NeedTerm? {} {
expr {
($b*$b1*$b2*$b12==0) ||
!($a/$b == $a1/$b1 && $a/$b == $a2/$b2 && $a/$b == $a12/$b12)
}
}
method Done? {} {
expr {$b==0 && $b1==0 && $b2==0 && $b12==0}
}
method Produce {} {
# Until we've drained both continued fractions...
while {$cf1 ne "" || $cf2 ne ""} {
if {[my NeedTerm?]} {
my Ingress
} else {
yield [my Egress]
}
}
# Drain our internal state
while {![my Done?]} {
yield [my Egress]
}
}
}
Demonstrating:
set op [[NG2 new 0 1 1 0 0 0 0 1] operands [R2CF new 1/2] [R2CF new 22/7]]
printcf "\[3;7\] + \[0;2\]" $op
set op [[NG2 new 1 0 0 0 0 0 0 1] operands [R2CF new 13/11] [R2CF new 22/7]]
printcf "\[1:5,2\] * \[3;7\]" $op
set op [[NG2 new 0 1 -1 0 0 0 0 1] operands [R2CF new 13/11] [R2CF new 22/7]]
printcf "\[1:5,2\] - \[3;7\]" $op
set op [[NG2 new 0 1 0 0 0 0 1 0] operands [R2CF new 484/49] [R2CF new 22/7]]
printcf "div test" $op
set op1 [[NG2 new 0 1 1 0 0 0 0 1] operands [R2CF new 2/7] [R2CF new 13/11]]
set op2 [[NG2 new 0 1 -1 0 0 0 0 1] operands [R2CF new 2/7] [R2CF new 13/11]]
set op3 [[NG2 new 1 0 0 0 0 0 0 1] operands $op1 $op2]
printcf "layered test" $op3
- Output:
[3;7] + [0;2] -> 3,1,1,1,4 [1:5,2] * [3;7]-> 3,1,2,2 [1:5,2] - [3;7]-> -2,25,1,2 div test -> 3,7 layered test -> -2,1,2,5,1,2,1,26,3
Wren
class MatrixNG {
construct new() {
_cfn = 0
_thisTerm = 0
_haveTerm = false
}
cfn { _cfn }
cfn=(v) { _cfn = v }
thisTerm { _thisTerm }
thisTerm=(v) { _thisTerm = v }
haveTerm { _haveTerm }
haveTerm=(v) { _haveTerm = v }
consumeTerm() {}
consumeTerm(n) {}
needTerm() {}
}
class NG4 is MatrixNG {
construct new(a1, a, b1, b) {
super()
_a1 = a1
_a = a
_b1 = b1
_b = b
_t = 0
}
needTerm() {
if (_b1 == 0 && _b == 0) return false
if (_b1 == 0 || _b == 0) return true
thisTerm = (_a / _b).truncate
if (thisTerm == (_a1 / _b1).truncate) {
_t = _a
_a = _b
_b = _t - _b * thisTerm
_t = _a1
_a1 = _b1
_b1 = _t - _b1 * thisTerm
haveTerm = true
return false
}
return true
}
consumeTerm() {
_a = _a1
_b = _b1
}
consumeTerm(n) {
_t = _a
_a = _a1
_a1 = _t + _a1 * n
_t = _b
_b = _b1
_b1 = _t + _b1 * n
}
}
class NG8 is MatrixNG {
construct new(a12, a1, a2, a, b12, b1, b2, b) {
super()
_a12 = a12
_a1 = a1
_a2 = a2
_a = a
_b12 = b12
_b1 = b1
_b2 = b2
_b = b
_t = 0
_ab = 0
_a1b1 = 0
_a2b2 = 0
_a12b12 = 0
}
chooseCFN() { ((_a1b1 - _ab).abs > (_a2b2 - _ab).abs) ? 0 : 1 }
needTerm() {
if (_b1 == 0 && _b == 0 && _b2 == 0 && _b12 == 0) return false
if (_b == 0) {
cfn = (_b2 == 0) ? 0 : 1
return true
} else _ab = _a/_b
if (_b2 == 0) {
cfn = 1
return true
} else _a2b2 = _a2/_b2
if (_b1 == 0) {
cfn = 0
return true
} else _a1b1 = _a1/_b1
if (_b12 == 0) {
cfn = chooseCFN()
return true
} else _a12b12 = _a12/_b12
thisTerm = _ab.truncate
if (thisTerm == _a1b1.truncate && thisTerm == _a2b2.truncate &&
thisTerm == _a12b12.truncate) {
_t = _a
_a = _b
_b = _t - _b * thisTerm
_t = _a1
_a1 = _b1
_b1 = _t - _b1 * thisTerm
_t = _a2
_a2 = _b2
_b2 = _t - _b2 * thisTerm
_t = _a12
_a12 = _b12
_b12 = _t - _b12 * thisTerm
haveTerm = true
return false
}
cfn = chooseCFN()
return true
}
consumeTerm() {
if (cfn == 0) {
_a = _a1
_a2 = _a12
_b = _b1
_b2 = _b12
} else {
_a = _a2
_a1 = _a12
_b = _b2
_b1 = _b12
}
}
consumeTerm(n) {
if (cfn == 0) {
_t = _a
_a = _a1
_a1 = _t + _a1 * n
_t = _a2
_a2 = _a12
_a12 = _t + _a12 * n
_t = _b
_b = _b1
_b1 = _t + _b1 * n
_t = _b2
_b2 = _b12
_b12 = _t + _b12 * n
} else {
_t = _a
_a = _a2
_a2 = _t + _a2 * n
_t = _a1
_a1 = _a12
_a12 = _t + _a12 * n
_t = _b
_b = _b2
_b2 = _t + _b2 * n
_t = _b1
_b1 = _b12
_b12 = _t + _b12 * n
}
}
}
class ContinuedFraction {
nextTerm() {}
moreTerms() {}
}
class R2cf is ContinuedFraction {
construct new(n1, n2) {
_n1 = n1
_n2 = n2
}
nextTerm() {
var thisTerm = (_n1/_n2).truncate
var t2 = _n2
_n2 = _n1 - thisTerm * _n2
_n1 = t2
return thisTerm
}
moreTerms() { _n2.abs > 0 }
}
class NG is ContinuedFraction {
construct new(ng, n1) {
_ng = ng
_n = [n1]
}
construct new(ng, n1, n2) {
_ng = ng
_n = [n1, n2]
}
nextTerm() {
_ng.haveTerm = false
return _ng.thisTerm
}
moreTerms() {
while (_ng.needTerm()) {
if (_n[_ng.cfn].moreTerms()) {
_ng.consumeTerm(_n[_ng.cfn].nextTerm())
} else {
_ng.consumeTerm()
}
}
return _ng.haveTerm
}
}
var test = Fn.new { |desc, cfs|
System.print("TESTING -> %(desc)")
for (cf in cfs) {
while (cf.moreTerms()) System.write("%(cf.nextTerm()) ")
System.print()
}
System.print()
}
var a = NG8.new(0, 1, 1, 0, 0, 0, 0, 1)
var n2 = R2cf.new(22, 7)
var n1 = R2cf.new(1, 2)
var a3 = NG4.new(2, 1, 0, 2)
var n3 = R2cf.new(22, 7)
test.call("[3;7] + [0;2]", [NG.new(a, n1, n2), NG.new(a3, n3)])
var b = NG8.new(1, 0, 0, 0, 0, 0, 0, 1)
var b1 = R2cf.new(13, 11)
var b2 = R2cf.new(22, 7)
test.call("[1;5,2] * [3;7]", [NG.new(b, b1, b2), R2cf.new(286, 77)])
var c = NG8.new(0, 1, -1, 0, 0, 0, 0, 1)
var c1 = R2cf.new(13, 11)
var c2 = R2cf.new(22, 7)
test.call("[1;5,2] - [3;7]", [NG.new(c, c1, c2), R2cf.new(-151, 77)])
var d = NG8.new(0, 1, 0, 0, 0, 0, 1, 0)
var d1 = R2cf.new(22 * 22, 7 * 7)
var d2 = R2cf.new(22,7)
test.call("Divide [] by [3;7]", [NG.new(d, d1, d2)])
var na = NG8.new(0, 1, 1, 0, 0, 0, 0, 1)
var a1 = R2cf.new(2, 7)
var a2 = R2cf.new(13, 11)
var aa = NG.new(na, a1, a2)
var nb = NG8.new(0, 1, -1, 0, 0, 0, 0, 1)
var b3 = R2cf.new(2, 7)
var b4 = R2cf.new(13, 11)
var bb = NG.new(nb, b3, b4)
var nc = NG8.new(1, 0, 0, 0, 0, 0, 0, 1)
var desc = "([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])"
test.call(desc, [NG.new(nc, aa, bb), R2cf.new(-7797, 5929)])- Output:
TESTING -> [3;7] + [0;2] 3 1 1 1 4 3 1 1 1 4 TESTING -> [1;5,2] * [3;7] 3 1 2 2 3 1 2 2 TESTING -> [1;5,2] - [3;7] -1 -1 -24 -1 -2 -1 -1 -24 -1 -2 TESTING -> Divide [] by [3;7] 3 7 TESTING -> ([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) -1 -3 -5 -1 -2 -1 -26 -3 -1 -3 -5 -1 -2 -1 -26 -3