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
(Margin note: This program is a bug-fix of an ATS program on which I based Python code. It does not really matter, however, which program came first. So I am calling this a translation from Python.)
The following program uses GNU C extensions: 128-bit integers, and integer operations that detect overflow. I add the 128-bit integers as a new g0int(int128knd) type. The overflow-detecting operations I call +!, -!, and *!. There are no multiple-precision integers or rationals. Rather than detect numbers "getting too large", I catch integer overflows. I use the most negative 128-bit value to represent "no finite term".
I have not included any code to weed out large terms that show up where "infinities" belong.
(Margin note: Sometimes infinite sequences get truncated due to overflow, though I have not seen this happen very near the beginning of a continued fraction. You can give a "maximum number of terms to print" argument to this program, to see the phenomenon for yourself. It happens relatively quickly with 128-bit integers.)
(*------------------------------------------------------------------*)
#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"
%{#
#include <stdint.h>
#include <limits.h>
#include <float.h>
#include <math.h>
%}
#define NIL list_nil ()
#define :: list_cons
exception gint_overflow of ()
(*------------------------------------------------------------------*)
extern fn {tk : tkind}
g0int_neginf :
() -<> g0int tk
extern fn {tk : tkind}
g0int_add_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_sub_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
extern fn {tk : tkind}
g0int_mul_overflow_exn :
(g0int tk, g0int tk) -< !exn > g0int tk
infixl ( + ) +!
infixl ( - ) -!
infixl ( * ) *!
overload neginf with g0int_neginf
overload +! with g0int_add_overflow_exn
overload -! with g0int_sub_overflow_exn
overload *! with g0int_mul_overflow_exn
(*------------------------------------------------------------------*)
(* 128-bit integers. *)
%{#
/* The most negative int128 will be treated as "neginf" or "negative
infinity". For our purposes the sign will not matter, though. */
#define neginf_int128() (((__int128) 1) << 127)
#define neg_c(x) (-(x))
#define add_c(x, y) ((x) + (y))
#define sub_c(x, y) ((x) - (y))
#define mul_c(x, y) ((x) * (y))
#define div_c(x, y) ((x) / (y))
#define mod_c(x, y) ((x) % (y))
#define eq_c(x, y) (((x) == (y)) ? atsbool_true : atsbool_false)
#define neq_c(x, y) (((x) != (y)) ? atsbool_true : atsbool_false)
#define lt_c(x, y) (((x) < (y)) ? atsbool_true : atsbool_false)
#define lte_c(x, y) (((x) <= (y)) ? atsbool_true : atsbool_false)
#define gt_c(x, y) (((x) > (y)) ? atsbool_true : atsbool_false)
#define gte_c(x, y) (((x) >= (y)) ? atsbool_true : atsbool_false)
/* GNU extensions for detection of integer overflow. */
#define add_overflow(x, y, pz) \
(__builtin_add_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define sub_overflow(x, y, pz) \
(__builtin_sub_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
#define mul_overflow(x, y, pz) \
(__builtin_mul_overflow ((x), (y), (pz)) ? \
atsbool_true : atsbool_false)
%}
tkindef int128_kind = "__int128" (* A GNU extension. *)
stadef int128knd = int128_kind
typedef int128_0 = g0int int128knd
typedef int128_1 (i : int) = g1int (int128knd, i)
stadef int128 = int128_1 // 2nd-select
stadef int128 = int128_0 // 1st-select
stadef Int128 = [i : int] int128_1 i
extern fn g0int_neginf_int128 : () -<> int128 = "mac#neginf_int128"
extern fn g0int_neg_int128 : int128 -<> int128 = "mac#neg_c"
extern fn g0int_add_int128 : (int128, int128) -<> int128 = "mac#add_c"
extern fn g0int_sub_int128 : (int128, int128) -<> int128 = "mac#sub_c"
extern fn g0int_mul_int128 : (int128, int128) -<> int128 = "mac#mul_c"
extern fn g0int_div_int128 : (int128, int128) -<> int128 = "mac#div_c"
extern fn g0int_mod_int128 : (int128, int128) -<> int128 = "mac#mod_c"
extern fn g0int_eq_int128 : (int128, int128) -<> bool = "mac#eq_c"
extern fn g0int_neq_int128 : (int128, int128) -<> bool = "mac#neq_c"
extern fn g0int_lt_int128 : (int128, int128) -<> bool = "mac#lt_c"
extern fn g0int_lte_int128 : (int128, int128) -<> bool = "mac#lte_c"
extern fn g0int_gt_int128 : (int128, int128) -<> bool = "mac#gt_c"
extern fn g0int_gte_int128 : (int128, int128) -<> bool = "mac#gte_c"
implement g0int_neginf<int128knd> () = g0int_neginf_int128 ()
implement g0int_neg<int128knd> x = g0int_neg_int128 x
implement g0int_add<int128knd> (x, y) = g0int_add_int128 (x, y)
implement g0int_sub<int128knd> (x, y) = g0int_sub_int128 (x, y)
implement g0int_mul<int128knd> (x, y) = g0int_mul_int128 (x, y)
implement g0int_div<int128knd> (x, y) = g0int_div_int128 (x, y)
implement g0int_mod<int128knd> (x, y) = g0int_mod_int128 (x, y)
implement g0int_eq<int128knd> (x, y) = g0int_eq_int128 (x, y)
implement g0int_neq<int128knd> (x, y) = g0int_neq_int128 (x, y)
implement g0int_lt<int128knd> (x, y) = g0int_lt_int128 (x, y)
implement g0int_lte<int128knd> (x, y) = g0int_lte_int128 (x, y)
implement g0int_gt<int128knd> (x, y) = g0int_gt_int128 (x, y)
implement g0int_gte<int128knd> (x, y) = g0int_gte_int128 (x, y)
implement g0int2int<intknd,int128knd> i = $UN.cast i
implement g0int2int<int128knd,intknd> i = $UN.cast i
implement g0int2float<int128knd,ldblknd> i = $UN.cast i
implement g0int_iseqz<int128knd> x = (x = g0i2i 0)
implement g0int_isneqz<int128knd> x = (x <> g0i2i 0)
implement g0int_isltz<int128knd> x = (x < g0i2i 0)
implement g0int_isltez<int128knd> x = (x <= g0i2i 0)
implement g0int_isgtz<int128knd> x = (x > g0i2i 0)
implement g0int_isgtez<int128knd> x = (x >= g0i2i 0)
implement g0int_abs<int128knd> x = (if isltz x then ~x else x)
local
extern fn
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
(*------------------------------------------------------------------*)
(* The continued fraction for a rational number or integer. *)
typedef ratnum = @(integer, integer)
fn
r2cf (ratnum : ratnum) : cf =
cf_make
let
val ratnum_ref : ref ratnum = ref ratnum
in
lam () =<cloref1>
let
val @(n, d) = !ratnum_ref
in
if iseqz d then
neginf<integerknd> ()
else
let
val @(q, r) = n divrem d
in
!ratnum_ref := @(d, r);
q
end
end
end
fn
i2cf (intnum : integer) : cf =
r2cf @(intnum, g0i2i 1)
(*------------------------------------------------------------------*)
(* Application of a homographic function to a continued fraction. *)
typedef ng4 = @(integer, integer, integer, integer)
fn
apply_ng4 (ng4 : ng4, x : cf) : cf =
cf_make
let
val state : ref @(ng4, Size_t) = ref @(ng4, i2sz 0)
in
lam () =<cloref1>
let
fnx
loop (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
if (iseqz b1) * (iseqz b) then
neginf<integerknd> ()
else if (iseqz b1) + (iseqz b) then
absorb_term (a1, a, b1, b, i)
else
let
val @(q, r) = a divrem b
and @(q1, r1) = a1 divrem b1
in
if q1 <> q then
absorb_term (a1, a, b1, b, i)
else
begin
!state :=
@(@(b1, b, r1, r), i);
q
end
end
and
absorb_term (a1 : integer,
a : integer,
b1 : integer,
b : integer,
i : Size_t)
: integer =
let
val term = x[i]
in
if term <> neginf<integerknd> () then
loop (a + (a1 * term), a1,
b + (b1 * term), b1, succ i)
else
loop (a1, a1, b1, b1, succ i)
end
val @(@(a1, a, b1, b), i) = !state
in
loop (a1, a, b1, b, i)
end
end
(*------------------------------------------------------------------*)
(* Some basic operations involving only one continued fraction. *)
fn
cf_neg (x : cf) : cf =
apply_ng4 (@(g0i2i ~1, g0i2i 0, g0i2i 0, g0i2i 1), x)
fn
cf_add_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(d, n, g0i2i 0, d), x)
end
fn
ratnum_add_cf (x : ratnum, y : cf) : cf =
cf_add_ratnum (y, x)
fn
cf_add_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(y, g0i2i 1))
fn
integer_add_cf (x : integer, y : cf) : cf =
cf_add_ratnum (y, @(x, g0i2i 1))
fn
cf_sub_ratnum (x : cf, y : ratnum) : cf =
cf_add_ratnum (x, @(~(y.0), (y.1)))
fn
ratnum_sub_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(~d, n, g0i2i 0, d), y)
end
fn
cf_sub_integer (x : cf, y : integer) : cf =
cf_add_ratnum (x, @(~y, g0i2i 1))
fn
integer_sub_cf (x : integer, y : cf) : cf =
ratnum_sub_cf (@(x, g0i2i 1), y)
fn
cf_mul_ratnum (x : cf, y : ratnum) : cf =
let
val @(n, d) = y
in
apply_ng4 (@(n, g0i2i 0, g0i2i 0, d), x)
end
fn
ratnum_mul_cf (x : ratnum, y : cf) : cf =
cf_mul_ratnum (y, x)
fn
cf_mul_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(y, g0i2i 1))
fn
integer_mul_cf (x : integer, y : cf) : cf =
cf_mul_ratnum (y, @(x, g0i2i 1))
fn
cf_div_ratnum (x : cf, y : ratnum) : cf =
cf_mul_ratnum (x, @(y.1, y.0))
fn
ratnum_div_cf (x : ratnum, y : cf) : cf =
let
val @(n, d) = x
in
apply_ng4 (@(g0i2i 0, n, d, g0i2i 0), y)
end
fn
cf_div_integer (x : cf, y : integer) : cf =
cf_mul_ratnum (x, @(g0i2i 1, y))
fn
integer_div_cf (x : integer, y : cf) : cf =
ratnum_div_cf (@(x, g0i2i 1), y)
overload ~ with cf_neg
overload + with cf_add_ratnum
overload + with ratnum_add_cf
overload + with cf_add_integer
overload + with integer_add_cf
overload - with cf_sub_ratnum
overload - with ratnum_sub_cf
overload - with cf_sub_integer
overload - with integer_sub_cf
overload * with cf_mul_ratnum
overload * with ratnum_mul_cf
overload * with cf_mul_integer
overload * with integer_mul_cf
overload / with cf_div_ratnum
overload / with ratnum_div_cf
overload / with cf_div_integer
overload / with integer_div_cf
(*------------------------------------------------------------------*)
(* Application of a bihomographic function to a continued fraction. *)
typedef ng8 = @(integer, integer, integer, integer,
integer, integer, integer, integer)
fn
apply_ng8 (ng8 : ng8, x : cf, y : cf) : cf =
cf_make
let
typedef state = '(ng8,
Size_t -<cloref1> integer,
Size_t -<cloref1> integer,
Size_t,
Size_t)
fn xsource (i : Size_t) :<cloref1> integer = x[i]
fn ysource (i : Size_t) :<cloref1> integer = y[i]
fn neginf_source (i : Size_t) :<cloref1> integer =
neginf<integerknd> ()
val state_ref : ref state =
ref '(ng8, xsource, ysource, i2sz 0, i2sz 0)
in
lam () =<cloref1>
let
fnx
loop () : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
xsource, ysource, i, j) = !state_ref
in
if iseqz b12 * iseqz b1 * iseqz b2 * iseqz b then
neginf<integerknd> ()
else if iseqz b * iseqz b2 then
absorb_x_term ()
else if iseqz b + iseqz b2 then
absorb_y_term ()
else if iseqz b1 then
absorb_x_term ()
else
let
val @(q, r) = a divrem b
and @(q1, r1) = a1 divrem b1
and @(q2, r2) = a2 divrem b2
and @(q12, r12) =
(if isneqz b12 then
a12 divrem b12
else
@(neginf (), neginf ())) : @(integer, integer)
in
if (isneqz b12)
* (q = q1) * (q = q2) * (q = q12) then
begin (* Output a term. *)
!state_ref :=
'(@(b12, b1, b2, b, r12, r1, r2, r),
xsource, ysource, i, j);
q
end
else
let
(* Put numerators over a common denominator,
then compare the numerators. *)
val n = a *! b1 *! b2
and n1 = a1 *! b *! b2
and n2 = a2 *! b *! b1
in
if abs (n1 -! n) > abs (n2 -! n) then
absorb_x_term ()
else
absorb_y_term ()
end
end
end
and
absorb_x_term () : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
xsource, ysource, i, j) = !state_ref
val term = xsource i
in
if term <> neginf<integerknd> () then
begin
try
let
val a12 = a2 +! (a12 *! term)
and a1 = a +! (a1 *! term)
and a2 = a12
and a = a1
and b12 = b2 +! (b12 *! term)
and b1 = b +! (b1 *! term)
and b2 = b12
and b = b1
in
!state_ref :=
'(@(a12, a1, a2, a, b12, b1, b2, b),
xsource, ysource, succ i, j);
(* Be aware: this is not a tail recursion! *)
loop ()
end
with
| ~ gint_overflow () =>
begin
(* Replace the sources with ones that return no
terms. (You have to replace BOTH sources.) *)
!state_ref :=
'(@(a12, a1, a12, a1, b12, b1, b12, b1),
neginf_source, neginf_source, succ i, j);
loop ()
end
end
else
begin
!state_ref :=
'(@(a12, a1, a12, a1, b12, b1, b12, b1),
xsource, ysource, succ i, j);
loop ()
end
end
and
absorb_y_term () : integer =
let
val '(@(a12, a1, a2, a, b12, b1, b2, b),
xsource, ysource, i, j) = !state_ref
val term = ysource j
in
if term <> neginf<integerknd> () then
begin
try
let
val a12 = a1 +! (a12 *! term)
and a1 = a12
and a2 = a +! (a2 *! term)
and a = a2
and b12 = b1 +! (b12 *! term)
and b1 = b12
and b2 = b +! (b2 *! term)
and b = b2
in
!state_ref :=
'(@(a12, a1, a2, a, b12, b1, b2, b),
xsource, ysource, i, succ j);
(* Be aware: this is not a tail recursion! *)
loop ()
end
with
| ~ gint_overflow () =>
begin
(* Replace the sources with ones that return no
terms. (You have to replace BOTH sources.) *)
!state_ref :=
'(@(a12, a12, a2, a2, b12, b12, b2, b2),
neginf_source, neginf_source, i, succ j);
loop ()
end
end
else
begin
!state_ref :=
'(@(a12, a12, a2, a2, b12, b12, b2, b2),
xsource, ysource, i, succ j);
loop ()
end
end
in
loop ()
end
end
(*------------------------------------------------------------------*)
(* Some basic operations on two continued fractions. *)
fn
cf_add_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_sub_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i ~1, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_mul_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 1), x, y)
fn
cf_div_cf (x : cf, y : cf) : cf =
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 1, g0i2i 0), x, y)
overload + with cf_add_cf
overload - with cf_sub_cf
overload * with cf_mul_cf
overload / with cf_div_cf
(*------------------------------------------------------------------*)
typedef charptr = $extype"char *"
fn
show_with_note (expression : string,
cf : cf,
note : string)
: void =
if note = "" then
ignoret ($extfcall (int, "printf", "%19s => %s\n",
$UN.cast{charptr} expression,
$UN.cast{charptr} (cf2string cf)))
else
ignoret ($extfcall (int, "printf", "%19s => %-46s %s\n",
$UN.cast{charptr} expression,
$UN.cast{charptr} (cf2string cf),
$UN.cast{charptr} note))
fn
show_without_note (expression : string, cf : cf) : void =
show_with_note (expression, cf, "")
overload show with show_with_note
overload show with show_without_note
(*------------------------------------------------------------------*)
val golden_ratio = cf_make (lam () =<cloref1> (g0i2i 1) : integer)
val silver_ratio = cf_make (lam () =<cloref1> (g0i2i 2) : integer)
val sqrt2 = silver_ratio - g0i2i 1
val frac_13_11 = r2cf @(g0i2i 13, g0i2i 11)
val frac_22_7 = r2cf @(g0i2i 22, g0i2i 7)
val one = i2cf (g0i2i 1)
val two = i2cf (g0i2i 2)
val three = i2cf (g0i2i 3)
val four = i2cf (g0i2i 4)
implement
main (argc, argv) =
begin
if 2 <= argc then
let
val n = g0string2uint (argv_get_at (argv, 1)) : uint
in
(* Set the maximum number of terms to print. *)
!cf2string_default_max_terms := g1u2u (g1ofg0 n)
end;
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "1 + sqrt(2)");
show ("sqrt(2)", sqrt2);
show ("13/11", frac_13_11);
show ("22/7", frac_22_7);
show ("one", one);
show ("two", two);
show ("three", three);
show ("four", four);
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (@(g0i2i 0, g0i2i 1, g0i2i 0, g0i2i 0,
g0i2i 0, g0i2i 0, g0i2i 2, g0i2i 0),
silver_ratio, sqrt2),
"method 1");
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (@(g0i2i 1, g0i2i 0, g0i2i 0, g0i2i 1,
g0i2i 0, g0i2i 0, g0i2i 0, g0i2i 8),
silver_ratio, silver_ratio),
"method 2");
show ("(1 + 1/sqrt(2))/2", (g0i2i 1 + (one / sqrt2)) / two,
"method 3");
show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2);
show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2);
show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2);
show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2);
0
end
(*------------------------------------------------------------------*)- Output:
$ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW bivariate-continued-fraction-task-memoizing.dats -lgc && ./a.out
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 1 + sqrt(2)
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
one => [1]
two => [2]
three => [3]
four => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2;7530688524100]
sqrt(2) / sqrt(2) => [1]
C
You will need Boehm GC and the GNU Multiple Precision Library.
(Actually you can leave out the Boehm GC parts. The program leaks memory, but harmlessly. Also, you can leave out the C23 attribute specifiers.)
/*------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <gc/gc.h> /* Boehm GC. */
#include <gmp.h> /* GNU Multiple Precision. */
void *
my_malloc (size_t size)
{
void *p = GC_MALLOC (size);
if (p == NULL)
{
fprintf (stderr, "Memory allocation failed.\n");
exit (1);
}
return p;
}
void *
my_realloc (void *p, size_t size)
{
void *q = GC_REALLOC (p, size);
if (q == NULL)
{
fprintf (stderr, "Memory allocation failed.\n");
exit (1);
}
return q;
}
void
my_free (void *p)
{
GC_FREE (p);
}
/*------------------------------------------------------------------*/
/* Some helper functions. */
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
#define MAX(x, y) (((x) > (y)) ? (x) : (y))
char *
string_append1 (const char *s1)
{
size_t n1 = strlen (s1);
char *s = my_malloc ((n1 + 1) * sizeof (char));
s[n1] = '\0';
memcpy (s, s1, n1);
return s;
}
char *
string_append2 (const char *s1, const char *s2)
{
size_t n1 = strlen (s1);
size_t n2 = strlen (s2);
char *s = my_malloc ((n1 + n2 + 1) * sizeof (char));
s[n1 + n2] = '\0';
memcpy (s, s1, n1);
memcpy (s + n1, s2, n2);
return s;
}
char *
string_append3 (const char *s1, const char *s2, const char *s3)
{
size_t n1 = strlen (s1);
size_t n2 = strlen (s2);
size_t n3 = strlen (s3);
char *s = my_malloc ((n1 + n2 + n3 + 1) * sizeof (char));
s[n1 + n2 + n3] = '\0';
memcpy (s, s1, n1);
memcpy (s + n1, s2, n2);
memcpy (s + n1 + n2, s3, n3);
return s;
}
char *
string_repeat (size_t n, const char *s)
{
/* This brute force implementation will suffice. */
char *t = "";
for (size_t i = 0; i != n; i += 1)
t = string_append2 (t, s);
return t;
}
/*------------------------------------------------------------------*/
typedef mpz_t *cf_func (size_t i, void *env);
struct cf
{
bool terminated; /* Are there no more terms? */
size_t m; /* The number of terms memoized. */
size_t n; /* The size of memoization storage. */
mpz_t **memo; /* Memoization storage. */
cf_func *func; /* A function that produces terms. */
void *env; /* An environment for func. */
};
typedef struct cf *cf_t;
cf_t
make_cf (cf_func *func, void *env)
{
cf_t cf = my_malloc (sizeof (struct cf));
cf->terminated = false;
cf->m = 0;
cf->n = 32;
cf->memo = my_malloc (cf->n * sizeof (mpz_t *));
cf->func = func;
cf->env = env;
return cf;
}
void
resize_cf (cf_t cf, size_t minimum)
{
/* Ensure there is at least twice the minimum storage. */
size_t size = 2 * minimum;
if (cf->n < size)
{
cf->memo = my_realloc (cf->memo, size * sizeof (mpz_t *));
cf->n = size;
}
}
void
update_cf (cf_t cf, size_t needed)
{
/* Ensure there are at least a certain number of finite terms
memoized (or else that all of them are memoized). */
if (!cf->terminated && cf->m < needed)
{
if (cf->n < needed)
resize_cf (cf, needed);
while (!cf->terminated && cf->m != needed)
{
cf->memo[cf->m] = cf->func (cf->m, cf->env);
cf->m += 1;
}
}
}
mpz_t *
cf_ref (cf_t cf, size_t i)
{
/* Get the ith term, or a NULL pointer if there is no finite ith
term. */
update_cf (cf, i + 1);
return (i < cf->m) ? cf->memo[i] : NULL;
}
size_t default_max_terms = 20;
char *
cf2string (cf_t cf, size_t max_terms)
{
if (max_terms == 0)
max_terms = default_max_terms;
size_t i = 0;
char *s = string_append1 ("[");
bool done = false;
while (!done)
{
mpz_t *term = cf_ref (cf, i);
if (term == NULL)
{
s = string_append2 (s, "]");
done = true;
}
else if (i == max_terms)
{
s = string_append2 (s, ",...]");
done = true;
}
else
{
static const char *separators[3] = { "", ";", "," };
const char *separator = separators[(i <= 1) ? i : 2];
const char *term_str = mpz_get_str (NULL, 10, *term);
s = string_append3 (s, separator, term_str);
i += 1;
}
}
return s;
}
/*------------------------------------------------------------------*/
cf_t golden_ratio;
cf_t silver_ratio;
mpz_t *
return_constant ([[maybe_unused]] size_t i, void *env)
{
mpz_t *term = my_malloc (sizeof (mpz_t));
mpz_init_set (*term, *((mpz_t *) env));
return term;
}
cf_t
make_cf_with_constant_terms (int term_si)
{
mpz_t *env = my_malloc (sizeof (mpz_t));
mpz_init_set_si (*env, term_si);
return make_cf (return_constant, env);
}
/*------------------------------------------------------------------*/
cf_t sqrt2;
mpz_t *
return_sqrt2_term (size_t i, [[maybe_unused]] void *env)
{
mpz_t *term = my_malloc (sizeof (mpz_t));
mpz_init_set_si (*term, (i == 0) ? 1 : 2);
return term;
}
cf_t
make_cf_sqrt2 (void)
{
return make_cf (return_sqrt2_term, NULL);
}
/*------------------------------------------------------------------*/
cf_t frac_13_11;
cf_t frac_22_7;
cf_t one;
cf_t two;
cf_t three;
cf_t four;
mpz_t *
return_rational_term ([[maybe_unused]] size_t i, void *env)
{
mpz_t *frac = env;
mpz_t *term = NULL;
if (mpz_sgn (frac[1]) != 0)
{
term = my_malloc (sizeof (mpz_t));
mpz_init (*term);
mpz_t r;
mpz_init (r);
mpz_fdiv_qr (*term, r, frac[0], frac[1]);
mpz_set (frac[0], frac[1]);
mpz_set (frac[1], r);
}
return term;
}
cf_t
make_cf_rational (int numerator_si, int denominator_si)
{
mpz_t *env = my_malloc (2 * sizeof (mpz_t));
mpz_init_set_si (env[0], numerator_si);
mpz_init_set_si (env[1], denominator_si);
return make_cf (return_rational_term, env);
}
cf_t
make_cf_integer (int integer_si)
{
return make_cf_rational (integer_si, 1);
}
/*------------------------------------------------------------------*/
/* Thresholds. */
mpz_t number_that_is_too_big;
mpz_t practically_infinite;
struct ng8_env
{
mpz_t ng[8];
cf_t x;
cf_t y;
size_t ix;
size_t iy;
bool overflow;
};
typedef struct ng8_env *ng8_env_t;
enum ng8_index
{
ng8a12 = 0,
ng8a1 = 1,
ng8a2 = 2,
ng8a = 3,
ng8b12 = 4,
ng8b1 = 5,
ng8b2 = 6,
ng8b = 7
};
static bool
ng8_too_big (const mpz_t ng[8])
{
bool too_big = false;
int i = 0;
while (!too_big && i != 8)
{
too_big = (0 <= mpz_cmpabs (ng[i], number_that_is_too_big));
i += 1;
}
return too_big;
}
static bool
treat_ng8_term_as_infinite (const mpz_t term)
{
return (0 <= mpz_cmpabs (term, practically_infinite));
}
static void
a_plus_bc (mpz_t result, const mpz_t a, const mpz_t b,
const mpz_t c)
{
mpz_set (result, a);
mpz_addmul (result, b, c);
}
static void
abc (mpz_t result, const mpz_t a, const mpz_t b, const mpz_t c)
{
mpz_mul (result, a, b);
mpz_mul (result, result, c);
}
static void
absorb_x_term (ng8_env_t env)
{
mpz_t tmp[8];
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->overflow) ? cf_ref (env->x, env->ix) : NULL;
env->ix += 1;
mpz_set (env->ng[ng8a2], tmp[ng8a12]);
mpz_set (env->ng[ng8a], tmp[ng8a1]);
mpz_set (env->ng[ng8b2], tmp[ng8b12]);
mpz_set (env->ng[ng8b], tmp[ng8b1]);
if (term != NULL)
{
a_plus_bc (env->ng[ng8a12], tmp[ng8a2], tmp[ng8a12], *term);
a_plus_bc (env->ng[ng8a1], tmp[ng8a], tmp[ng8a1], *term);
a_plus_bc (env->ng[ng8b12], tmp[ng8b2], tmp[ng8b12], *term);
a_plus_bc (env->ng[ng8b1], tmp[ng8b], tmp[ng8b1], *term);
if (ng8_too_big (env->ng))
{
env->overflow = true;
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a1], tmp[ng8a1]);
mpz_set (env->ng[ng8b12], tmp[ng8b12]);
mpz_set (env->ng[ng8b1], tmp[ng8b1]);
}
}
}
static void
absorb_y_term (ng8_env_t env)
{
mpz_t tmp[8];
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->overflow) ? cf_ref (env->y, env->iy) : NULL;
env->iy += 1;
mpz_set (env->ng[ng8a1], tmp[ng8a12]);
mpz_set (env->ng[ng8a], tmp[ng8a2]);
mpz_set (env->ng[ng8b1], tmp[ng8b12]);
mpz_set (env->ng[ng8b], tmp[ng8b2]);
if (term != NULL)
{
a_plus_bc (env->ng[ng8a12], tmp[ng8a1], tmp[ng8a12], *term);
a_plus_bc (env->ng[ng8a2], tmp[ng8a], tmp[ng8a2], *term);
a_plus_bc (env->ng[ng8b12], tmp[ng8b1], tmp[ng8b12], *term);
a_plus_bc (env->ng[ng8b2], tmp[ng8b], tmp[ng8b2], *term);
if (ng8_too_big (env->ng))
{
env->overflow = true;
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a2], tmp[ng8a2]);
mpz_set (env->ng[ng8b12], tmp[ng8b12]);
mpz_set (env->ng[ng8b2], tmp[ng8b2]);
}
}
}
mpz_t *
return_ng8_term ([[maybe_unused]] size_t i, void *env)
{
/* The McCabe complexity of this function is high. Please be careful
if modifying the code. */
ng8_env_t p = env;
mpz_t *term = NULL;
bool done = false;
while (!done)
{
const bool b12_zero = (mpz_sgn (p->ng[ng8b12]) == 0);
const bool b1_zero = (mpz_sgn (p->ng[ng8b1]) == 0);
const bool b2_zero = (mpz_sgn (p->ng[ng8b2]) == 0);
const bool b_zero = (mpz_sgn (p->ng[ng8b]) == 0);
if (b_zero && b1_zero && b2_zero && b12_zero)
done = true; /* There are no more terms. */
else if (b_zero && b2_zero)
absorb_x_term (p);
else if (b_zero || b2_zero)
absorb_y_term (p);
else if (b1_zero)
absorb_x_term (p);
else
{
mpz_t q, r;
mpz_inits (q, r, NULL);
mpz_t q1, r1;
mpz_inits (q1, r1, NULL);
mpz_t q2, r2;
mpz_inits (q2, r2, NULL);
mpz_t q12, r12;
mpz_inits (q12, r12, NULL);
mpz_fdiv_qr (q, r, p->ng[ng8a], p->ng[ng8b]);
mpz_fdiv_qr (q1, r1, p->ng[ng8a1], p->ng[ng8b1]);
mpz_fdiv_qr (q2, r2, p->ng[ng8a2], p->ng[ng8b2]);
if (!b12_zero)
mpz_fdiv_qr (q12, r12, p->ng[ng8a12], p->ng[ng8b12]);
if (!b12_zero
&& mpz_cmp (q, q1) == 0
&& mpz_cmp (q, q2) == 0
&& mpz_cmp (q, q12) == 0)
{
// Output a term.
mpz_set (p->ng[ng8a12], p->ng[ng8b12]);
mpz_set (p->ng[ng8a1], p->ng[ng8b1]);
mpz_set (p->ng[ng8a2], p->ng[ng8b2]);
mpz_set (p->ng[ng8a], p->ng[ng8b]);
mpz_set (p->ng[ng8b12], r12);
mpz_set (p->ng[ng8b1], r1);
mpz_set (p->ng[ng8b2], r2);
mpz_set (p->ng[ng8b], r);
if (!treat_ng8_term_as_infinite (q))
{
term = my_malloc (sizeof (mpz_t));
mpz_init_set (*term, q);
}
done = true;
}
else
{
/* Rather than compare fractions, we will put the
numerators over a common denominator of b*b1*b2, and
then compare the new numerators. */
mpz_t n, n1, n2, n1_diff, n2_diff;
mpz_inits (n, n1, n2, n1_diff, n2_diff, NULL);
abc (n, p->ng[ng8a], p->ng[ng8b1], p->ng[ng8b2]);
abc (n1, p->ng[ng8a1], p->ng[ng8b], p->ng[ng8b2]);
abc (n2, p->ng[ng8a2], p->ng[ng8b], p->ng[ng8b1]);
mpz_sub (n1_diff, n1, n);
mpz_sub (n2_diff, n2, n);
if (mpz_cmpabs (n1_diff, n2_diff) > 0)
absorb_x_term (p);
else
absorb_y_term (p);
}
}
}
return term;
}
cf_t
make_cf_ng8 (int ng[8], cf_t x, cf_t y)
{
ng8_env_t env = my_malloc (sizeof (struct ng8_env));
for (int i = 0; i != 8; i += 1)
mpz_init_set_si (env->ng[i], ng[i]);
env->x = x;
env->y = y;
env->ix = 0;
env->iy = 0;
env->overflow = false;
return make_cf (return_ng8_term, env);
}
/*------------------------------------------------------------------*/
static void *
gmp_malloc (size_t alloc_size)
{
return my_malloc (alloc_size);
}
static void *
gmp_realloc (void *p,
[[maybe_unused]] size_t old_size,
size_t new_size)
{
return my_realloc (p, new_size);
}
static void
gmp_free (void *p, [[maybe_unused]] size_t size)
{
/* There is no need for us to explicitly free memory, and
performance might even suffer if we do. On the other hand, maybe
GMP will free memory that otherwise would have been passed over
for collection. */
my_free (p); /* <-- optional */
}
void
show (const char *expression, cf_t cf, const char *note)
{
size_t nexpr = strlen (expression);
char *padding = string_repeat (MAX (19, nexpr + 1) - nexpr, " ");
char *line = string_append3 (padding, expression, " => ");
char *cfstr = cf2string (cf, 0);
line = string_append2 (line, cfstr);
if (note != NULL)
{
size_t ncfstr = strlen (cfstr);
padding = string_repeat (MAX (48, ncfstr + 1) - ncfstr, " ");
line = string_append3 (line, padding, note);
}
puts (line);
}
int ng8_add[8] = { 0, 1, 1, 0, 0, 0, 0, 1 };
int ng8_sub[8] = { 0, 1, -1, 0, 0, 0, 0, 1 };
int ng8_mul[8] = { 1, 0, 0, 0, 0, 0, 0, 1 };
int ng8_div[8] = { 0, 1, 0, 0, 0, 0, 1, 0 };
int
main (void)
{
GC_INIT ();
/* GMP has to be told to use Boehm GC as its allocator. */
mp_set_memory_functions (gmp_malloc, gmp_realloc, gmp_free);
/* Initialize thresholds, to values chosen merely for
demonstration. */
mpz_init_set_si (number_that_is_too_big, 1);
mpz_mul_2exp (number_that_is_too_big, number_that_is_too_big,
512); /* 2**512 */
mpz_init_set_si (practically_infinite, 1);
mpz_mul_2exp (practically_infinite, practically_infinite,
64); /* 2**64 */
/* Initialize global continued fractions. */
golden_ratio = make_cf_with_constant_terms (1);
silver_ratio = make_cf_with_constant_terms (2);
sqrt2 = make_cf_sqrt2 ();
frac_13_11 = make_cf_rational (13, 11);
frac_22_7 = make_cf_rational (22, 7);
one = make_cf_integer (1);
two = make_cf_integer (2);
three = make_cf_integer (3);
four = make_cf_integer (4);
/* Divide the silver ratio by 2 times the square root of 2. */
int ng8_method1[8] = { 0, 1, 0, 0, 0, 0, 2, 0 };
cf_t method1 = make_cf_ng8 (ng8_method1, silver_ratio, sqrt2);
/* Add 1/8 to 1/8th of the square of the silver ratio. */
int ng8_method2[8] = { 1, 0, 0, 1, 0, 0, 0, 8 };
cf_t method2 = make_cf_ng8 (ng8_method2, silver_ratio,
silver_ratio);
/* Thrice divide the silver ratio by the square root of 2. */
cf_t method3 = make_cf_ng8 (ng8_div, silver_ratio, sqrt2);
method3 = make_cf_ng8 (ng8_div, method3, sqrt2);
method3 = make_cf_ng8 (ng8_div, method3, sqrt2);
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "1 + sqrt(2)");
show ("sqrt(2)", sqrt2, NULL);
show ("13/11", frac_13_11, NULL);
show ("22/7", frac_22_7, NULL);
show ("one", one, NULL);
show ("two", two, NULL);
show ("three", three, NULL);
show ("four", four, NULL);
show ("(1 + 1/sqrt(2))/2", method1, "method 1");
show ("(1 + 1/sqrt(2))/2", method2, "method 2");
show ("(1 + 1/sqrt(2))/2", method3, "method 3");
show ("sqrt(2) + sqrt(2)", make_cf_ng8 (ng8_add, sqrt2, sqrt2),
NULL);
show ("sqrt(2) - sqrt(2)", make_cf_ng8 (ng8_sub, sqrt2, sqrt2),
NULL);
show ("sqrt(2) * sqrt(2)", make_cf_ng8 (ng8_mul, sqrt2, sqrt2),
NULL);
show ("sqrt(2) / sqrt(2)", make_cf_ng8 (ng8_div, sqrt2, sqrt2),
NULL);
return 0;
}
/*------------------------------------------------------------------*/
- Output:
$ gcc -std=gnu2x -Wall -Wextra -g bivariate-continued-fraction-task-gmp.c -lgmp -lgc && ./a.out
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 1 + sqrt(2)
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
one => [1]
two => [2]
three => [3]
four => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
C++
Uses matrixNG, NG_4 and NG from Continued_fraction/Arithmetic/G(matrix_NG,_Contined_Fraction_N)#C++, and r2cf from Continued_fraction/Arithmetic/Construct_from_rational_number#C++
/* Implement matrix NG
Nigel Galloway, February 12., 2013
*/
class NG_8 : public matrixNG {
private: int a12, a1, a2, a, b12, b1, b2, b, t;
double ab, a1b1, a2b2, a12b12;
const int chooseCFN(){return fabs(a1b1-ab) > fabs(a2b2-ab)? 0 : 1;}
const bool needTerm() {
if (b1==0 and b==0 and b2==0 and b12==0) return false;
if (b==0){cfn = b2==0? 0:1; return true;} else ab = ((double)a)/b;
if (b2==0){cfn = 1; return true;} else a2b2 = ((double)a2)/b2;
if (b1==0){cfn = 0; return true;} else a1b1 = ((double)a1)/b1;
if (b12==0){cfn = chooseCFN(); return true;} else a12b12 = ((double)a12)/b12;
thisTerm = (int)ab;
if (thisTerm==(int)a1b1 and thisTerm==(int)a2b2 and thisTerm==(int)a12b12){
t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm; t=a2; a2=b2; b2=t-b2*thisTerm; t=a12; a12=b12; b12=t-b12*thisTerm;
haveTerm = true; return false;
}
cfn = chooseCFN();
return true;
}
void consumeTerm(){if(cfn==0){a=a1; a2=a12; b=b1; b2=b12;} else{a=a2; a1=a12; b=b2; b1=b12;}}
void consumeTerm(int n){
if(cfn==0){t=a; a=a1; a1=t+a1*n; t=a2; a2=a12; a12=t+a12*n; t=b; b=b1; b1=t+b1*n; t=b2; b2=b12; b12=t+b12*n;}
else{t=a; a=a2; a2=t+a2*n; t=a1; a1=a12; a12=t+a12*n; t=b; b=b2; b2=t+b2*n; t=b1; b1=b12; b12=t+b12*n;}
}
public:
NG_8(int a12, int a1, int a2, int a, int b12, int b1, int b2, int b): a12(a12), a1(a1), a2(a2), a(a), b12(b12), b1(b1), b2(b2), b(b){
}};
Testing
[3;7] + [0;2]
int main() {
NG_8 a(0,1,1,0,0,0,0,1);
r2cf n2(22,7);
r2cf n1(1,2);
for(NG n(&a, &n1, &n2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
NG_4 a3(2,1,0,2);
r2cf n3(22,7);
for(NG n(&a3, &n3); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 1 1 1 4 3 1 1 1 4
[1:5,2] * [3;7]
int main() {
NG_8 b(1,0,0,0,0,0,0,1);
r2cf b1(13,11);
r2cf b2(22,7);
for(NG n(&b, &b1, &b2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(NG n(&a, &b2, &b1); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(286,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 1 2 2 3 1 2 2
[1:5,2] - [3;7]
int main() {
NG_8 c(0,1,-1,0,0,0,0,1);
r2cf c1(13,11);
r2cf c2(22,7);
for(NG n(&c, &c1, &c2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(-151,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
-1 -1 -24 -1 -2 -1 -1 -24 -1 -2
Divide [] by [3;7]
int main() {
NG_8 d(0,1,0,0,0,0,1,0);
r2cf d1(22*22,7*7);
r2cf d2(22,7);
for(NG n(&d, &d1, &d2); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
- Output:
3 7
([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])
int main() {
r2cf a1(2,7);
r2cf a2(13,11);
NG_8 na(0,1,1,0,0,0,0,1);
NG A(&na, &a1, &a2); //[0;3,2] + [1;5,2]
r2cf b1(2,7);
r2cf b2(13,11);
NG_8 nb(0,1,-1,0,0,0,0,1);
NG B(&nb, &b1, &b2); //[0;3,2] - [1;5,2]
NG_8 nc(1,0,0,0,0,0,0,1); //A*B
for(NG n(&nc, &A, &B); n.moreTerms(); std::cout << n.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(2,7); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(13,11); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
for(r2cf cf(-7797,5929); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
std::cout << std::endl;
return 0;
}
Common Lisp
(defstruct (cf (:conc-name %%cf-)
(:constructor make-cf (generator)))
"continued fraction"
(generator nil :type function)
(terminated-p nil :type boolean)
(memo (make-array '(32)) :type (array integer))
(memo-count 0 :type fixnum))
(defstruct (ng8 (:constructor ng8 (a12 a1 a2 a
b12 b1 b2 b)))
"coefficients of a bihomographic function"
(a12 0 :type integer)
(a1 0 :type integer)
(a2 0 :type integer)
(a 0 :type integer)
(b12 0 :type integer)
(b1 0 :type integer)
(b2 0 :type integer)
(b 0 :type integer))
(defun cf-ref (cf i)
"Return the ith term, or nil if there is none."
(declare (cf cf) (fixnum i))
(defun get-more-terms (needed)
(declare (fixnum needed))
(loop while (not (%%cf-terminated-p cf))
while (< (%%cf-memo-count cf) needed)
do (let ((term (funcall (%%cf-generator cf))))
(cond (term (let ((memo (%%cf-memo cf))
(m (%%cf-memo-count cf)))
(setf (aref memo m) term)
(setf (%%cf-memo-count cf) (1+ m))))
(t (setf (%%cf-terminated-p cf) t))))))
(defun update (needed)
(declare (fixnum needed))
(cond ((%%cf-terminated-p cf) (progn))
((<= needed (%%cf-memo-count cf)) (progn))
((<= needed (array-dimension (%%cf-memo cf) 0))
(get-more-terms needed))
(t (let* ((n1 (+ needed needed))
(memo1 (make-array (list n1))))
(loop for i from 0 upto (1- (%%cf-memo-count cf))
do (setf (aref memo1 i) (aref (%%cf-memo cf) i)))
(setf (%%cf-memo cf) memo1)
(get-more-terms needed)))))
(update (1+ i))
(the (or integer null) (and (< i (%%cf-memo-count cf))
(aref (%%cf-memo cf) i))))
(defparameter *cf-max-terms* 20
"Default term-count limit for cf->string.")
(defun cf->string (cf &optional (max-terms *cf-max-terms*))
"Make a readable string from a continued fraction."
(declare (cf cf))
(loop with i = 0
with s = "["
do (let ((term (cf-ref cf i)))
(cond ((not term)
(return (concatenate 'string s "]")))
((= i max-terms)
(return (concatenate 'string s ",...]")))
(t (let ((separator (case i
((0) "")
((1) ";")
(t ",")))
(term-str (format nil "~A" term)))
(setf i (1+ i))
(setf s (concatenate 'string s separator
term-str))))))))
(defun integer->cf (num)
"Transform an integer to a continued fraction."
(declare (integer num))
(let ((terminated-p nil))
(declare (boolean terminated-p))
(make-cf #'(lambda ()
(and (not terminated-p)
(progn (setf terminated-p t)
num))))))
(defun ratio->cf (num)
"Transform a ratio to a continued fraction."
(declare (ratio num))
(let ((n (the integer (numerator num)))
(d (the integer (denominator num))))
(make-cf #'(lambda ()
(and (not (zerop d))
(multiple-value-bind (q r) (floor n d)
(setf n d)
(setf d r)
q))))))
;; Thresholds chosen merely for demonstration.
(defparameter number-that-is-too-big (expt 2 512))
(defparameter practically-infinite (expt 2 64))
(defun num-too-big-p (num)
(declare (integer num))
(>= (abs num) (abs number-that-is-too-big)))
(defun ng8-too-big-p (ng)
(declare (ng8 ng))
(or (num-too-big-p (ng8-a12 ng))
(num-too-big-p (ng8-a1 ng))
(num-too-big-p (ng8-a2 ng))
(num-too-big-p (ng8-a ng))
(num-too-big-p (ng8-b12 ng))
(num-too-big-p (ng8-b1 ng))
(num-too-big-p (ng8-b2 ng))
(num-too-big-p (ng8-b ng))))
(defun treat-as-infinite-p (term)
(declare (integer term))
(>= (abs term) (abs practically-infinite)))
(defun quotient (u v)
(declare (integer u v))
(if (zerop v)
(list nil nil)
(multiple-value-list (floor u v))))
(defmacro absorb-x-term (ng xsource ix)
`(let ((a12 (ng8-a12 ,ng))
(a1 (ng8-a1 ,ng))
(a2 (ng8-a2 ,ng))
(a (ng8-a ,ng))
(b12 (ng8-b12 ,ng))
(b1 (ng8-b1 ,ng))
(b2 (ng8-b2 ,ng))
(b (ng8-b ,ng))
(term (funcall ,xsource ,ix)))
(setf ,ix (1+ ,ix))
(if term
(let ((ng^ (ng8 (+ a2 (* a12 term))
(+ a (* a1 term)) a12 a1
(+ b2 (* b12 term))
(+ b (* b1 term)) b12 b1)))
(if (not (ng8-too-big-p ng^))
(setf ,ng ng^)
(progn (setf ,ng (ng8 a12 a1 a12 a1 b12 b1 b12 b1))
;; Replace the x source with one that never
;; returns a term.
(setf ,xsource #'(lambda (i) nil)))))
(setf ,ng (ng8 a12 a1 a12 a1 b12 b1 b12 b1)))))
(defmacro absorb-y-term (ng ysource iy)
`(let ((a12 (ng8-a12 ,ng))
(a1 (ng8-a1 ,ng))
(a2 (ng8-a2 ,ng))
(a (ng8-a ,ng))
(b12 (ng8-b12 ,ng))
(b1 (ng8-b1 ,ng))
(b2 (ng8-b2 ,ng))
(b (ng8-b ,ng))
(term (funcall ,ysource ,iy)))
(setf ,iy (1+ ,iy))
(if term
(let ((ng^ (ng8 (+ a1 (* a12 term)) a12
(+ a (* a2 term)) a2
(+ b1 (* b12 term)) b12
(+ b (* b2 term)) b2)))
(if (not (ng8-too-big-p ng^))
(setf ,ng ng^)
(progn (setf ,ng (ng8 a12 a12 a2 a2 b12 b12 b2 b2))
;; Replace the y source with one that never
;; returns a term.
(setf ysource #'(lambda (i) nil)))))
(setf ,ng (ng8 a12 a12 a2 a2 b12 b12 b2 b2)))))
(defun apply-ng8 (ng8 x y)
(declare (ng8 ng8))
(let ((ng ng8)
(xsource #'(lambda (i) (cf-ref x i)))
(ysource #'(lambda (i) (cf-ref y i)))
(ix 0)
(iy 0))
(flet
((main ()
(loop
with absorb
for bzero = (zerop (ng8-b ng))
for b1zero = (zerop (ng8-b1 ng))
for b2zero = (zerop (ng8-b2 ng))
for b12zero = (zerop (ng8-b12 ng))
do (multiple-value-bind (q r q1 r1 q2 r2 q12 r12)
(values-list
`(,@(quotient (ng8-a ng) (ng8-b ng))
,@(quotient (ng8-a1 ng) (ng8-b1 ng))
,@(quotient (ng8-a2 ng) (ng8-b2 ng))
,@(quotient (ng8-a12 ng) (ng8-b12 ng))))
(cond
((and bzero b1zero b2zero b12zero) (return nil))
((and bzero b2zero) (setf absorb 'x))
((or bzero b2zero) (setf absorb 'y))
(b1zero (setf absorb 'x))
((and (not b12zero) (= q q1 q2 q12))
;;
;; Output a term.
;;
(setf ng (ng8 (ng8-b12 ng) (ng8-b1 ng)
(ng8-b2 ng) (ng8-b ng)
r12 r1 r2 r))
(return (and (not (treat-as-infinite-p q)) q)))
(t
;;
;; Rather than compare fractions, we will put the
;; numerators over a common denominator of
;; b*b1*b2, and then compare the new numerators.
;;
(let ((n (* (ng8-a ng) (ng8-b1 ng) (ng8-b2 ng)))
(n1 (* (ng8-a1 ng) (ng8-b ng) (ng8-b2 ng)))
(n2 (* (ng8-a2 ng) (ng8-b ng) (ng8-b1 ng))))
(if (> (abs (- n1 n)) (abs (- n2 n)))
(setf absorb 'x)
(setf absorb 'y))))))
when (eq absorb 'x)
do (absorb-x-term ng xsource ix)
when (eq absorb 'y)
do (absorb-y-term ng ysource iy))))
(make-cf #'main))))
(defun show (expression cf &optional (note ""))
(format t "~A => ~A~A~%" expression (cf->string cf) note))
(defvar golden-ratio (make-cf #'(lambda () 1)))
(defvar silver-ratio (make-cf #'(lambda () 2)))
(defvar sqrt2 (make-cf (let ((next-term 1))
#'(lambda ()
(let ((term next-term))
(setf next-term 2)
term)))))
(defvar frac13/11 (ratio->cf 13/11))
(defvar frac22/7 (ratio->cf 22/7))
(defvar one (integer->cf 1))
(defvar two (integer->cf 2))
(defvar three (integer->cf 3))
(defvar four (integer->cf 4))
(defun cf+ (x y) (apply-ng8 (ng8 0 1 1 0 0 0 0 1) x y))
(defun cf- (x y) (apply-ng8 (ng8 0 1 -1 0 0 0 0 1) x y))
(defun cf* (x y) (apply-ng8 (ng8 1 0 0 0 0 0 0 1) x y))
(defun cf/ (x y) (apply-ng8 (ng8 0 1 0 0 0 0 1 0) x y))
(show " golden ratio" golden-ratio)
(show " silver ratio" silver-ratio)
(show " sqrt(2)" sqrt2)
(show " 13/11" frac13/11)
(show " 22/7" frac22/7)
(show " 1" one)
(show " 2" two)
(show " 3" three)
(show " 4" four)
(show " (1 + 1/sqrt(2))/2" (cf/ silver-ratio
(cf* sqrt2 (cf* sqrt2 sqrt2)))
" method 1")
(show " (1 + 1/sqrt(2))/2" (apply-ng8 (ng8 1 0 0 1 0 0 0 8)
silver-ratio
silver-ratio)
" method 2")
(show " (1 + 1/sqrt(2))/2" (cf/ (cf/ (cf/ silver-ratio sqrt2)
sqrt2)
sqrt2)
" method 3")
(show " sqrt(2) + sqrt(2)" (cf+ sqrt2 sqrt2))
(show " sqrt(2) - sqrt(2)" (cf- sqrt2 sqrt2))
(show " sqrt(2) * sqrt(2)" (cf* sqrt2 sqrt2))
(show " sqrt(2) / sqrt(2)" (cf/ sqrt2 sqrt2))
- Output:
$ sbcl --script bivariate-continued-fraction-task.lisp
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
D
//--------------------------------------------------------------------
import std.algorithm;
import std.bigint;
import std.conv;
import std.math;
import std.stdio;
import std.string;
import std.typecons;
//--------------------------------------------------------------------
class CF // Continued fraction.
{
alias Term = BigInt; // The type for terms.
alias Index = size_t; // The type for indexing terms.
protected bool terminated; // Are there no more terms?
protected size_t m; // The number of terms memoized.
private Term[] memo; // Memoization storage.
static Index maxTerms = 20; // Maximum number of terms in the
// string representation.
this ()
{
terminated = false;
m = 0;
memo.length = 32;
}
protected Nullable!Term generate ()
{
// Return terms for zero. To get different terms, override this
// method.
auto retval = (Term(0)).nullable;
if (m != 0)
retval.nullify();
return retval;
}
public Nullable!Term opIndex (Index i)
{
void update (size_t needed)
{
// Ensure all finite terms with indices 0 <= i < needed are
// memoized.
if (!terminated && m < needed)
{
if (memo.length < needed)
// To reduce the frequency of reallocation, increase the
// space to twice what might be needed right now.
memo.length = 2 * needed;
while (m != needed && !terminated)
{
auto term = generate ();
if (!term.isNull())
{
memo[m] = term.get();
m += 1;
}
else
terminated = true;
}
}
}
update (i + 1);
Nullable!Term retval;
if (i < m)
retval = memo[i].nullable;
return retval;
}
public override string toString ()
{
static string[3] separators = ["", ";", ","];
string s = "[";
Index i = 0;
bool done = false;
while (!done)
{
auto term = this[i];
if (term.isNull())
{
s ~= "]";
done = true;
}
else if (i == maxTerms)
{
s ~= ",...]";
done = true;
}
else
{
s ~= separators[(i <= 1) ? i : 2];
s ~= to!string (term.get());
i += 1;
}
}
return s;
}
public CF opBinary(string op : "+") (CF other)
{
return new cfNG8 (ng8_add, this, other);
}
public CF opBinary(string op : "-") (CF other)
{
return new cfNG8 (ng8_sub, this, other);
}
public CF opBinary(string op : "*") (CF other)
{
return new cfNG8 (ng8_mul, this, other);
}
public CF opBinary(string op : "/") (CF other)
{
return new cfNG8 (ng8_div, this, other);
}
};
//--------------------------------------------------------------------
class cfIndexed : CF // Continued fraction with an index-to-term map.
{
alias Mapper = Nullable!Term delegate (Index);
protected Mapper map;
this (Mapper map)
{
this.map = map;
}
protected override Nullable!Term generate ()
{
return map (m);
}
}
__gshared goldenRatio =
new cfIndexed ((i) => CF.Term(1).nullable);
__gshared silverRatio =
new cfIndexed ((i) => CF.Term(2).nullable);
__gshared sqrt2 =
new cfIndexed ((i) => CF.Term(min (i + 1, 2)).nullable);
//--------------------------------------------------------------------
class cfRational : CF // CF for a rational number.
{
private Term n;
private Term d;
this (Term numer, Term denom = Term(1))
{
n = numer;
d = denom;
}
protected override Nullable!Term generate ()
{
Nullable!Term term;
if (d != 0)
{
auto q = n / d;
auto r = n % d;
n = d;
d = r;
term = q.nullable;
}
return term;
}
}
__gshared frac_13_11 = new cfRational (CF.Term(13), CF.Term(11));
__gshared frac_22_7 = new cfRational (CF.Term(22), CF.Term(7));
__gshared one = new cfRational (CF.Term(1));
__gshared two = new cfRational (CF.Term(2));
__gshared three = new cfRational (CF.Term(3));
__gshared four = new cfRational (CF.Term(4));
//--------------------------------------------------------------------
class NG8 // Bihomographic function.
{
public CF.Term a12, a1, a2, a;
public CF.Term b12, b1, b2, b;
this (CF.Term a12, CF.Term a1, CF.Term a2, CF.Term a,
CF.Term b12, CF.Term b1, CF.Term b2, CF.Term b)
{
this.a12 = a12;
this.a1 = a1;
this.a2 = a2;
this.a = a;
this.b12 = b12;
this.b1 = b1;
this.b2 = b2;
this.b = b;
}
this (long a12, long a1, long a2, long a,
long b12, long b1, long b2, long b)
{
this.a12 = a12;
this.a1 = a1;
this.a2 = a2;
this.a = a;
this.b12 = b12;
this.b1 = b1;
this.b2 = b2;
this.b = b;
}
this (NG8 other)
{
this.a12 = other.a12;
this.a1 = other.a1;
this.a2 = other.a2;
this.a = other.a;
this.b12 = other.b12;
this.b1 = other.b1;
this.b2 = other.b2;
this.b = other.b;
}
}
class cfNG8 : CF // CF that is a bihomographic function of other CF.
{
private NG8 ng;
private CF x;
private CF y;
private Index ix;
private Index iy;
private bool overflow;
//
// Thresholds chosen merely for demonstration.
//
static number_that_is_too_big = // 2 ** 512
BigInt ("13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096");
static practically_infinite = // 2 ** 64
BigInt ("18446744073709551616");
this (NG8 ng, CF x, CF y)
{
this.ng = new NG8 (ng);
this.x = x;
this.y = y;
ix = 0;
iy = 0;
overflow = false;
}
protected override Nullable!Term generate ()
{
// The McCabe complexity of this function is high. Please be
// careful if modifying the code.
Nullable!Term term;
bool done = false;
while (!done)
{
bool bz = (ng.b == 0);
bool b1z = (ng.b1 == 0);
bool b2z = (ng.b2 == 0);
bool b12z = (ng.b12 == 0);
if (bz && b1z && b2z && b12z)
done = true; // There are no more terms.
else if (bz && b2z)
absorb_x_term ();
else if (bz || b2z)
absorb_y_term ();
else if (b1z)
absorb_x_term ();
else
{
Term q, r;
Term q1, r1;
Term q2, r2;
Term q12, r12;
divMod (ng.a, ng.b, q, r);
divMod (ng.a1, ng.b1, q1, r1);
divMod (ng.a2, ng.b2, q2, r2);
if (ng.b12 != 0)
divMod (ng.a12, ng.b12, q12, r12);
if (!b12z && q == q1 && q == q2 && q == q12)
{
// Output a term.
ng = new NG8 (ng.b12, ng.b1, ng.b2, ng.b,
r12, r1, r2, r);
if (!treat_as_infinite (q))
term = q.nullable;
done = true;
}
else
{
//
// Rather than compare fractions, we will put the
// numerators over a common denominator of b*b1*b2,
// and then compare the new numerators.
//
Term n = ng.a * ng.b1 * ng.b2;
Term n1 = ng.a1 * ng.b * ng.b2;
Term n2 = ng.a2 * ng.b * ng.b1;
if (abs (n1 - n) > abs (n2 - n))
absorb_x_term ();
else
absorb_y_term ();
}
}
}
return term;
}
private void absorb_x_term ()
{
Nullable!Term term;
if (!overflow)
term = x[ix];
ix += 1;
if (!term.isNull())
{
auto t = term.get();
auto new_ng = new NG8 (ng.a2 + (ng.a12 * t),
ng.a + (ng.a1 * t),
ng.a12, ng.a1,
ng.b2 + (ng.b12 * t),
ng.b + (ng.b1 * t),
ng.b12, ng.b1);
if (!too_big (new_ng))
ng = new_ng;
else
{
ng = new NG8 (ng.a12, ng.a1, ng.a12, ng.a1,
ng.b12, ng.b1, ng.b12, ng.b1);
overflow = true;
}
}
else
ng = new NG8 (ng.a12, ng.a1, ng.a12, ng.a1,
ng.b12, ng.b1, ng.b12, ng.b1);
}
private void absorb_y_term ()
{
Nullable!Term term;
if (!overflow)
term = y[iy];
iy += 1;
if (!term.isNull())
{
auto t = term.get();
auto new_ng = new NG8 (ng.a1 + (ng.a12 * t), ng.a12,
ng.a + (ng.a2 * t), ng.a2,
ng.b1 + (ng.b12 * t), ng.b12,
ng.b + (ng.b2 * t), ng.b2);
if (!too_big (new_ng))
ng = new_ng;
else
{
ng = new NG8 (ng.a12, ng.a12, ng.a2, ng.a2,
ng.b12, ng.b12, ng.b2, ng.b2);
overflow = true;
}
}
else
ng = new NG8 (ng.a12, ng.a12, ng.a2, ng.a2,
ng.b12, ng.b12, ng.b2, ng.b2);
}
private bool too_big (NG8 ng)
{
// Stop computing if a number reaches the threshold.
return (too_big (ng.a12) || too_big (ng.a1) ||
too_big (ng.a2) || too_big (ng.a) ||
too_big (ng.b12) || too_big (ng.b1) ||
too_big (ng.b2) || too_big (ng.b));
}
private bool too_big (Term u)
{
return (abs (u) >= abs (number_that_is_too_big));
}
private bool treat_as_infinite (Term u)
{
return (abs(u) >= abs (practically_infinite));
}
}
__gshared NG8 ng8_add = new NG8 (0, 1, 1, 0, 0, 0, 0, 1);
__gshared NG8 ng8_sub = new NG8 (0, 1, -1, 0, 0, 0, 0, 1);
__gshared NG8 ng8_mul = new NG8 (1, 0, 0, 0, 0, 0, 0, 1 );
__gshared NG8 ng8_div = new NG8 (0, 1, 0, 0, 0, 0, 1, 0);
//--------------------------------------------------------------------
void
show (string expression, CF cf, string note = "")
{
auto line = rightJustify (expression, 19) ~ " => ";
auto cf_str = to!string (cf);
if (note == "")
line ~= cf_str;
else
line ~= leftJustify (cf_str, 48) ~ note;
writeln (line);
}
int
main (char[][] args)
{
show ("golden ratio", goldenRatio, "(1 + sqrt(5))/2");
show ("silver ratio", silverRatio, "1 + sqrt(2)");
show ("sqrt(2)", sqrt2);
show ("13/11", frac_13_11);
show ("22/7", frac_22_7);
show ("one", one);
show ("two", two);
show ("three", three);
show ("four", four);
show ("(1 + 1/sqrt(2))/2",
new cfNG8 (new NG8 (0, 1, 0, 0, 0, 0, 2, 0),
silverRatio, sqrt2),
"method 1");
show ("(1 + 1/sqrt(2))/2",
new cfNG8 (new NG8 (1, 0, 0, 1, 0, 0, 0, 8),
silverRatio, silverRatio),
"method 2");
show ("(1 + 1/sqrt(2))/2", (one + (one / sqrt2)) / two,
"method 3");
show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2);
show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2);
show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2);
show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2);
return 0;
}
//--------------------------------------------------------------------
- Output:
$ gdc -g -Wall -Wextra bivariate_continued_fraction_task_dlang.d && ./a.out
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (1 + sqrt(5))/2
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] 1 + sqrt(2)
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
one => [1]
two => [2]
three => [3]
four => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
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)
Python
#!/bin/env python3
#---------------------------------------------------------------------
class ContinuedFraction:
def __init__ (self, func, env = None):
self._terminated = False # Are there no more terms?
self._memo = [] # The memoized terms.
self._func = func # A function that produces terms.
self._env = env # An environment for _func.
def __getitem__ (self, i):
while not self._terminated and len(self._memo) <= i:
term = self._func (i, self._env) # i may be ignored.
if term is None:
self._terminated = True
else:
self._memo.append (term)
return (self._memo[i] if i < len(self._memo) else None)
class NG8:
def __init__ (self, a12, a1, a2, a, b12, b1, b2, b):
self.a12 = a12
self.a1 = a1
self.a2 = a2
self.a = a
self.b12 = b12
self.b1 = b1
self.b2 = b2
self.b = b
def coefs (self):
return (self.a12, self.a1, self.a2, self.a,
self.b12, self.b1, self.b2, self.b)
def __call__ (self, x, y):
return apply_ng8 (self.coefs(), x, y)
#---------------------------------------------------------------------
def cf2string (cf, max_terms = 20):
i = 0
s = "["
done = False
while not done:
term = cf[i]
if term is None:
s += "]"
done = True
elif i == max_terms:
s += ",...]"
done = True
else:
if i == 0:
separator = ""
elif i == 1:
separator = ";"
else:
separator = ","
s += separator
s += str (term)
i += 1
return s
def r2cf (n, d): # Rational to continued fraction.
def func (i, env):
n = env[0]
d = env[1]
if d == 0 :
retval = None
else:
q = n // d
r = n - (q * d)
env[0] = d
env[1] = r
retval = q
return retval
return ContinuedFraction (func, [n, d])
def i2cf (i): # Integer to continued fraction.
return r2cf (i, 1)
def apply_ng8 (ng8_tuple, x, y):
# Thresholds chosen merely for demonstration.
number_that_is_too_big = 2 ** 512
practically_infinite = 2 ** 64
ng = 0
ix = 1
iy = 2
overflow = 3
def too_big (values):
# Stop computing if a number reaches the threshold.
return any ((abs (x) >= abs (number_that_is_too_big)
for x in values))
def treat_as_infinite (x):
return (abs (x) >= abs (practically_infinite))
def func (i, env):
def absorb_x_term ():
(a12, a1, a2, a, b12, b1, b2, b) = env[ng]
if env[overflow]:
term = None
else:
term = x[env[ix]]
env[ix] += 1
if term is not None:
new_ng = (a2 + (a12 * term), a + (a1 * term), a12, a1,
b2 + (b12 * term), b + (b1 * term), b12, b1)
if not too_big (new_ng):
env[ng] = new_ng
else:
env[ng] = (a12, a1, a12, a1, b12, b1, b12, b1)
env[overflow] = True
else:
env[ng] = (a12, a1, a12, a1, b12, b1, b12, b1)
return
def absorb_y_term ():
(a12, a1, a2, a, b12, b1, b2, b) = env[ng]
if env[overflow]:
term = None
else:
term = y[env[iy]]
env[iy] += 1
if term is not None:
new_ng = (a1 + (a12 * term), a12, a + (a2 * term), a2,
b1 + (b12 * term), b12, b + (b2 * term), b2)
if not too_big (new_ng):
env[ng] = new_ng
else:
env[ng] = (a12, a12, a2, a2, b12, b12, b2, b2)
env[overflow] = True
else:
env[ng] = (a12, a12, a2, a2, b12, b12, b2, b2)
return
done = False
while not done:
(a12, a1, a2, a, b12, b1, b2, b) = env[ng]
if b == 0 and b1 == 0 and b2 == 0 and b12 == 0:
# There are no more terms.
retval = None
done = True
elif b == 0 and b2 == 0:
absorb_x_term ()
elif b == 0 or b2 == 0:
absorb_y_term ()
elif b1 == 0:
absorb_x_term ()
else:
q = a // b
q1 = a1 // b1
q2 = a2 // b2
if b12 != 0:
q12 = a12 // b12
if b12 != 0 and q == q1 and q == q2 and q == q12:
# Output a term. Notice the resemblance to r2cf.
env[0] = (b12, b1, b2, b, # Divisors.
a12 - (b12 * q), # Remainder.
a1 - (b1 * q), # Remainder.
a2 - (b2 * q), # Remainder.
a - (b * q)) # Remainder.
retval = (None if treat_as_infinite (q) else q)
done = True
else:
# Rather than compare fractions, we will put the
# numerators over a common denominator of b*b1*b2,
# and then compare the new numerators.
n = a * b1 * b2
n1 = a1 * b * b2
n2 = a2 * b * b1
if abs (n1 - n) > abs (n2 - n):
absorb_x_term ()
else:
absorb_y_term ()
return retval
return ContinuedFraction (func, [ng8_tuple, 0, 0, False])
#---------------------------------------------------------------------
golden_ratio = ContinuedFraction (lambda i, env : 1) # (1 + sqrt(5))/2
silver_ratio = ContinuedFraction (lambda i, env : 2) # 1 + sqrt(2)
sqrt2 = ContinuedFraction (lambda i, env : 1 if i == 0 else 2)
frac_13_11 = r2cf (13, 11)
frac_22_7 = r2cf (22, 7)
one = i2cf (1)
two = i2cf (2)
three = i2cf (3)
four = i2cf (4)
cf_add = NG8 (0, 1, 1, 0, 0, 0, 0, 1)
cf_sub = NG8 (0, 1, -1, 0, 0, 0, 0, 1)
cf_mul = NG8 (1, 0, 0, 0, 0, 0, 0, 1)
cf_div = NG8 (0, 1, 0, 0, 0, 0, 1, 0)
print (" golden ratio => ", cf2string (golden_ratio))
print (" silver ratio => ", cf2string (silver_ratio))
print (" sqrt(2) => ", cf2string (sqrt2))
print (" 13/11 => ", cf2string (frac_13_11))
print (" 22/7 => ", cf2string (frac_22_7))
print (" 1 => ", cf2string (one))
print (" 2 => ", cf2string (two))
print (" 3 => ", cf2string (three))
print (" 4 => ", cf2string (four))
print (" (1 + 1/sqrt(2))/2 => ",
cf2string (cf_div (cf_add (one, cf_div (one, sqrt2)), two)),
" method 1")
print (" (1 + 1/sqrt(2))/2 => ",
cf2string (NG8 (1, 0, 0, 1, 0, 0, 0, 8) (silver_ratio,
silver_ratio)),
" method 2")
print (" (1 + 1/sqrt(2))/2 => ",
cf2string (cf_div (cf_div (cf_div (silver_ratio, sqrt2),
sqrt2),
sqrt2)),
" method 3")
print (" sqrt(2) + sqrt(2) => ", cf2string (cf_add (sqrt2, sqrt2)))
print (" sqrt(2) - sqrt(2) => ", cf2string (cf_sub (sqrt2, sqrt2)))
print (" sqrt(2) * sqrt(2) => ", cf2string (cf_mul (sqrt2, sqrt2)))
print (" sqrt(2) / sqrt(2) => ", cf2string (cf_div (sqrt2, sqrt2)))
#---------------------------------------------------------------------
- Output:
$ python3 bivariate_continued_fraction_task.py
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
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
Scheme
The program is for R7RS Scheme. To compile it with CHICKEN, you need the r7rs and srfi-41 eggs.
(cond-expand
(r7rs)
(chicken (import (r7rs))))
;;;-------------------------------------------------------------------
(define-library (cf)
(export make-cf
cf?
*cf-max-terms*
cf->string
number->cf
ng8->procedure)
(import (scheme base)
(scheme case-lambda))
(begin
(define-record-type <cf>
(%%make-cf terminated? ;; No more terms?
m ;; No. of terms memoized.
memo ;; Memoized terms.
gen) ;; Term-generating thunk.
cf?
(terminated? cf-terminated? set-cf-terminated?!)
(m cf-m set-cf-m!)
(memo cf-memo set-cf-memo!)
(gen cf-gen set-cf-gen!))
(define (make-cf gen) ; Make a new continued fraction.
(%%make-cf #f 0 (make-vector 32) gen))
(define (cf-ref cf i) ; Get the ith term, or #f if there is none.
(define (get-more-terms! needed)
(do () ((or (cf-terminated? cf) (= (cf-m cf) needed)))
(let ((term ((cf-gen cf))))
(if term
(begin
(vector-set! (cf-memo cf) (cf-m cf) term)
(set-cf-m! cf (+ (cf-m cf) 1)))
(set-cf-terminated?! cf #t)))))
(define (update! needed)
(cond ((cf-terminated? cf) (begin))
((<= needed (cf-m cf)) (begin))
((<= needed (vector-length (cf-memo cf)))
(get-more-terms! needed))
(else ;; Increase the storage space for memoization.
(let* ((n1 (+ needed needed))
(memo1 (make-vector n1)))
(vector-copy! memo1 0 (cf-memo cf) 0 (cf-m cf))
(set-cf-memo! cf memo1))
(get-more-terms! needed))))
(update! (+ i 1))
(and (< i (cf-m cf))
(vector-ref (cf-memo cf) i)))
(define *cf-max-terms* ; Default term-count limit for cf->string.
(make-parameter 20))
(define cf->string ; Make a string for a continued fraction.
(case-lambda
((cf) (cf->string cf (*cf-max-terms*)))
((cf max-terms)
(let loop ((i 0)
(s "["))
(let ((term (cf-ref cf i)))
(cond ((not term) (string-append s "]"))
((= i max-terms) (string-append s ",...]"))
(else
(let ((separator (case i
((0) "")
((1) ";")
(else ",")))
(term-str (number->string term)))
(loop (+ i 1) (string-append s separator
term-str))))))))))
(define (number->cf num) ; Convert a number to a continued fraction.
(let ((num (exact num)))
(let ((n (numerator num))
(d (denominator num)))
(make-cf
(lambda ()
(and (not (zero? d))
(let-values (((q r) (floor/ n d)))
(set! n d)
(set! d r)
q)))))))
(define (%%divide a b)
(if (zero? b)
(values #f #f)
(floor/ a b)))
(define (ng8->procedure ng8)
;; Thresholds chosen merely for demonstration.
(define number-that-is-too-big (expt 2 512))
(define practically-infinite (expt 2 64))
(define too-big? ; Stop computing if a no. reaches a threshold.
(lambda (values)
(cond ((null? values) #f)
((>= (abs (car values))
(abs number-that-is-too-big)) #t)
(else (too-big? (cdr values))))))
(define (treat-as-infinite? term)
(>= (abs term) (abs practically-infinite)))
(lambda (x y)
(define ng ng8)
(define xsource (lambda (i) (cf-ref x i)))
(define ysource (lambda (i) (cf-ref y i)))
(define ix 0)
(define iy 0)
;; The procedures "main", "compare-quotients",
;; "absorb-x-term", and "absorb-y-term" form a mutually
;; tail-recursive set. In standard Scheme, such an arrangement
;; requires no special notations, and WILL NOT blow up the
;; stack.
(define (main) ; "main" is the term-generating thunk.
(define-values (a12 a1 a2 a b12 b1 b2 b)
(apply values ng))
(define bz? (zero? b))
(define b1z? (zero? b1))
(define b2z? (zero? b2))
(define b12z? (zero? b12))
(cond
((and bz? b1z? b2z? b12z?) #f)
((and bz? b2z?) (absorb-x-term))
((or bz? b2z?) (absorb-y-term))
(b1z? (absorb-x-term))
(else
(let-values (((q r) (%%divide a b))
((q1 r1) (%%divide a1 b1))
((q2 r2) (%%divide a2 b2))
((q12 r12) (%%divide a12 b12)))
(if (and (not b12z?) (= q q1 q2 q12))
(output-a-term q b12 b1 b2 b r12 r1 r2 r)
(compare-quotients a2 a1 a b2 b1 b))))))
(define (compare-quotients a2 a1 a b2 b1 b)
(let ((n (* a b1 b2))
(n1 (* a1 b b2))
(n2 (* a2 b b1)))
(if (> (abs (- n1 n)) (abs (- n2 n)))
(absorb-x-term)
(absorb-y-term))))
(define (absorb-x-term)
(define-values (a12 a1 a2 a b12 b1 b2 b)
(apply values ng))
(define term (xsource ix))
(set! ix (+ ix 1))
(if term
(let ((new-ng (list (+ a2 (* a12 term))
(+ a (* a1 term))
a12 a1
(+ b2 (* b12 term))
(+ b (* b1 term))
b12 b1)))
(if (not (too-big? new-ng))
(set! ng new-ng)
(begin
;; Replace the x source with one that returns no
;; terms.
(set! xsource (lambda (i) #f))
(set! ng (list a12 a1 a12 a1 b12 b1 b12 b1)))))
(set! ng (list a12 a1 a12 a1 b12 b1 b12 b1)))
(main))
(define (absorb-y-term)
(define-values (a12 a1 a2 a b12 b1 b2 b)
(apply values ng))
(define term (ysource iy))
(set! iy (+ iy 1))
(if term
(let ((new-ng (list (+ a1 (* a12 term)) a12
(+ a (* a2 term)) a2
(+ b1 (* b12 term)) b12
(+ b (* b2 term)) b2)))
(if (not (too-big? new-ng))
(set! ng new-ng)
(begin
;; Replace the y source with one that returns no
;; terms.
(set! ysource (lambda (i) #f))
(set! ng (list a12 a12 a2 a2 b12 b12 b2 b2)))))
(set! ng (list a12 a12 a2 a2 b12 b12 b2 b2)))
(main))
(define (output-a-term q b12 b1 b2 b r12 r1 r2 r)
(let ((new-ng (list b12 b1 b2 b r12 r1 r2 r))
(term (and (not (treat-as-infinite? q)) q)))
(set! ng new-ng)
term))
(make-cf main))) ;; end procedure ng8->procedure
)) ;; end library (cf)
;;;-------------------------------------------------------------------
(import (scheme base))
(import (scheme case-lambda))
(import (scheme write))
(import (cf))
(define golden-ratio (make-cf (lambda () 1)))
(define silver-ratio (make-cf (lambda () 2)))
(define sqrt2 (make-cf (let ((next-term 1))
(lambda ()
(let ((term next-term))
(set! next-term 2)
term)))))
(define frac13/11 (number->cf 13/11))
(define frac22/7 (number->cf 22/7))
(define one (number->cf 1))
(define two (number->cf 2))
(define three (number->cf 3))
(define four (number->cf 4))
(define cf+ (ng8->procedure '(0 1 1 0 0 0 0 1)))
(define cf- (ng8->procedure '(0 1 -1 0 0 0 0 1)))
(define cf* (ng8->procedure '(1 0 0 0 0 0 0 1)))
(define cf/ (ng8->procedure '(0 1 0 0 0 0 1 0)))
(define show
(case-lambda
((expression cf note)
(display expression)
(display " => ")
(display (cf->string cf))
(display note)
(newline))
((expression cf)
(show expression cf ""))))
(show " golden ratio" golden-ratio)
(show " silver ratio" silver-ratio)
(show " sqrt(2)" sqrt2)
(show " 13/11" frac13/11)
(show " 22/7" frac22/7)
(show " 1" one)
(show " 2" two)
(show " 3" three)
(show " 4" four)
(show " (1 + 1/sqrt(2))/2" (cf/ (cf+ one (cf/ one sqrt2)) two)
" method 1")
(show " (1 + 1/sqrt(2))/2" ((ng8->procedure '(1 0 0 1 0 0 0 8))
silver-ratio silver-ratio)
" method 2")
(show " (1 + 1/sqrt(2))/2" (cf/ (cf/ (cf/ silver-ratio sqrt2)
sqrt2)
sqrt2)
" method 3")
(show " sqrt(2) + sqrt(2)" (cf+ sqrt2 sqrt2))
(show " sqrt(2) - sqrt(2)" (cf- sqrt2 sqrt2))
(show " sqrt(2) * sqrt(2)" (cf* sqrt2 sqrt2))
(show " sqrt(2) / sqrt(2)" (cf/ sqrt2 sqrt2))
;;;-------------------------------------------------------------------
- Output:
$ gosh bivariate-continued-fraction-task-frompy.scm
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...]
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
13/11 => [1;5,2]
22/7 => [3;7]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 1
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 2
(1 + 1/sqrt(2))/2 => [0;1,5,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] method 3
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt(2) - sqrt(2) => [0]
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
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