Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2): Difference between revisions
Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2) (view source)
Revision as of 17:20, 20 November 2023
, 5 months ago→{{header|Wren}}: Changed to Wren S/H
m (→{{header|Wren}}: Changed to Wren S/H) |
|||
| (39 intermediate revisions by 3 users not shown) | |||
Line 1:
{{
This task performs the basic mathematical functions on 2 continued fractions. This requires the full version of matrix NG:
: <math>\begin{bmatrix}
Line 59:
When performing arithmetic operation on two potentially infinite continued fractions it is possible to generate a rational number. eg <math>\sqrt{2}</math> * <math>\sqrt{2}</math> 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.
=={{header|Ada}}==
{{trans|Python}}
{{works with|GCC|12.2.1}}
<syntaxhighlight lang="ada">
pragma ada_2022; -- When big_integers were introduced.
with ada.numerics.big_numbers.big_integers;
use ada.numerics.big_numbers.big_integers;
with ada.strings; use ada.strings;
with ada.strings.fixed; use ada.strings.fixed;
with ada.strings.unbounded; use ada.strings.unbounded;
with ada.text_io; use ada.text_io;
procedure BIVARIATE_CONTINUED_FRACTION_TASK is
package CONTINUED_FRACTIONS is
type memoization_storage is array (natural range <>) of big_integer;
type memoization_access is access memoization_storage;
type continued_fraction_record is abstract tagged
record
terminated : boolean := false; -- Are there no more terms?
memo_count : natural := 0; -- How many terms are memoized?
memo : memoization_access -- Memoized terms.
:= new memoization_storage (0 .. 31);
end record;
procedure generate_term (cf : in out continued_fraction_record;
output_exists : out boolean;
output : out big_integer) is abstract;
type continued_fraction is access all
continued_fraction_record'class; -- The 'class notation is important.
function term_exists (cf : in continued_fraction;
i : in natural)
return boolean;
function get_term (cf : in continued_fraction;
i : in natural)
return big_integer
with pre => i < cf.memo_count;
function cf2string (cf : in continued_fraction;
max_terms : in positive := 20)
return unbounded_string;
end CONTINUED_FRACTIONS;
package body CONTINUED_FRACTIONS is
function term_exists (cf : in continued_fraction;
i : in natural)
return boolean is
procedure resize_if_necessary is
memo1 : memoization_access;
begin
if cf.memo'length <= i then
memo1 := new memoization_storage(0 .. 2 * (i + 1));
for i in 0 .. cf.memo_count - 1 loop
memo1(i) := cf.memo(i);
end loop;
cf.memo := memo1;
end if;
end;
exists : boolean;
term : big_integer;
begin
if i < cf.memo_count then
exists := true;
elsif cf.terminated then
exists := false;
else
resize_if_necessary;
while cf.memo_count <= i and not cf.terminated loop
generate_term (cf.all, exists, term);
if exists then
cf.memo(cf.memo_count) := term;
cf.memo_count := cf.memo_count + 1;
else
cf.terminated := true;
end if;
end loop;
exists := term_exists (cf, i);
end if;
return exists;
end;
function get_term (cf : in continued_fraction;
i : in natural)
return big_integer is
begin
return cf.memo(i);
end;
function cf2string (cf : in continued_fraction;
max_terms : in positive := 20)
return unbounded_string is
s : unbounded_string := null_unbounded_string;
done : boolean;
i : natural;
term : big_integer;
begin
s := s & "[";
i := 0;
done := false;
while not done loop
if not term_exists (cf, i) then
s := s & "]";
done := true;
elsif i = max_terms then
s := s & ",...]";
done := true;
else
if i = 1 then
s := s & ";";
elsif i /= 0 then
s := s & ",";
end if;
term := get_term (cf, i);
s := s & trim (term'image, left);
i := i + 1;
end if;
end loop;
return s;
end;
end CONTINUED_FRACTIONS;
package CONSTANT_TERM_CONTINUED_FRACTIONS is
use CONTINUED_FRACTIONS;
type constant_term_continued_fraction_record is
new continued_fraction_record with
record
term : big_integer;
end record;
type constant_term_continued_fraction is access all
constant_term_continued_fraction_record;
function constant_term_cf (term : in big_integer)
return continued_fraction;
overriding
procedure generate_term (cf : in out constant_term_continued_fraction_record;
output_exists : out boolean;
output : out big_integer);
end CONSTANT_TERM_CONTINUED_FRACTIONS;
package body CONSTANT_TERM_CONTINUED_FRACTIONS is
function constant_term_cf (term : in big_integer)
return continued_fraction is
cf : constant_term_continued_fraction;
begin
cf := new constant_term_continued_fraction_record;
cf.term := term;
return continued_fraction (cf);
end;
overriding
procedure generate_term (cf : in out constant_term_continued_fraction_record;
output_exists : out boolean;
output : out big_integer) is
begin
output_exists := true;
output := cf.term;
end;
end CONSTANT_TERM_CONTINUED_FRACTIONS;
package INTEGER_CONTINUED_FRACTIONS is
use CONTINUED_FRACTIONS;
type integer_continued_fraction_record is
new continued_fraction_record with
record
term : big_integer;
done : boolean := false;
end record;
type integer_continued_fraction is access all
integer_continued_fraction_record;
function i2cf (term : in big_integer)
return continued_fraction;
overriding
procedure generate_term (cf : in out integer_continued_fraction_record;
output_exists : out boolean;
output : out big_integer);
end INTEGER_CONTINUED_FRACTIONS;
package body INTEGER_CONTINUED_FRACTIONS is
function i2cf (term : in big_integer)
return continued_fraction is
cf : integer_continued_fraction;
begin
cf := new integer_continued_fraction_record;
cf.term := term;
return continued_fraction (cf);
end;
overriding
procedure generate_term (cf : in out integer_continued_fraction_record;
output_exists : out boolean;
output : out big_integer) is
begin
output_exists := not (cf.done);
cf.done := true;
if output_exists then
output := cf.term;
end if;
end;
end INTEGER_CONTINUED_FRACTIONS;
package NG8_CONTINUED_FRACTIONS is
use CONTINUED_FRACTIONS;
stopping_processing_threshold : big_integer := 2 ** 512;
infinitization_threshold : big_integer := 2 ** 64;
type ng8_continued_fraction_record is
new continued_fraction_record with
record
a12, a1, a2, a : big_integer;
b12, b1, b2, b : big_integer;
x, y : continued_fraction;
ix, iy : natural;
xoverflow : boolean;
yoverflow : boolean;
end record;
type ng8_continued_fraction is access all
ng8_continued_fraction_record;
function apply_ng8 (a12, a1, a2, a : in big_integer;
b12, b1, b2, b : in big_integer;
x, y : in continued_fraction)
return continued_fraction;
-- Addition.
function "+" (x, y : in continued_fraction)
return continued_fraction;
function "+" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "+" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;
-- Keeping the same sign. (Effectively clones x as an
-- ng8_continued_fraction.)
function "+" (x : in continued_fraction)
return continued_fraction;
-- Subtraction.
function "-" (x, y : in continued_fraction)
return continued_fraction;
function "-" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "-" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;
-- Negation.
function "-" (x : in continued_fraction)
return continued_fraction;
-- Multiplication.
function "*" (x, y : in continued_fraction)
return continued_fraction;
function "*" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "*" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;
-- Division.
function "/" (x, y : in continued_fraction)
return continued_fraction;
function "/" (x : in continued_fraction;
y : in big_integer)
return continued_fraction;
function "/" (x : in big_integer;
y : in continued_fraction)
return continued_fraction;
-- A rational number as a continued fraction. The terms are
-- memoized, so this implementation will not be as inefficient as
-- one might suppose.
function r2cf (n, d : in big_integer)
return continued_fraction;
overriding
procedure generate_term (cf : in out ng8_continued_fraction_record;
output_exists : out boolean;
output : out big_integer);
end NG8_CONTINUED_FRACTIONS;
package body NG8_CONTINUED_FRACTIONS is
use CONTINUED_FRACTIONS;
use CONSTANT_TERM_CONTINUED_FRACTIONS;
-- An arbitrary infinite source of non-zero finite terms.
forever_cf : continued_fraction := constant_term_cf (1234);
function apply_ng8 (a12, a1, a2, a : in big_integer;
b12, b1, b2, b : in big_integer;
x, y : in continued_fraction)
return continued_fraction is
cf : ng8_continued_fraction;
begin
cf := new ng8_continued_fraction_record;
cf.a12 := a12;
cf.a1 := a1;
cf.a2 := a2;
cf.a := a;
cf.b12 := b12;
cf.b1 := b1;
cf.b2 := b2;
cf.b := b;
cf.x := x;
cf.y := y;
cf.ix := 0;
cf.iy := 0;
cf.xoverflow := false;
cf.yoverflow := false;
return continued_fraction (cf);
end;
function "+" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1, x, y);
end;
function "+" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, y, 0, 0, 0, 1, x, forever_cf);
end;
function "+" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 1, x, 0, 0, 0, 1, forever_cf, y);
end;
function "+" (x : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 1, 0, 0, 0, 0, 1, forever_cf, x);
end;
function "-" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1, x, y);
end;
function "-" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, -y, 0, 0, 0, 1, x, forever_cf);
end;
function "-" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, -1, x, 0, 0, 0, 1, forever_cf, y);
end;
function "-" (x : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, -1, 0, 0, 0, 0, 1, forever_cf, x);
end;
function "*" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1, x, y);
end;
function "*" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, y, 0, 0, 0, 0, 0, 1, x, forever_cf);
end;
function "*" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, x, 0, 0, 0, 0, 1, forever_cf, y);
end;
function "/" (x, y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0, x, y);
end;
function "/" (x : in continued_fraction;
y : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 1, 0, 0, 0, 0, 0, y, x, forever_cf);
end;
function "/" (x : in big_integer;
y : in continued_fraction)
return continued_fraction is
begin
return apply_ng8 (0, 0, 0, x, 0, 0, 1, 0, forever_cf, y);
end;
function r2cf (n, d : in big_integer)
return continued_fraction is
begin
return apply_ng8 (0, 0, 0, n, 0, 0, 0, d, forever_cf, forever_cf);
end;
procedure possibly_infinitize_output (q : in big_integer;
output_exists : out boolean;
output : out big_integer) is
begin
output_exists := abs (q) < abs (infinitization_threshold);
if output_exists then
output := q;
end if;
end;
procedure divide (a, b : in big_integer;
q, r : out big_integer) is
begin
if b /= 0 then
q := a / b;
r := a rem b;
end if;
end;
function too_big (num : big_integer)
return boolean is
begin
return (abs (stopping_processing_threshold) <= abs (num));
end;
function any_too_big (a, b, c, d, e, f, g, h : in big_integer)
return boolean is
begin
return (too_big (a) or else
too_big (b) or else
too_big (c) or else
too_big (d) or else
too_big (e) or else
too_big (f) or else
too_big (g) or else
too_big (h));
end;
overriding
procedure generate_term (cf : in out ng8_continued_fraction_record;
output_exists : out boolean;
output : out big_integer) is
a12, a1, a2, a : big_integer;
b12, b1, b2, b : big_integer;
q12, q1, q2, q : big_integer;
r12, r1, r2, r : big_integer;
absorb_y_instead_of_x : boolean;
done : boolean;
function all_b_are_zero
return boolean is
begin
return (b12 = 0 and b1 = 0 and b2 = 0 and b = 0);
end;
function all_q_are_equal
return boolean is
begin
return (q = q1 and q = q2 and q = q12);
end;
procedure compare_fractions is
n1, n2, n : big_integer;
begin
-- Rather than compare fractions, we will put the numerators over
-- a common denominator of b*b1*b2, and then compare the new
-- numerators.
n1 := a1 * b2 * b;
n2 := a2 * b1 * b;
n := a * b1 * b2;
absorb_y_instead_of_x := (abs (n1 - n) <= abs (n2 - n));
end;
procedure absorb_x_term is
term : big_integer;
new_a12, new_a1, new_a2, new_a : big_integer;
new_b12, new_b1, new_b2, new_b : big_integer;
begin
new_a2 := a12;
new_a := a1;
new_b2 := b12;
new_b := b1;
if not cf.xoverflow and then term_exists (cf.x, cf.ix) then
term := get_term (cf.x, cf.ix);
new_a12 := a2 + (a12 * term);
new_a1 := a + (a1 * term);
new_b12 := b2 + (b12 * term);
new_b1 := b + (b1 * term);
if any_too_big (new_a12, new_a1, new_a2, new_a,
new_b12, new_b1, new_b2, new_b) then
cf.xoverflow := true;
new_a12 := a12;
new_a1 := a1;
new_b12 := b12;
new_b1 := b1;
else
cf.ix := cf.ix + 1;
end if;
else
new_a12 := a12;
new_a1 := a1;
new_b12 := b12;
new_b1 := b1;
end if;
a12 := new_a12;
a1 := new_a1;
a2 := new_a2;
a := new_a;
b12 := new_b12;
b1 := new_b1;
b2 := new_b2;
b := new_b;
end;
procedure absorb_y_term is
term : big_integer;
new_a12, new_a1, new_a2, new_a : big_integer;
new_b12, new_b1, new_b2, new_b : big_integer;
begin
new_a1 := a12;
new_a := a2;
new_b1 := b12;
new_b := b2;
if not cf.yoverflow and then term_exists (cf.y, cf.iy) then
term := get_term (cf.y, cf.iy);
new_a12 := a1 + (a12 * term);
new_a2 := a + (a2 * term);
new_b12 := b1 + (b12 * term);
new_b2 := b + (b2 * term);
if any_too_big (new_a12, new_a1, new_a2, new_a,
new_b12, new_b1, new_b2, new_b) then
cf.yoverflow := true;
new_a12 := a12;
new_a2 := a2;
new_b12 := b12;
new_b2 := b2;
else
cf.iy := cf.iy + 1;
end if;
else
new_a12 := a12;
new_a2 := a2;
new_b12 := b12;
new_b2 := b2;
end if;
a12 := new_a12;
a1 := new_a1;
a2 := new_a2;
a := new_a;
b12 := new_b12;
b1 := new_b1;
b2 := new_b2;
b := new_b;
end;
procedure absorb_term is
begin
if absorb_y_instead_of_x then
absorb_y_term;
else
absorb_x_term;
end if;
end;
begin
a12 := cf.a12;
a1 := cf.a1;
a2 := cf.a2;
a := cf.a;
b12 := cf.b12;
b1 := cf.b1;
b2 := cf.b2;
b := cf.b;
done := false;
while not done loop
absorb_y_instead_of_x := false;
if all_b_are_zero then
-- There are no more terms.
output_exists := false;
done := true;
elsif b2 = 0 and b = 0 then
null;
elsif b2 = 0 or b = 0 then
absorb_y_instead_of_x := true;
elsif b1 = 0 then
null;
else
divide (a12, b12, q12, r12);
divide (a1, b1, q1, r1);
divide (a2, b2, q2, r2);
divide (a, b, q, r);
if b12 /= 0 and then all_q_are_equal then
-- Output a term.
cf.a12 := b12;
cf.a1 := b1;
cf.a2 := b2;
cf.a := b;
cf.b12 := r12;
cf.b1 := r1;
cf.b2 := r2;
cf.b := r;
possibly_infinitize_output (q, output_exists, output);
done := true;
else
compare_fractions;
end if;
end if;
absorb_term;
end loop;
end;
end NG8_CONTINUED_FRACTIONS;
use CONTINUED_FRACTIONS;
use CONSTANT_TERM_CONTINUED_FRACTIONS;
use INTEGER_CONTINUED_FRACTIONS;
use NG8_CONTINUED_FRACTIONS;
procedure show (expression : in string;
cf : in continued_fraction;
note : in string := "") is
expr : string := 19 * ' ';
contfrac : string := 48 * ' ';
begin
move (source => expression,
target => expr,
justify => right);
put (expr);
put (" => ");
if note = "" then
put_line (to_string (cf2string (cf)));
else
move (source => to_string (cf2string (cf)),
target => contfrac,
justify => left);
put (contfrac);
put_line (note);
end if;
end;
golden_ratio : continued_fraction := constant_term_cf (1);
silver_ratio : continued_fraction := constant_term_cf (2);
one : continued_fraction := i2cf (1);
two : continued_fraction := i2cf (2);
three : continued_fraction := i2cf (3);
four : continued_fraction := i2cf (4);
sqrt2 : continued_fraction := silver_ratio - 1;
begin
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "(1 + sqrt(2))");
show ("sqrt(2)", sqrt2, "silver ratio minus 1");
show ("13/11", r2cf (13, 11));
show ("22/7", r2cf (22, 7), "approximately pi");
show ("1", one);
show ("2", two);
show ("3", three);
show ("4", four);
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (0, 1, 0, 0, 0, 0, 2, 0, silver_ratio, sqrt2),
"method 1");
show ("(1 + 1/sqrt(2))/2",
apply_ng8 (1, 0, 0, 1, 0, 0, 0, 8, silver_ratio, silver_ratio),
"method 2");
show ("(1 + 1/sqrt(2))/2",
(one / 2) * (one + (1 / sqrt2)),
"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);
end BIVARIATE_CONTINUED_FRACTION_TASK;
-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end:
</syntaxhighlight>
{{out}}
<pre>$ gnatmake -q -g bivariate_continued_fraction_task.adb && ./bivariate_continued_fraction_task
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,...] silver ratio minus 1
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
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]
</pre>
=={{header|ATS}}==
{{trans|Python}}
===Using 128-bit integers===
(''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.)
Line 1,069 ⟶ 1,819:
sqrt(2) * sqrt(2) => [2;7530688524100]
sqrt(2) / sqrt(2) => [1]
</pre>
===Using multiple precision numbers===
{{trans|Standard ML}}
{{libheader|ats2-xprelude}}
{{libheader|GMP}}
For this program you need [https://sourceforge.net/p/chemoelectric/ats2-xprelude ats2-xprelude].
The program closely follows the Standard ML code, so one can compare the two languages with each other. They have similar syntaxes, but are very different. Notice, for instance, that ATS has overloads, whereas SML does not. (SML has signatures with respective namespaces.) In ATS, a function is not a closure unless you explicitly make it one, whereas in SML no special notation is needed. And so on.
ATS is translated to C, and its functions (except closures) are essentially just C functions. One can write Arduino code and Linux kernel modules in ATS, because ATS is, in some sense, an elaborate way to write C. Nevertheless, there is enough similarity between ATS and Standard ML to easily translate the SML code for this Rosetta Code task to ATS.
I have broken the program into three files, to demonstrate what an ATS program might look like, if it were broken into separately compiled parts.
The first file is an "interface" specification for a <code>continued_fraction</code> data type. The file is called <code>continued_fraction.sats</code>:
<syntaxhighlight lang="ats">
(* "Static" file. (Exported declarations.) *)
(* To set up a predictable name-mangling scheme: *)
#define ATS_PACKNAME "rosetta-code.continued_fraction"
(* Load declarations from ats2-xprelude: *)
#include "xprelude/HATS/xprelude_sats.hats"
staload "xprelude/SATS/exrat.sats"
(* A term_generator thunk generates terms, which a continued_fraction
data structure memoizes. The internals of continued_fraction are
not exposed here. It is an abstract type, the size of (but not the
same type as) a pointer. SIDE NOTE: In ATS2, we get a conventional
function, rather than a closure, unless we say explicitly that we
want a closure; "cloref1" means a particular kind of closure one
often uses when linking the program with Boehm GC. *)
typedef term_generator = () -<cloref1> Option exrat
abstype continued_fraction = ptr
(* Create a continued fraction. *)
fn continued_fraction_make : term_generator -> continued_fraction
(* Does the indexed term exist? *)
fn term_exists : (continued_fraction, intGte 0) -> bool
(* Retrieve the indexed term. Raise an exception if there is no such
term. The precedence of the overload must exceed that of an
overload that is in the prelude. (To see what I mean, try removing
the "of 1".) *)
val get_term_exn : (continued_fraction, intGte 0) -> exrat
overload [] with get_term_exn of 1
(* Use a continued_fraction as a term_generator thunk. *)
fn continued_fraction_to_term_generator :
continued_fraction -> term_generator
(* Get a human-readable string. *)
val default_max_terms : ref (intGte 1)
fn continued_fraction_to_string_given_max_terms :
(continued_fraction, intGte 1) -> string
fn continued_fraction_to_string_default_max_terms :
continued_fraction -> string
overload continued_fraction_to_string with
continued_fraction_to_string_given_max_terms
overload continued_fraction_to_string with
continued_fraction_to_string_default_max_terms
overload cf2string with continued_fraction_to_string
(* Make a continued_fraction for an integer. *)
fn int_to_continued_fraction : int -> continued_fraction
overload i2cf with int_to_continued_fraction
(* Make a continued_fraction for a rational number. *)
fn exrat_to_continued_fraction : exrat -> continued_fraction
fn rational_to_continued_fraction :
(int, [d : int | d != 0] int d) -> continued_fraction
overload r2cf with exrat_to_continued_fraction
overload r2cf with rational_to_continued_fraction
(* Make a continued_fraction with one term repeated forever. *)
fn continued_fraction_make_constant_term : int -> continued_fraction
overload constant_term_cf with continued_fraction_make_constant_term
(* Make a continued fraction via binary arithmetic operations. (I have
not bothered here to implement ng4, although one likely would wish
to have ng4 as well.) *)
(* The @() denotes an unboxed tuple. A boxed tuple is written '() and
would be put in the heap. An unboxed tuple may also be written
without the @-sign, but then the compiler might confuse it with,
for instance, an argument list. (ATS2 has conventional argument
lists that are distinct from tuples, and supports
call-by-reference, where an argument is mutable.) *)
typedef ng8 = @(exrat, exrat, exrat, exrat,
exrat, exrat, exrat, exrat)
typedef continued_fraction_binary_op_cloref =
(continued_fraction, continued_fraction) -<cloref1> continued_fraction
(* ng8_make_int takes ONE argument, which is a tuple. *)
val ng8_make_int : @(int, int, int, int, int, int, int, int) -> ng8
val ng8_stopping_processing_threshold : ref exrat
val ng8_infinitization_threshold : ref exrat
val ng8_apply : ng8 -> continued_fraction_binary_op_cloref
val ng8_apply_add : continued_fraction_binary_op_cloref
val ng8_apply_sub : continued_fraction_binary_op_cloref
val ng8_apply_mul : continued_fraction_binary_op_cloref
val ng8_apply_div : continued_fraction_binary_op_cloref
(* The following two are regular functions, not closures. They are
translated by the ATS compiler into ordinary C functions. *)
fn ng8_apply_neg : continued_fraction -> continued_fraction
fn ng8_apply_pow : (continued_fraction, int) -> continued_fraction
overload + with ng8_apply_add
overload - with ng8_apply_sub
overload * with ng8_apply_mul
overload / with ng8_apply_div
overload ~ with ng8_apply_neg
overload ** with ng8_apply_pow
(* Miscellanous continued fractions. *)
val zero : continued_fraction
val one : continued_fraction
val two : continued_fraction
val three : continued_fraction
val four : continued_fraction
//
val one_fourth : continued_fraction
val one_third : continued_fraction
val one_half : continued_fraction
val two_thirds : continued_fraction
val three_fourths : continued_fraction
//
val golden_ratio : continued_fraction
val silver_ratio : continued_fraction
val sqrt2 : continued_fraction
val sqrt5 : continued_fraction
</syntaxhighlight>
The second file is an implementation of the stuff declared in the first file. The second file is called <code>continued_fraction.dats</code>:
<syntaxhighlight lang="ats">
(* "Dynamic" file. (Implementations.) *)
(* To set up a predictable name-mangling scheme: *)
#define ATS_PACKNAME "rosetta-code.continued_fraction"
(* Load templates from the ATS2 prelude: *)
#include "share/atspre_staload.hats"
(* Load declarations and templates from ats2-xprelude: *)
#include "xprelude/HATS/xprelude.hats"
staload "xprelude/SATS/exrat.sats"
staload _ = "xprelude/DATS/exrat.dats"
(* Load the declarations for this package: *)
staload "continued_fraction.sats"
typedef cf_record =
(* A cf_record is an unboxed record, denoted by @{}. A boxed record
would be written '{} and would be placed in the heap. Either way,
it is an immutable record. For a mutable record, we would have to
use vtypedef to make it a LINEAR type. *)
@{
terminated = bool, (* Is the generator exhausted? *)
memo_count = size_t, (* How many terms are memoized? *)
(* An arrszref is an array with runtime bounds checking. An
arrszref is less efficient than an arrayref, but will not force
us to use dependent types for the indices. *)
memo = arrszref exrat, (* Memoized terms. *)
generate = term_generator (* The source of terms. *)
}
(* The actual type of a continued_fraction is a MUTABLE reference to
the (immutable) cf_record. Within this file, we may also call the
type cf_t. *)
typedef cf_t = ref cf_record
assume continued_fraction = cf_t
implement
continued_fraction_make generator =
let
val record : cf_record =
@{
terminated = false,
memo_count = i2sz 0,
memo = arrszref_make_elt<exrat> (i2sz 32, exrat_make (0, 1)),
generate = generator
}
in
ref record
end
(* "fn" means a non-recursive function. A function that might be
recursive is written "fun" (or sometimes "fnx"). Incidentally: it
is common to see the recursions put into nested functions, with the
function a programmer is supposed to call being non-recursive. This
is often a matter of style. (By the way, in a "*.sats" file there
is no distinction between "fn" and "fun" that I know of.) *)
fn
resize_if_necessary (cf : cf_t, i : size_t) : void =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if size memo <= i then
let
val new_size = i2sz 2 * (succ i)
val new_memo =
arrszref_make_elt<exrat> (new_size, exrat_make (0, 1))
val new_record : cf_record =
@{
terminated = terminated,
memo_count = memo_count,
memo = new_memo,
generate = generate
}
var i : size_t (* A C-style automatic variable. *)
in
(* A C-style for-loop. *)
for (i := i2sz 0; i <> memo_count; i := succ i)
new_memo[i] := memo[i];
!cf := new_record
end
end
fn
update_terms (cf : cf_t, i : size_t) : void =
let
fun
loop () : void =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if terminated then
()
else if i < memo_count then
()
else
case generate () of
| None () =>
let
val new_record =
@{
terminated = true,
memo_count = memo_count,
memo = memo,
generate = generate
}
in
!cf := new_record
end
| Some term =>
(* "begin-end" is a synonym for "()". *)
begin
memo[memo_count] := term;
let
val new_record =
@{
terminated = false,
memo_count = succ memo_count,
memo = memo,
generate = generate
}
in
!cf := new_record;
loop ()
end
end
end
in
loop ()
end
implement
term_exists (cf, i) =
let
val i = i2sz i
fun
loop () =
let
val @{
terminated = terminated,
memo_count = memo_count,
memo = memo,
generate = generate
} = !cf
in
if i < memo_count then
true
else if terminated then
false
else
begin
resize_if_necessary (cf, i);
update_terms (cf, i);
loop ()
end
end
in
loop ()
end
implement
get_term_exn (cf, i) =
if i2sz i < (!cf).memo_count then
let
val memo = (!cf).memo
in
memo[i]
end
else
$raise IllegalArgExn "get_term_exn:out_of_bounds"
implement
continued_fraction_to_term_generator cf =
let
val i : ref (intGte 0) = ref 0
in
lam () =<cloref1>
let
val j = !i
in
if term_exists (cf, j) then
begin
!i := succ j;
Some (cf[j])
end
else
None ()
end
end
implement default_max_terms = ref 20
implement
continued_fraction_to_string_given_max_terms (cf, max_terms) =
let
fun
loop (i : intGte 0, accum : string) : string =
if ~term_exists (cf, i) then
(* The return value of string_append is a LINEAR, MUTABLE
strptr, which we cast to a nonlinear, immutable string.
(One could introduce one's own shorthands, though.) *)
strptr2string (string_append (accum, "]"))
else if i = max_terms then
strptr2string (string_append (accum, ",...]"))
else
let
val separator =
if i = 0 then
""
else if i = 1 then
";"
else
","
and term_string = tostring_val<exrat> cf[i]
in
loop (succ i,
strptr2string (string_append (accum, separator,
term_string)))
end
in
loop (0, "[")
end
implement
continued_fraction_to_string_default_max_terms cf =
let
val max_terms = !default_max_terms
in
continued_fraction_to_string_given_max_terms (cf, max_terms)
end
implement
int_to_continued_fraction i =
let
val done : ref bool = ref false
val i = (g0i2f i) : exrat
in
continued_fraction_make
(lam () =<cloref1>
if !done then
None ()
else
begin
!done := true;
Some i
end)
end
implement
exrat_to_continued_fraction num =
let
val done : ref bool = ref false
val num : ref exrat = ref num
in
continued_fraction_make
(lam () =<cloref1>
if !done then
None ()
else
let
val q = floor !num
val r = !num - q
in
if iseqz r then
!done := true
else
!num := reciprocal r;
Some q
end)
end
implement
rational_to_continued_fraction (numer, denom) =
exrat_to_continued_fraction (exrat_make (numer, denom))
implement
continued_fraction_make_constant_term i =
let
val i = (g0i2f i) : exrat
in
continued_fraction_make (lam () =<cloref1> Some i)
end
implement
ng8_make_int tuple =
let
val @(a12, a1, a2, a, b12, b1, b2, b) = tuple
fn f (i : int) : exrat = exrat_make (i, 1)
in
@(f a12, f a1, f a2, f a, f b12, f b1, f b2, f b)
end
implement ng8_stopping_processing_threshold =
ref (exrat_make (2, 1) ** 512)
implement ng8_infinitization_threshold =
ref (exrat_make (2, 1) ** 64)
fn
too_big (term : exrat) : bool =
abs (term) >= abs (!ng8_stopping_processing_threshold)
fn
any_too_big (ng : ng8) : bool =
(* "orelse" may also be (and usually is) written "||", as in C.
The "orelse" notation resembles that of Standard ML.
Non-shortcircuiting OR also exists, and can be written "+". *)
case+ ng of (* <-- the + sign means all cases must have a match. *)
| @(a, b, c, d, e, f, g, h) =>
too_big (a) orelse too_big (b) orelse
too_big (c) orelse too_big (d) orelse
too_big (e) orelse too_big (f) orelse
too_big (g) orelse too_big (h)
fn
infinitize (term : exrat) : Option exrat =
if abs (term) >= abs (!ng8_infinitization_threshold) then
None ()
else
Some term
val no_terms_source : term_generator =
lam () =<cloref1> None ()
fn
divide (a : exrat, b : exrat) : @(exrat, exrat) =
if iseqz b then
@(exrat_make (0, 1), exrat_make (0, 1))
else
(* Do integer division of the numerators of a and b. The following
particular function does floor division if the divisor is
positive, ceiling division if the divisor is negative. Thus the
remainder is never negative. *)
exrat_numerator_euclid_division (a, b)
implement
ng8_apply ng =
lam (x, y) =>
let
val ng : ref ng8 = ref ng
and xsource : ref term_generator =
ref (continued_fraction_to_term_generator x)
and ysource : ref term_generator =
ref (continued_fraction_to_term_generator y)
fn
all_b_are_zero () : bool =
let
val @(_, _, _, _, b12, b1, b2, b) = !ng
in
(* Instead of the Standard ML-like notation "andalso", one
may (and usually does) use the C-like notation
"&&". There is also non-shortcircuiting AND, written
"*". *)
iseqz b andalso
iseqz b2 andalso
iseqz b1 andalso
iseqz b12
end
fn
all_four_equal (a : exrat, b : exrat,
c : exrat, d : exrat) : bool =
a = b && a = c && a = d
fn
absorb_x_term () =
let
val @(a12, a1, a2, a, b12, b1, b2, b) = !ng
in
case (!xsource) () of
| Some term =>
let
val new_ng = (a2 + (a12 * term),
a + (a1 * term), a12, a1,
b2 + (b12 * term),
b + (b1 * term), b12, b1)
in
if any_too_big new_ng then
(* Pretend all further x terms are infinite. *)
(!ng := @(a12, a1, a12, a1, b12, b1, b12, b1);
!xsource := no_terms_source)
else
!ng := new_ng
end
| None () =>
!ng := @(a12, a1, a12, a1, b12, b1, b12, b1)
end
fn
absorb_y_term () =
let
val @(a12, a1, a2, a, b12, b1, b2, b) = !ng
in
case (!ysource) () of
| Some term =>
let
val new_ng = (a1 + (a12 * term), a12,
a + (a2 * term), a2,
b1 + (b12 * term), b12,
b + (b2 * term), b2)
in
if any_too_big new_ng then
(* Pretend all further y terms are infinite. *)
(!ng := @(a12, a12, a2, a2, b12, b12, b2, b2);
!ysource := no_terms_source)
else
!ng := new_ng
end
| None () =>
!ng := @(a12, a12, a2, a2, b12, b12, b2, b2)
end
fun
loop () =
(* ATS2 can do mutual recursion with proper tail calls, but,
to stay closer to the Standard ML code, here I use only
single tail recursion. To do mutual recursion with proper
tail calls, one says "fnx" instead of "fun". *)
if all_b_are_zero () then
None () (* There are no more terms to output. *)
else
let
val @(_, _, _, _, b12, b1, b2, b) = !ng
in
if iseqz b andalso iseqz b2 then
(absorb_x_term (); loop ())
else if iseqz b orelse iseqz b2 then
(absorb_y_term (); loop ())
else if iseqz b1 then
(absorb_x_term (); loop ())
else
let
val @(a12, a1, a2, a, _, _, _, _) = !ng
val @(q12, r12) = divide (a12, b12)
and @(q1, r1) = divide (a1, b1)
and @(q2, r2) = divide (a2, b2)
and @(q, r) = divide (a, b)
in
if isneqz b12 andalso
all_four_equal (q12, q1, q2, q) then
(!ng := (b12, b1, b2, b, r12, r1, r2, r);
(* Return a term--or, if a magnitude threshold is
reached, return no more terms . *)
infinitize q)
else
let
(* Put a1, a2, and a over a common denominator and
compare some magnitudes. (SIDE NOTE: We are
representing big integers as EXACT rationals
with denominator one, so in fact could have put
a1, a2, and a over their respective
denominators and compared the
fractions. However, I have retained the
phrasing of the Standard ML program.) *)
val n1 = a1 * b2 * b
and n2 = a2 * b1 * b
and n = a * b1 * b2
in
if abs (n1 - n) > abs (n2 - n) then
(absorb_x_term (); loop ())
else
(absorb_y_term (); loop ())
end
end
end
in
continued_fraction_make (lam () =<cloref1> loop ())
end
(* A macro definition: *)
macdef make_op (tuple) = ng8_apply (ng8_make_int ,(tuple))
implement ng8_apply_add = make_op @(0, 1, 1, 0, 0, 0, 0, 1)
implement ng8_apply_sub = make_op @(0, 1, ~1, 0, 0, 0, 0, 1)
implement ng8_apply_mul = make_op @(1, 0, 0, 0, 0, 0, 0, 1)
implement ng8_apply_div = make_op @(0, 1, 0, 0, 0, 0, 1, 0)
(* Here the closure is "wrapped" in an ordinary function. *)
val _ng8_apply_neg = make_op @(0, 0, ~1, 0, 0, 0, 0, 1)
implement ng8_apply_neg cf = _ng8_apply_neg (cf, cf)
val _reciprocal = make_op @(0, 0, 0, 1, 0, 1, 0, 0)
implement
ng8_apply_pow (cf, i) =
let
macdef reciprocal cf = _reciprocal (,(cf), ,(cf))
fun
loop (x : continued_fraction,
n : int,
accum : continued_fraction) : continued_fraction =
if 1 < n then
let
val nhalf = n / 2
and xsquare = x * x
in
if nhalf + nhalf <> n then
loop (xsquare, nhalf, accum * x)
else
loop (xsquare, nhalf, accum)
end
else if n = 1 then
accum * x
else
accum
in
if 0 <= i then
loop (cf, i, one)
else
reciprocal (loop (cf, ~i, one))
end
implement zero = i2cf 0
implement one = i2cf 1
implement two = i2cf 2
implement three = i2cf 3
implement four = i2cf 4
implement one_fourth = r2cf (1, 4)
implement one_third = r2cf (1, 3)
implement one_half = r2cf (1, 2)
implement two_thirds = r2cf (2, 3)
implement three_fourths = r2cf (3, 4)
implement golden_ratio = constant_term_cf 1
implement silver_ratio = constant_term_cf 2
implement sqrt2 = silver_ratio - one
implement sqrt5 = (two * golden_ratio) - one
</syntaxhighlight>
The third file is the main program. It is called <code>continued-fraction-task.dats</code>:
<syntaxhighlight lang="ats">
(* Main program. *)
(*
Install ats2-xprelude, being sure to enable GMP support:
https://sourceforge.net/p/chemoelectric/ats2-xprelude
If you have it installed already, there might have been bugfixes
since. So try updating.
Then, to compile the program:
patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW \
$(pkg-config --cflags ats2-xprelude) \
$(pkg-config --variable=PATSCCFLAGS ats2-xprelude) \
continued-fraction-task.dats continued_fraction.{s,d}ats \
$(pkg-config --libs ats2-xprelude) -lgc -lm
*)
(* Load templates from the ATS2 prelude: *)
#include "share/atspre_staload.hats"
staload "continued_fraction.sats" (* Programmer access to exported stuff. *)
dynload "continued_fraction.dats" (* Initialize the "val". *)
(* Load declarations and templates from ats2-xprelude: *)
#include "xprelude/HATS/xprelude.hats"
staload "xprelude/SATS/exrat.sats"
staload _ = "xprelude/DATS/exrat.dats"
fn
make_pad (n : size_t) : string =
let
val n = g1ofg0 n
prval () = lemma_g1uint_param n
implement string_tabulate$fopr<> _ = ' '
in
strnptr2string (string_tabulate<> n)
end
fn
show_with_note (expression : string,
cf : continued_fraction,
note : string) : void =
let
val cf_str = cf2string cf
val expr_sz = strlen expression
and cf_sz = strlen cf_str
and note_sz = strlen note
val expr_pad_sz = max (i2sz 19 - expr_sz, i2sz 0)
and cf_pad_sz =
if iseqz note_sz then
i2sz 0
else
max (i2sz 48 - cf_sz, i2sz 0)
val expr_pad = make_pad expr_pad_sz
and cf_pad = make_pad cf_pad_sz
in
println! (expr_pad, expression, " => ",
cf_str, cf_pad, note)
end
fn
show_without_note (expression : string,
cf : continued_fraction) : void =
show_with_note (expression, cf, "")
overload show with show_with_note
overload show with show_without_note
implement
main0 () = (* A main that takes no arguments and returns 0. *)
begin
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", silver_ratio, "(1 + sqrt(2))");
show ("sqrt2", sqrt2);
show ("sqrt5", sqrt5);
show ("1/4", one_fourth);
show ("1/3", one_third);
show ("1/2", one_half);
show ("2/3", two_thirds);
show ("3/4", three_fourths);
show ("13/11", r2cf (13, 11));
show ("22/7", r2cf (22, 7), "approximately pi");
show ("0", zero);
show ("1", one);
show ("2", two);
show ("3", three);
show ("4", four);
show ("4 + 3", four + three);
show ("4 - 3", four - three);
show ("4 * 3", four * three);
show ("4 / 3", four / three);
show ("4 ** 3", four ** 3);
show ("4 ** (-3)", four ** (~3));
show ("negative 4", ~four);
show ("(1 + 1/sqrt(2))/2",
(one + (one / sqrt2)) / two, "method 1");
show ("(1 + 1/sqrt(2))/2",
silver_ratio * (sqrt2 ** (~3)), "method 2");
show ("(1 + 1/sqrt(2))/2",
((silver_ratio ** 2) + one) / (four * two), "method 3");
show ("sqrt2 + sqrt2", sqrt2 + sqrt2);
show ("sqrt2 - sqrt2", sqrt2 - sqrt2);
show ("sqrt2 * sqrt2", sqrt2 * sqrt2);
show ("sqrt2 / sqrt2", sqrt2 / sqrt2);
end
</syntaxhighlight>
{{out}}
To compile the program, you might try something like the following (assuming you have Boehm GC):
<pre>patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW $(pkg-config --cflags ats2-xprelude) $(pkg-config --variable=PATSCCFLAGS ats2-xprelude) continued-fraction-task.dats continued_fraction.{s,d}ats $(pkg-config --libs ats2-xprelude) -lgc -lm</pre>
You have to specify some C language standard, because patscc defaults to C99.
Then run the program by typing <pre>./a.out</pre>
The output should resemble that of the Standard ML program from which the ATS was translated. Minus signs might look different:
<pre>
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))
sqrt2 => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt5 => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...]
1/4 => [0;4]
1/3 => [0;3]
1/2 => [0;2]
2/3 => [0;1,2]
3/4 => [0;1,3]
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
0 => [0]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
4 + 3 => [7]
4 - 3 => [1]
4 * 3 => [12]
4 / 3 => [1;3]
4 ** 3 => [64]
4 ** (-3) => [0;64]
negative 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
sqrt2 + sqrt2 => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt2 - sqrt2 => [0]
sqrt2 * sqrt2 => [2]
sqrt2 / sqrt2 => [1]
</pre>
Line 1,371 ⟶ 2,965:
size_t ix;
size_t iy;
bool
bool yoverflow;
};
Line 1,428 ⟶ 3,023:
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->
env->ix += 1;
mpz_set (env->ng[ng8a2], tmp[ng8a12]);
Line 1,442 ⟶ 3,037:
if (ng8_too_big (env->ng))
{
env->
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a1], tmp[ng8a1]);
Line 1,457 ⟶ 3,052:
for (int i = 0; i != 8; i += 1)
mpz_init_set (tmp[i], env->ng[i]);
mpz_t *term = (!env->
env->iy += 1;
mpz_set (env->ng[ng8a1], tmp[ng8a12]);
Line 1,471 ⟶ 3,066:
if (ng8_too_big (env->ng))
{
env->
mpz_set (env->ng[ng8a12], tmp[ng8a12]);
mpz_set (env->ng[ng8a2], tmp[ng8a2]);
Line 1,577 ⟶ 3,172:
env->ix = 0;
env->iy = 0;
env->
env->yoverflow = false;
return make_cf (return_ng8_term, env);
}
Line 2,393 ⟶ 3,989:
private Index ix;
private Index iy;
private bool
private bool yoverflow;
//
Line 2,410 ⟶ 4,007:
ix = 0;
iy = 0;
yoverflow = false;
}
Line 2,479 ⟶ 4,077:
{
Nullable!Term term;
if (!
term = x[ix];
ix += 1;
Line 2,497 ⟶ 4,095:
ng = new NG8 (ng.a12, ng.a1, ng.a12, ng.a1,
ng.b12, ng.b1, ng.b12, ng.b1);
}
}
Line 2,508 ⟶ 4,106:
{
Nullable!Term term;
if (!
term = y[iy];
iy += 1;
Line 2,524 ⟶ 4,122:
ng = new NG8 (ng.a12, ng.a12, ng.a2, ng.a2,
ng.b12, ng.b12, ng.b2, ng.b2);
}
}
Line 2,622 ⟶ 4,220:
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
</pre>
=={{header|Fortran}}==
{{trans|Python}}
{{trans|ObjectIcon}}
This program includes a primitive module for multiple-precision integer arithmetic. It is adequate for the task.
Parts of the program might assume two's-complement representation of signed integers. The requirement that integers be two's-complement seems to me unlikely ever to become a part of Fortran standards (even though it will be required in future C standards).
<syntaxhighlight lang="fortran">
!---------------------------------------------------------------------
module big_integers ! Big (but not very big) integers.
! NOTE: I assume that iachar and achar do not alter the most
! significant bit.
use, intrinsic :: iso_fortran_env, only: int16
implicit none
private
public :: big_integer
public :: integer2big
public :: string2big
public :: big2string
public :: big_sgn
public :: big_cmp, big_cmpabs
public :: big_neg, big_abs
public :: big_addabs, big_add
public :: big_subabs, big_sub
public :: big_mul ! One might also include a big_muladd.
public :: big_divrem ! One could also include big_div and big_rem.
public :: operator(+)
public :: operator(-)
public :: operator(*)
type :: big_integer
! The representation is sign-magnitude. The radix is 256, which
! is not speed-efficient, but which seemed relatively easy to
! work with if one were writing in standard Fortran (and assuming
! iachar and achar were "8-bit clean").
logical :: sign = .false. ! .false. => +sign, .true. => -sign.
character, allocatable :: bytes(:)
end type big_integer
character, parameter :: zero = achar (0)
character, parameter :: one = achar (1)
! An integer type capable of holding an unsigned 8-bit value.
integer, parameter :: bytekind = int16
interface operator(+)
module procedure big_add
end interface
interface operator(-)
module procedure big_neg
module procedure big_sub
end interface
interface operator(*)
module procedure big_mul
end interface
contains
elemental function logical2byte (bool) result (byte)
logical, intent(in) :: bool
character :: byte
if (bool) then
byte = one
else
byte = zero
end if
end function logical2byte
elemental function logical2i (bool) result (i)
logical, intent(in) :: bool
integer :: i
if (bool) then
i = 1
else
i = 0
end if
end function logical2i
elemental function byte2i (c) result (i)
character, intent(in) :: c
integer :: i
i = iachar (c)
end function byte2i
elemental function i2byte (i) result (c)
integer, intent(in) :: i
character :: c
c = achar (i)
end function i2byte
elemental function byte2bk (c) result (i)
character, intent(in) :: c
integer(bytekind) :: i
i = iachar (c, kind = bytekind)
end function byte2bk
elemental function bk2byte (i) result (c)
integer(bytekind), intent(in) :: i
character :: c
c = achar (i)
end function bk2byte
elemental function bk2i (i) result (j)
integer(bytekind), intent(in) :: i
integer :: j
j = int (i)
end function bk2i
elemental function i2bk (i) result (j)
integer, intent(in) :: i
integer(bytekind) :: j
j = int (iand (i, 255), kind = bytekind)
end function i2bk
! Left shift of the least significant 8 bits of a bytekind integer.
elemental function lshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ishft (ibits (a, 0, 8 - i), i)
end function lshftbk
! Right shift of the least significant 8 bits of a bytekind integer.
elemental function rshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ibits (a, i, 8 - i)
end function rshftbk
! Left shift an integer.
elemental function lshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, i)
end function lshfti
! Right shift an integer.
elemental function rshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, -i)
end function rshfti
function integer2big (i) result (a)
integer, intent(in) :: i
type(big_integer), allocatable :: a
!
! To write a more efficient implementation of this procedure is
! left as an exercise for the reader.
!
character(len = 100) :: buffer
write (buffer, '(I0)') i
a = string2big (trim (buffer))
end function integer2big
function string2big (s) result (a)
character(len = *), intent(in) :: s
type(big_integer), allocatable :: a
integer :: n, i, istart, iend
integer :: digit
if ((s(1:1) == '-') .or. s(1:1) == '+') then
istart = 2
else
istart = 1
end if
iend = len (s)
n = (iend - istart + 2) / 2
allocate (a)
allocate (a%bytes(n))
a%bytes = zero
do i = istart, iend
digit = ichar (s(i:i)) - ichar ('0')
if (digit < 0 .or. 9 < digit) error stop
a = short_multiplication (a, 10)
a = short_addition (a, digit)
end do
a%sign = (s(1:1) == '-')
call normalize (a)
end function string2big
function big2string (a) result (s)
type(big_integer), intent(in) :: a
character(len = :), allocatable :: s
type(big_integer), allocatable :: q
integer :: r
integer :: sgn
sgn = big_sgn (a)
if (sgn == 0) then
s = '0'
else
q = a
s = ''
do while (big_sgn (q) /= 0)
call short_division (q, 10, q, r)
s = achar (r + ichar ('0')) // s
end do
if (sgn < 0) s = '-' // s
end if
end function big2string
function big_sgn (a) result (sgn)
type(big_integer), intent(in) :: a
integer :: sgn
integer :: n, i
n = size (a%bytes)
i = 1
sgn = 1234
do while (sgn == 1234)
if (i == n + 1) then
sgn = 0
else if (a%bytes(i) /= zero) then
if (a%sign) then
sgn = -1
else
sgn = 1
end if
else
i = i + 1
end if
end do
end function big_sgn
function big_cmp (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp
if (a%sign) then
if (b%sign) then
cmp = -big_cmpabs (a, b)
else
cmp = -1
end if
else
if (b%sign) then
cmp = 1
else
cmp = big_cmpabs (a, b)
end if
end if
end function big_cmp
function big_cmpabs (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp
integer :: n, i
integer :: ia, ib
cmp = 1234
n = max (size (a%bytes), size (b%bytes))
i = n
do while (cmp == 1234)
if (i == 0) then
cmp = 0
else
ia = byteval (a, i)
ib = byteval (b, i)
if (ia < ib) then
cmp = -1
else if (ia > ib) then
cmp = 1
else
i = i - 1
end if
end if
end do
end function big_cmpabs
function big_neg (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .not. c%sign
end function big_neg
function big_abs (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .false.
end function big_abs
function big_add (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c
logical :: sign
if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_addabs (a, b)
sign = .true.
else ! a <= 0, b >= 0
c = big_subabs (a, b)
sign = .not. c%sign
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_subabs (a, b)
sign = c%sign
else ! a >= 0, b >= 0
c = big_addabs (a, b)
sign = .false.
end if
end if
c%sign = sign
end function big_add
function big_sub (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c
logical :: sign
if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_subabs (a, b)
sign = .not. c%sign
else ! a <= 0, b >= 0
c = big_addabs (a, b)
sign = .true.
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_addabs (a, b)
sign = .false.
else ! a >= 0, b >= 0
c = big_subabs (a, b)
sign = c%sign
end if
end if
c%sign = sign
end function big_sub
function big_addabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
! Compute abs(a) + abs(b).
integer :: n, nc, i
logical :: carry
type(big_integer), allocatable :: tmp
n = max (size (a%bytes), size (b%bytes))
nc = n + 1
allocate(tmp)
allocate(tmp%bytes(nc))
call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function big_addabs
function big_subabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
! Compute abs(a) - abs(b). The result is signed.
integer :: n, i
logical :: carry
type(big_integer), allocatable :: tmp
n = max (size (a%bytes), size (b%bytes))
allocate(tmp)
allocate(tmp%bytes(n))
if (big_cmpabs (a, b) >= 0) then
tmp%sign = .false.
call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
else
tmp%sign = .true.
call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry)
end do
end if
call normalize (tmp)
c = tmp
end function big_subabs
function big_mul (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
!
! This is Knuth, Volume 2, Algorithm 4.3.1M.
!
integer :: na, nb, nc
integer :: i, j
integer :: ia, ib, ic
integer :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nb = size (b%bytes)
nc = na + nb + 1
allocate (tmp)
allocate (tmp%bytes(nc))
tmp%bytes = zero
j = 1
do j = 1, nb
ib = byte2i (b%bytes(j))
if (ib /= 0) then
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = byte2i (tmp%bytes(i + j - 1))
ic = (ia * ib) + ic + carry
tmp%bytes(i + j - 1) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(na + j) = i2byte (carry)
end if
end do
tmp%sign = (a%sign .neqv. b%sign)
call normalize (tmp)
c = tmp
end function big_mul
subroutine big_divrem (a, b, q, r)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable, intent(inout) :: q, r
!
! Division with a remainder that is never negative. Equivalently,
! this is floor division if the divisor is positive, and ceiling
! division if the divisor is negative.
!
! See Raymond T. Boute, "The Euclidean definition of the functions
! div and mod", ACM Transactions on Programming Languages and
! Systems, Volume 14, Issue 2, pp. 127-144.
! https://doi.org/10.1145/128861.128862
!
call nonnegative_division (a, b, .true., .true., q, r)
if (a%sign) then
if (big_sgn (r) /= 0) then
q = short_addition (q, 1)
r = big_sub (big_abs (b), r)
end if
q%sign = .not. b%sign
else
q%sign = b%sign
end if
end subroutine big_divrem
function short_addition (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c
! Compute abs(a) + b.
integer :: na, nc, i
logical :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nc = na + 1
allocate(tmp)
allocate(tmp%bytes(nc))
call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry)
do i = 2, na
call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function short_addition
function short_multiplication (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c
integer :: i, na, nc
integer :: ia, ic
integer :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nc = na + 1
allocate (tmp)
allocate (tmp%bytes(nc))
tmp%sign = a%sign
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = (ia * b) + carry
tmp%bytes(i) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(nc) = i2byte (carry)
call normalize (tmp)
c = tmp
end function short_multiplication
! Division without regard to signs.
subroutine nonnegative_division (a, b, want_q, want_r, q, r)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_q, want_r
type(big_integer), intent(inout), allocatable :: q, r
integer :: na, nb
integer :: remainder
na = size (a%bytes)
nb = size (b%bytes)
! It is an error if b has "significant" zero-bytes or is equal to
! zero.
if (b%bytes(nb) == zero) error stop
if (nb == 1) then
if (want_q) then
call short_division (a, byte2i (b%bytes(1)), q, remainder)
else
block
type(big_integer), allocatable :: bit_bucket
call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder)
end block
end if
if (want_r) then
if (allocated (r)) deallocate (r)
allocate (r)
allocate (r%bytes(1))
r%bytes(1) = i2byte (remainder)
end if
else
if (na >= nb) then
call long_division (a, b, want_q, want_r, q, r)
else
if (want_q) q = string2big ("0")
if (want_r) r = a
end if
end if
end subroutine nonnegative_division
subroutine short_division (a, b, q, r)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), intent(inout), allocatable :: q
integer, intent(inout) :: r
!
! This is Knuth, Volume 2, Exercise 4.3.1.16.
!
! The divisor is assumed to be positive.
!
integer :: n, i
integer :: ia, ib, iq
type(big_integer), allocatable :: tmp
ib = b
n = size (a%bytes)
allocate (tmp)
allocate (tmp%bytes(n))
r = 0
do i = n, 1, -1
ia = (256 * r) + byte2i (a%bytes(i))
iq = ia / ib
r = mod (ia, ib)
tmp%bytes(i) = i2byte (iq)
end do
tmp%sign = a%sign
call normalize (tmp)
q = tmp
end subroutine short_division
subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_quotient, want_remainder
type(big_integer), intent(inout), allocatable :: quotient
type(big_integer), intent(inout), allocatable :: remainder
!
! This is Knuth, Volume 2, Algorithm 4.3.1D.
!
! We do not deal here with the signs of the inputs and outputs.
!
! It is assumed size(a%bytes) >= size(b%bytes), and that b has no
! leading zero-bytes and is at least two bytes long. If b is one
! byte long and nonzero, use short division.
!
integer :: na, nb, m, n
integer :: num_lz, num_nonlz
integer :: j
integer :: qhat
logical :: carry
!
! We will NOT be working with VERY large numbers, and so it will
! be safe to put temporary storage on the stack. (Note: your
! Fortran might put this storage in a heap instead of the stack.)
!
! v = b, normalized to put its most significant 1-bit all the
! way left.
!
! u = a, shifted left by the same amount as b.
!
! q = the quotient.
!
! The remainder, although shifted left, will end up in u.
!
integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes))
integer(bytekind) :: v(0:size (b%bytes) - 1)
integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes))
na = size (a%bytes)
nb = size (b%bytes)
n = nb
m = na - nb
! In the most significant byte of the divisor, find the number of
! leading zero bits, and the number of bits after that.
block
integer(bytekind) :: tmp
tmp = byte2bk (b%bytes(n))
num_nonlz = bit_size (tmp) - leadz (tmp)
num_lz = 8 - num_nonlz
end block
call normalize_v (b%bytes) ! Make the most significant bit of v be one.
call normalize_u (a%bytes) ! Shifted by the same amount as v.
! Assure ourselves that the most significant bit of v is a one.
if (.not. btest (v(n - 1), 7)) error stop
do j = m, 0, -1
call calculate_qhat (qhat)
call multiply_and_subtract (carry)
q(j) = i2bk (qhat)
if (carry) call add_back
end do
if (want_quotient) then
if (allocated (quotient)) deallocate (quotient)
allocate (quotient)
allocate (quotient%bytes(m + 1))
quotient%bytes = bk2byte (q)
call normalize (quotient)
end if
if (want_remainder) then
if (allocated (remainder)) deallocate (remainder)
allocate (remainder)
allocate (remainder%bytes(n))
call unnormalize_u (remainder%bytes)
call normalize (remainder)
end if
contains
subroutine normalize_v (b_bytes)
character, intent(in) :: b_bytes(n)
!
! Normalize v so its most significant bit is a one. Any
! normalization factor that achieves this goal will suffice; we
! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: btmp(0:n - 1)
btmp = byte2bk (b_bytes)
v(0) = lshftbk (btmp(0), num_lz)
do i = 1, n - 1
v(i) = ior (lshftbk (btmp(i), num_lz), &
& rshftbk (btmp(i - 1), num_nonlz))
end do
end subroutine normalize_v
subroutine normalize_u (a_bytes)
character, intent(in) :: a_bytes(m + n)
!
! Shift a leftwards to get u. Shift by as much as b was shifted
! to get v.
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: atmp(0:m + n - 1)
atmp = byte2bk (a_bytes)
u(0) = lshftbk (atmp(0), num_lz)
do i = 1, m + n - 1
u(i) = ior (lshftbk (atmp(i), num_lz), &
& rshftbk (atmp(i - 1), num_nonlz))
end do
u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz)
end subroutine normalize_u
subroutine unnormalize_u (r_bytes)
character, intent(out) :: r_bytes(n)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: rtmp(0:n - 1)
do i = 0, n - 1
rtmp(i) = ior (rshftbk (u(i), num_lz), &
& lshftbk (u(i + 1), num_nonlz))
end do
rtmp(n - 1) = rshftbk (u(n - 1), num_lz)
r_bytes = bk2byte (rtmp)
end subroutine unnormalize_u
subroutine calculate_qhat (qhat)
integer, intent(out) :: qhat
integer :: itmp, rhat
logical :: adjust
itmp = ior (lshfti (bk2i (u(j + n)), 8), &
& bk2i (u(j + n - 1)))
qhat = itmp / bk2i (v(n - 1))
rhat = mod (itmp, bk2i (v(n - 1)))
adjust = .true.
do while (adjust)
if (rshfti (qhat, 8) /= 0) then
continue
else if (qhat * bk2i (v(n - 2)) &
& > ior (lshfti (rhat, 8), &
& bk2i (u(j + n - 2)))) then
continue
else
adjust = .false.
end if
if (adjust) then
qhat = qhat - 1
rhat = rhat + bk2i (v(n - 1))
if (rshfti (rhat, 8) == 0) then
adjust = .false.
end if
end if
end do
end subroutine calculate_qhat
subroutine multiply_and_subtract (carry)
logical, intent(out) :: carry
integer :: i
integer :: qhat_v
integer :: mul_carry, sub_carry
integer :: diff
mul_carry = 0
sub_carry = 0
do i = 0, n
! Multiplication.
qhat_v = mul_carry
if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i)))
mul_carry = rshfti (qhat_v, 8)
qhat_v = iand (qhat_v, 255)
! Subtraction.
diff = bk2i (u(j + i)) - qhat_v + sub_carry
sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1.
u(j + i) = i2bk (diff)
end do
carry = (sub_carry /= 0)
end subroutine multiply_and_subtract
subroutine add_back
integer :: i, carry, sum
q(j) = q(j) - 1_bytekind
carry = 0
do i = 0, n - 1
sum = bk2i (u(j + i)) + bk2i (v(i)) + carry
carry = ishft (sum, -8)
u(j + i) = i2bk (sum)
end do
end subroutine add_back
end subroutine long_division
subroutine add_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out
integer :: ia, ib, ic
ia = byte2i (a)
if (carry_in) ia = ia + 1
ib = byte2i (b)
ic = ia + ib
c = i2byte (iand (ic, 255))
carry_out = (ic >= 256)
end subroutine add_bytes
subroutine sub_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out
integer :: ia, ib, ic
ia = byte2i (a)
ib = byte2i (b)
if (carry_in) ib = ib + 1
ic = ia - ib
carry_out = (ic < 0)
if (carry_out) ic = ic + 256
c = i2byte (iand (ic, 255))
end subroutine sub_bytes
function get_byte (a, i) result (byte)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
character :: byte
if (size (a%bytes) < i) then
byte = zero
else
byte = a%bytes(i)
end if
end function get_byte
function byteval (a, i) result (v)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
integer :: v
if (size (a%bytes) < i) then
v = 0
else
v = byte2i (a%bytes(i))
end if
end function byteval
subroutine normalize (a)
type(big_integer), intent(inout) :: a
logical :: done
integer :: i
character, allocatable :: fewer_bytes(:)
! Shorten to the minimum number of bytes.
i = size (a%bytes)
done = .false.
do while (.not. done)
if (i == 1) then
done = .true.
else if (a%bytes(i) /= zero) then
done = .true.
else
i = i - 1
end if
end do
if (i /= size (a%bytes)) then
allocate (fewer_bytes (i))
fewer_bytes = a%bytes(1:i)
call move_alloc (fewer_bytes, a%bytes)
end if
! If the magnitude is zero, then clear the sign bit.
if (size (a%bytes) == 1) then
if (a%bytes(1) == zero) then
a%sign = .false.
end if
end if
end subroutine normalize
end module big_integers
!---------------------------------------------------------------------
module continued_fractions
use, non_intrinsic :: big_integers
implicit none
private
public :: continued_fraction
public :: term_generator
public :: term_generator_procedure
public :: make_continued_fraction
public :: i2cf
public :: make_integer_continued_fraction
public :: make_integer_continued_fraction_from_integer
public :: constant_term_cf
public :: make_constant_term_continued_fraction
public :: make_constant_term_continued_fraction_from_integer
public :: apply_ng8
public :: apply_ng8_big_integers
public :: apply_ng8_integers
public :: ng8_coefficient_threshold
public :: ng8_term_threshold
public :: add_continued_fractions
public :: subtract_continued_fractions
public :: multiply_continued_fractions
public :: divide_continued_fractions
public :: cf2string
public :: continued_fraction_to_string_given_max_terms
public :: continued_fraction_to_string_with_default_max_terms
public :: default_continued_fraction_max_terms
type :: continued_fraction
class(continued_fraction_record), pointer, private :: p => null ()
contains
procedure, pass :: get_term => get_continued_fraction_term
procedure, pass :: term_exists => continued_fraction_term_exists
procedure, pass :: term => continued_fraction_term
procedure, pass :: to_string => continued_fraction_to_string_with_default_max_terms
procedure, pass :: add => add_continued_fractions
generic :: operator(+) => add
procedure, pass :: subtract => subtract_continued_fractions
generic :: operator(-) => subtract
procedure, pass :: multiply => multiply_continued_fractions
generic :: operator(*) => multiply
procedure, pass :: divide => divide_continued_fractions
generic :: operator(/) => divide
procedure, pass, private :: continued_fraction_make_new_ref
generic :: assignment(=) => continued_fraction_make_new_ref
final :: continued_fraction_final
end type continued_fraction
type :: continued_fraction_record
logical, private :: terminated = .false. ! No more terms?
integer, private :: m = 0 ! No. of terms memoized.
type(big_integer), private, allocatable :: memo(:) ! Memoized terms.
class(term_generator), pointer :: gen ! Where terms come from.
integer :: refcount = 0
contains
procedure, pass :: get_term => get_continued_fraction_record_term
procedure, pass :: term_exists => continued_fraction_record_term_exists
procedure, pass :: term => continued_fraction_record_term
final :: continued_fraction_record_final
end type continued_fraction_record
type, abstract :: term_generator
contains
procedure(term_generator_procedure), pass, deferred :: generate
end type term_generator
interface
subroutine term_generator_procedure (gen, term_exists, term)
import term_generator
import big_integer
class(term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
end subroutine term_generator_procedure
end interface
type, extends (term_generator) :: integer_term_generator
type(big_integer), allocatable :: term
logical :: no_more_terms = .false.
contains
procedure, pass :: generate => integer_term_generator_generate
end type integer_term_generator
type, extends (term_generator) :: constant_term_generator
type(big_integer), allocatable :: term
contains
procedure, pass :: generate => constant_term_generator_generate
end type constant_term_generator
type, extends (term_generator) :: ng8_term_generator
type(big_integer), allocatable :: a12, a1, a2, a
type(big_integer), allocatable :: b12, b1, b2, b
type(continued_fraction) :: x, y
integer :: ix = 0
integer :: iy = 0
logical :: x_overflow = .false.
logical :: y_overflow = .false.
contains
procedure, pass :: generate => ng8_term_generator_generate
end type ng8_term_generator
interface i2cf
module procedure make_integer_continued_fraction
module procedure make_integer_continued_fraction_from_integer
end interface i2cf
interface constant_term_cf
module procedure make_constant_term_continued_fraction
module procedure make_constant_term_continued_fraction_from_integer
end interface constant_term_cf
interface apply_ng8
module procedure apply_ng8_big_integers
module procedure apply_ng8_integers
end interface apply_ng8
interface cf2string
module procedure continued_fraction_to_string_given_max_terms
module procedure continued_fraction_to_string_with_default_max_terms
end interface cf2string
integer :: default_continued_fraction_max_terms = 20
type(big_integer), allocatable :: ng8_coefficient_threshold
type(big_integer), allocatable :: ng8_term_threshold
contains
subroutine continued_fraction_make_new_ref (dst, src)
class(continued_fraction), intent(inout) :: dst
class(continued_fraction), intent(in) :: src
if (associated (dst%p)) deallocate (dst%p)
dst%p => src%p
dst%p%refcount = dst%p%refcount + 1
end subroutine continued_fraction_make_new_ref
subroutine continued_fraction_final (cf)
type(continued_fraction), intent(inout) :: cf
cf%p%refcount = cf%p%refcount - 1
if (cf%p%refcount == 0) deallocate (cf%p)
end subroutine continued_fraction_final
function make_continued_fraction (gen) result (cf)
class(term_generator), pointer, intent(in) :: gen
type(continued_fraction) :: cf
allocate (cf%p)
allocate (cf%p%memo(0:31)) ! The starting size is arbitrary.
cf%p%gen => gen
cf%p%refcount = cf%p%refcount + 1
end function make_continued_fraction
subroutine continued_fraction_record_final (cfrec)
type(continued_fraction_record), intent(inout) :: cfrec
deallocate (cfrec%gen)
end subroutine continued_fraction_record_final
function make_integer_continued_fraction (bigint) result (cf)
type(big_integer), intent(in) :: bigint
type(continued_fraction) :: cf
class(integer_term_generator), pointer :: gen
allocate (gen)
gen%term = bigint
cf = make_continued_fraction (gen)
end function make_integer_continued_fraction
function make_integer_continued_fraction_from_integer (i) result (cf)
integer, intent(in) :: i
type(continued_fraction) :: cf
cf = make_integer_continued_fraction (integer2big (i))
end function make_integer_continued_fraction_from_integer
subroutine integer_term_generator_generate (gen, term_exists, term)
class(integer_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
term_exists = (.not. gen%no_more_terms)
if (term_exists) term = gen%term
gen%no_more_terms = .true.
end subroutine integer_term_generator_generate
function make_constant_term_continued_fraction (bigint) result (cf)
type(big_integer), intent(in) :: bigint
type(continued_fraction) :: cf
class(constant_term_generator), pointer :: gen
allocate (gen)
gen%term = bigint
cf = make_continued_fraction (gen)
end function make_constant_term_continued_fraction
function make_constant_term_continued_fraction_from_integer (i) result (cf)
integer, intent(in) :: i
type(continued_fraction) :: cf
cf = make_constant_term_continued_fraction (integer2big (i))
end function make_constant_term_continued_fraction_from_integer
subroutine constant_term_generator_generate (gen, term_exists, term)
class(constant_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
term_exists = .true.
if (term_exists) term = gen%term
end subroutine constant_term_generator_generate
function apply_ng8_big_integers (a12, a1, a2, a, &
& b12, b1, b2, b, x, y) result (cf)
type(big_integer), intent(in) :: a12, a1, a2, a
type(big_integer), intent(in) :: b12, b1, b2, b
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
class(ng8_term_generator), pointer :: gen
allocate (gen)
gen%a12 = a12; gen%a1 = a1; gen%a2 = a2; gen%a = a
gen%b12 = b12; gen%b1 = b1; gen%b2 = b2; gen%b = b
gen%x = x
gen%y = y
cf = make_continued_fraction (gen)
end function apply_ng8_big_integers
function apply_ng8_integers (a12, a1, a2, a, &
& b12, b1, b2, b, x, y) result (cf)
integer, intent(in) :: a12, a1, a2, a
integer, intent(in) :: b12, b1, b2, b
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8_big_integers (integer2big (a12), &
& integer2big (a1), &
& integer2big (a2), &
& integer2big (a), &
& integer2big (b12), &
& integer2big (b1), &
& integer2big (b2), &
& integer2big (b), x, y)
end function apply_ng8_integers
subroutine ng8_term_generator_generate (gen, term_exists, term)
class(ng8_term_generator), intent(inout) :: gen
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
logical :: done
logical :: b12z, b1z, b2z, bz
type(big_integer), allocatable :: q12, r12
type(big_integer), allocatable :: q1, r1
type(big_integer), allocatable :: q2, r2
type(big_integer), allocatable :: q, r
logical :: absorb_x, absorb_y, compare_fracs
done = .false.
do while (.not. done)
absorb_x = .false.
absorb_y = .false.
compare_fracs = .false.
b12z = (big_sgn (gen%b12) == 0)
b1z = (big_sgn (gen%b1) == 0)
b2z = (big_sgn (gen%b2) == 0)
bz = (big_sgn (gen%b) == 0)
if (b12z .and. b1z .and. b2z .and. bz) then
! There are no more terms.
term_exists = .false.
done = .true.
else if (b2z .and. bz) then
absorb_x = .true.
else if (b2z .or. bz) then
absorb_y = .true.
else if (b1z) then
absorb_x = .true.
else
call big_divrem (gen%a1, gen%b1, q1, r1)
call big_divrem (gen%a2, gen%b2, q2, r2)
call big_divrem (gen%a, gen%b, q, r)
if (.not. b12z) then
call big_divrem (gen%a12, gen%b12, q12, r12)
if (big_cmp (q, q1) /= 0) then
compare_fracs = .true.
else if (big_cmp (q, q2) /= 0) then
compare_fracs = .true.
else if (big_cmp (q, q12) /= 0) then
compare_fracs = .true.
else
call output_term
done = .true.
end if
end if
end if
if (compare_fracs) call compare_fractions (absorb_x, absorb_y)
if (absorb_x) call absorb_x_term
if (absorb_y) call absorb_y_term
end do
contains
subroutine output_term
gen%a12 = gen%b12; gen%a1 = gen%b1; gen%a2 = gen%b2; gen%a = gen%b
gen%b12 = r12; gen%b1 = r1; gen%b2 = r2; gen%b = r
term_exists = (.not. treat_as_infinite (q))
if (term_exists) term = q
end subroutine output_term
subroutine compare_fractions (absorb_x, absorb_y)
logical, intent(inout) :: absorb_x, absorb_y
! Rather than compare fractions, we will put the numerators over
! a common denominator of b1*b2*b, and then compare the new
! numerators.
type(big_integer), allocatable :: n1, n2, n
n1 = gen%a1 * gen%b2 * gen%b
n2 = gen%a2 * gen%b1 * gen%b
n = gen%a * gen%b1 * gen%b2
if (big_cmpabs (n1 - n, n2 - n) > 0) then
absorb_x = .true.
else
absorb_y = .true.
end if
end subroutine compare_fractions
subroutine absorb_x_term
logical :: term_exists
type(big_integer), allocatable :: term
type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a
type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b
if (gen%x_overflow) then
term_exists = .false.
else
term_exists = gen%x%term_exists(gen%ix)
end if
new_a2 = gen%a12
new_a = gen%a1
new_b2 = gen%b12
new_b = gen%b1
if (term_exists) then
term = gen%x%term(gen%ix)
new_a12 = gen%a2 + (gen%a12 * term)
new_a1 = gen%a + (gen%a1 * term)
new_b12 = gen%b2 + (gen%b12 * term)
new_b1 = gen%b + (gen%b1 * term)
if (any_too_big (new_a12, new_a1, new_a2, new_a, &
& new_b12, new_b1, new_b2, new_b)) then
gen%x_overflow = .true.
new_a12 = gen%a12
new_a1 = gen%a1
new_b12 = gen%b12
new_b1 = gen%b1
end if
else
new_a12 = gen%a12
new_a1 = gen%a1
new_b12 = gen%b12
new_b1 = gen%b1
end if
gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a
gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b
gen%ix = gen%ix + 1
end subroutine absorb_x_term
subroutine absorb_y_term
logical :: term_exists
type(big_integer), allocatable :: term
type(big_integer), allocatable :: new_a12, new_a1, new_a2, new_a
type(big_integer), allocatable :: new_b12, new_b1, new_b2, new_b
if (gen%y_overflow) then
term_exists = .false.
else
term_exists = gen%y%term_exists(gen%iy)
end if
new_a1 = gen%a12
new_a = gen%a2
new_b1 = gen%b12
new_b = gen%b2
if (term_exists) then
term = gen%y%term(gen%iy)
new_a12 = gen%a1 + (gen%a12 * term)
new_a2 = gen%a + (gen%a2 * term)
new_b12 = gen%b1 + (gen%b12 * term)
new_b2 = gen%b + (gen%b2 * term)
if (any_too_big (new_a12, new_a1, new_a2, new_a, &
& new_b12, new_b1, new_b2, new_b)) then
gen%y_overflow = .true.
new_a12 = gen%a12
new_a2 = gen%a2
new_b12 = gen%b12
new_b2 = gen%b2
end if
else
new_a12 = gen%a12
new_a2 = gen%a2
new_b12 = gen%b12
new_b2 = gen%b2
end if
gen%a12 = new_a12; gen%a1 = new_a1; gen%a2 = new_a2; gen%a = new_a
gen%b12 = new_b12; gen%b1 = new_b1; gen%b2 = new_b2; gen%b = new_b
gen%iy = gen%iy + 1
end subroutine absorb_y_term
function any_too_big (a, b, c, d, e, f, g, h) result (yes)
type(big_integer), intent(in) :: a, b, c, d, e, f, g, h
logical :: yes
if (too_big (a)) then
yes = .true.
else if (too_big (b)) then
yes = .true.
else if (too_big (c)) then
yes = .true.
else if (too_big (d)) then
yes = .true.
else if (too_big (e)) then
yes = .true.
else if (too_big (f)) then
yes = .true.
else if (too_big (g)) then
yes = .true.
else if (too_big (h)) then
yes = .true.
else
yes = .false.
end if
end function any_too_big
function too_big (coef) result (yes)
type(big_integer), intent(in) :: coef
logical :: yes
if (.not. allocated (ng8_coefficient_threshold)) then
! 2**512
ng8_coefficient_threshold = string2big ('1340780792994259709957&
&402499820584612747936582059239337772356144372176403007354&
&697680187429816690342769003185818648605085375388281194656&
&9946433649006084096')
end if
yes = (big_cmpabs (coef, ng8_coefficient_threshold) >= 0)
end function too_big
function treat_as_infinite (term) result (yes)
type(big_integer), intent(in) :: term
logical :: yes
if (.not. allocated (ng8_term_threshold)) then
! 2**64
ng8_term_threshold = string2big ('18446744073709551616')
end if
yes = (big_cmpabs (term, ng8_term_threshold) >= 0)
end function treat_as_infinite
end subroutine ng8_term_generator_generate
function add_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1, x, y)
end function add_continued_fractions
function subtract_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1, x, y)
end function subtract_continued_fractions
function multiply_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1, x, y)
end function multiply_continued_fractions
function divide_continued_fractions (x, y) result (cf)
class(continued_fraction), intent(in) :: x, y
type(continued_fraction) :: cf
cf = apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0, x, y)
end function divide_continued_fractions
subroutine get_continued_fraction_term (cf, i, term_exists, term)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
call get_continued_fraction_record_term (cf%p, i, term_exists, term)
end subroutine get_continued_fraction_term
subroutine get_continued_fraction_record_term (cfrec, i, term_exists, term)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
logical, intent(out) :: term_exists
type(big_integer), allocatable, intent(out) :: term
if (i < 0) error stop
call update (i + 1)
term_exists = (i < cfrec%m)
if (term_exists) term = cfrec%memo(i)
contains
subroutine update (needed)
integer :: needed
logical :: term_exists1
type(big_integer), allocatable :: term1
if (.not. cfrec%terminated .and. cfrec%m < needed) then
if (size (cfrec%memo) < needed) then
block ! Allocate more storage.
type(big_integer), allocatable :: memo1(:)
allocate (memo1(0 : (2 * needed) - 1))
memo1(0:cfrec%m - 1) = cfrec%memo(0:cfrec%m - 1)
call move_alloc (memo1, cfrec%memo)
end block
end if
do while (.not. cfrec%terminated .and. cfrec%m < needed)
call cfrec%gen%generate (term_exists1, term1)
if (term_exists1) then
cfrec%memo(cfrec%m) = term1
cfrec%m = cfrec%m + 1
else
cfrec%terminated = .true.
end if
end do
end if
end subroutine update
end subroutine get_continued_fraction_record_term
function continued_fraction_term_exists (cf, i) result (term_exists)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
logical :: term_exists
term_exists = continued_fraction_record_term_exists (cf%p, i)
end function continued_fraction_term_exists
function continued_fraction_record_term_exists (cfrec, i) result (term_exists)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
logical :: term_exists
type(big_integer), allocatable :: term
call get_continued_fraction_record_term (cfrec, i, term_exists, term)
end function continued_fraction_record_term_exists
function continued_fraction_term (cf, i) result (term)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: i
type(big_integer), allocatable :: term
term = continued_fraction_record_term (cf%p, i)
end function continued_fraction_term
function continued_fraction_record_term (cfrec, i) result (term)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: i
type(big_integer), allocatable :: term
logical :: term_exists
call get_continued_fraction_record_term (cfrec, i, term_exists, term)
if (.not. term_exists) error stop
end function continued_fraction_record_term
function continued_fraction_to_string_given_max_terms (cf, max_terms) result (s)
class(continued_fraction), intent(in) :: cf
integer, intent(in) :: max_terms
character(len = :), allocatable :: s
s = continued_fraction_record_to_string_given_max_terms (cf%p, max_terms)
end function continued_fraction_to_string_given_max_terms
function continued_fraction_record_to_string_given_max_terms (cfrec, max_terms) result (s)
class(continued_fraction_record), intent(inout) :: cfrec
integer, intent(in) :: max_terms
character(len = :), allocatable :: s
integer :: i
logical :: done
i = 0
s = '['
done = .false.
do while (.not. done)
if (.not. cfrec%term_exists(i)) then
s = s // "]"
done = .true.
else if (i == max_terms) then
s = s // ",...]"
done = .true.
else
select case (i)
case (0); continue
case (1); s = s // ";"
case default; s = s // ","
end select
s = s // big2string (cfrec%term(i))
i = i + 1
end if
end do
end function continued_fraction_record_to_string_given_max_terms
function continued_fraction_to_string_with_default_max_terms (cf) result (s)
class(continued_fraction), intent(in) :: cf
character(len = :), allocatable :: s
s = continued_fraction_record_to_string_with_default_max_terms (cf%p)
end function continued_fraction_to_string_with_default_max_terms
function continued_fraction_record_to_string_with_default_max_terms (cfrec) result (s)
class(continued_fraction_record), intent(inout) :: cfrec
character(len = :), allocatable :: s
integer :: max_terms
max_terms = max (default_continued_fraction_max_terms, 1)
s = continued_fraction_record_to_string_given_max_terms (cfrec, max_terms)
end function continued_fraction_record_to_string_with_default_max_terms
end module continued_fractions
!---------------------------------------------------------------------
program bivariate_continued_fraction_task
use, non_intrinsic :: big_integers
use, non_intrinsic :: continued_fractions
implicit none
type(continued_fraction) :: golden_ratio
type(continued_fraction) :: silver_ratio
type(continued_fraction) :: sqrt2
type(continued_fraction) :: one
type(continued_fraction) :: two
type(continued_fraction) :: three
type(continued_fraction) :: four
type(continued_fraction) :: method1
type(continued_fraction) :: method2
type(continued_fraction) :: method3
block
!
! Let us start with a test of the long division routine, with some
! numbers known to trigger a bug in the first version I
! posted. That version had a buggy "add_back" routine.
!
! (How I found such numbers is easy: I used random search.)
!
type(big_integer), allocatable :: a, b, q, r
a = string2big ("95292350834616415707142739736356097545484658215015733475&
&690528634954280668802285176284181482202905629004984123564335019024318905")
b = string2big ("63677949970178275389170357551071104191609722674550547056511830750")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("5286200801181288750950358142425694618335361315503743069535407838&
&1104411448878793976933118436177295215313131557463887718741957154")
b = string2big ("54401097470188014066128968444633185848791550678521")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("74352827755975214935544861176217106447734695144315262422&
&288346418457330023596489437154599028318030933202302606937951415862696330")
b = string2big ("291979433784649910583546698460221489986784915256036707914578&
&892106828527219639012712")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("7515839498480018152556264500298705705667515770181724145893&
&9544448656273749453586884931339958104923411455488844854494605760712247")
b = string2big ("8600698996698965932302079501896131441135807557744554902970070&
&402964318496325075877027770784963001")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
a = string2big ("13370595927769020368832742717678609885835798503146654175875&
&149837801359758893206045930442389897206420586502996797614097489470778")
b = string2big ("871343613388")
call big_divrem (a, b, q, r)
if (big_sgn (((b * q) + r) - a) /= 0) error stop
end block
golden_ratio = constant_term_cf (1)
silver_ratio = constant_term_cf (2)
one = i2cf (1)
two = i2cf (2)
three = i2cf (3)
four = i2cf (4)
sqrt2 = silver_ratio - one
method1 = apply_ng8 (0, 1, 0, 0, 0, 0, 2, 0, silver_ratio, sqrt2)
method2 = apply_ng8 (1, 0, 0, 1, 0, 0, 0, 8, silver_ratio, silver_ratio)
method3 = (one / two) * (one + (one / sqrt2))
call show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2")
call show ("silver ratio", silver_ratio, "(1 + sqrt(2))")
call show ("sqrt(2)", sqrt2, "silver ratio minus 1")
call show ("13/11", i2cf (13) / i2cf (11), "")
call show ("22/7", i2cf (22) / i2cf (7), "")
call show ("1", one, "")
call show ("2", two, "")
call show ("3", three, "")
call show ("4", four, "")
call show ("(1 + 1/sqrt(2))/2", method1, "method 1")
call show ("(1 + 1/sqrt(2))/2", method2, "method 2")
call show ("(1 + 1/sqrt(2))/2", method3, "method 3")
call show ("sqrt(2) + sqrt(2)", sqrt2 + sqrt2, "")
call show ("sqrt(2) - sqrt(2)", sqrt2 - sqrt2, "")
call show ("sqrt(2) * sqrt(2)", sqrt2 * sqrt2, "")
call show ("sqrt(2) / sqrt(2)", sqrt2 / sqrt2, "")
contains
subroutine show (expression, cf, note)
character(len = *), intent(in) :: expression
class(continued_fraction), intent(in) :: cf
character(len = *), intent(in) :: note
write (*, '(A19, " => ", A, T73, A)') &
& expression, cf%to_string(), note
end subroutine show
end program bivariate_continued_fraction_task
!---------------------------------------------------------------------
</syntaxhighlight>
{{out}}
<pre>$ gfortran -O2 -g -fbounds-check -Wall -Wextra bivariate_continued_fraction_task.f90 && ./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,...] silver ratio minus 1
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]
</pre>
Line 2,816 ⟶ 6,126:
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]
</pre>
=={{header|Haskell}}==
{{trans|Mercury}}
This Haskell follows the Mercury, in using infinitely long lazy lists to represent continued fractions. There are two kinds of terms: "infinite" and "finite integer".
<syntaxhighlight lang="Haskell">
----------------------------------------------------------------------
data Term = InfiniteTerm | IntegerTerm Integer
type ContinuedFraction = [Term] -- The list should be infinitely long.
type NG8 = (Integer, Integer, Integer, Integer,
Integer, Integer, Integer, Integer)
----------------------------------------------------------------------
cf2string (cf :: ContinuedFraction) =
loop 0 "[" cf
where loop i s lst =
case lst of {
(InfiniteTerm : _) -> s ++ "]" ;
(IntegerTerm value : tail) ->
(if i == 20 then
s ++ ",...]"
else
let {
sepStr =
case i of {
0 -> "";
1 -> ";";
_ -> ","
};
termStr = show value;
s1 = s ++ sepStr ++ termStr
}
in loop (i + 1) s1 tail)
}
----------------------------------------------------------------------
repeatingTerm (term :: Term) =
term : repeatingTerm term
infiniteContinuedFraction = repeatingTerm InfiniteTerm
i2cf (i :: Integer) =
-- Continued fraction representing an integer.
IntegerTerm i : infiniteContinuedFraction
r2cf (n :: Integer) (d :: Integer) =
-- Continued fraction representing a rational number.
let (q, r) = divMod n d in
(if r == 0 then
(IntegerTerm q : infiniteContinuedFraction)
else
(IntegerTerm q : r2cf d r))
----------------------------------------------------------------------
add_cf = apply_ng8 (0, 1, 1, 0, 0, 0, 0, 1)
sub_cf = apply_ng8 (0, 1, -1, 0, 0, 0, 0, 1)
mul_cf = apply_ng8 (1, 0, 0, 0, 0, 0, 0, 1)
div_cf = apply_ng8 (0, 1, 0, 0, 0, 0, 1, 0)
apply_ng8
(ng :: NG8)
(x :: ContinuedFraction)
(y :: ContinuedFraction) =
--
let (a12, a1, a2, a, b12, b1, b2, b) = ng in
if iseqz [b12, b1, b2, b] then
infiniteContinuedFraction -- No more finite terms to output.
else if iseqz [b2, b] then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else if atLeastOne_iseqz [b2, b] then
let (ng1, x1, y1) = absorb_y_term ng x y in
apply_ng8 ng1 x1 y1
else if iseqz [b1] then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else
let {
(q12, r12) = maybeDivide a12 b12;
(q1, r1) = maybeDivide a1 b1;
(q2, r2) = maybeDivide a2 b2;
(q, r) = maybeDivide a b
}
in
if not (iseqz [b12]) && q == q12 && q == q1 && q == q2 then
-- Output a term.
(if integerExceedsInfinitizingThreshold q then
infiniteContinuedFraction
else
let new_ng = (b12, b1, b2, b, r12, r1, r2, r) in
(IntegerTerm q : apply_ng8 new_ng x y))
else
-- Put a1, a2, and a over a common denominator and compare
-- some magnitudes.
let {
n1 = a1 * b2 * b;
n2 = a2 * b1 * b;
n = a * b1 * b2
}
in
(if abs (n1 - n) > abs (n2 - n) then
let (ng1, x1, y1) = absorb_x_term ng x y in
apply_ng8 ng1 x1 y1
else
let (ng1, x1, y1) = absorb_y_term ng x y in
apply_ng8 ng1 x1 y1)
absorb_x_term
((a12, a1, a2, a, b12, b1, b2, b) :: NG8)
(x :: ContinuedFraction)
(y :: ContinuedFraction) =
--
case x of {
(IntegerTerm n : xtail) -> (
let new_ng = (a2 + (a12 * n), a + (a1 * n), a12, a1,
b2 + (b12 * n), b + (b1 * n), b12, b1) in
if (ng8ExceedsProcessingThreshold new_ng) then
-- Pretend we have reached an infinite term.
((a12, a1, a12, a1, b12, b1, b12, b1),
infiniteContinuedFraction, y)
else
(new_ng, xtail, y)
);
(InfiniteTerm : _) ->
((a12, a1, a12, a1, b12, b1, b12, b1),
infiniteContinuedFraction, y)
}
absorb_y_term
((a12, a1, a2, a, b12, b1, b2, b) :: NG8)
(x :: ContinuedFraction)
(y :: ContinuedFraction) =
--
case y of {
(IntegerTerm n : ytail) -> (
let new_ng = (a1 + (a12 * n), a12, a + (a2 * n), a2,
b1 + (b12 * n), b12, b + (b2 * n), b2) in
if (ng8ExceedsProcessingThreshold new_ng) then
-- Pretend we have reached an infinite term.
((a12, a12, a2, a2, b12, b12, b2, b2),
x, infiniteContinuedFraction)
else
(new_ng, x, ytail)
);
(InfiniteTerm : _) ->
((a12, a12, a2, a2, b12, b12, b2, b2),
x, infiniteContinuedFraction)
}
ng8ExceedsProcessingThreshold (a12, a1, a2, a,
b12, b1, b2, b) =
(integerExceedsProcessingThreshold a12 ||
integerExceedsProcessingThreshold a1 ||
integerExceedsProcessingThreshold a2 ||
integerExceedsProcessingThreshold a ||
integerExceedsProcessingThreshold b12 ||
integerExceedsProcessingThreshold b1 ||
integerExceedsProcessingThreshold b2 ||
integerExceedsProcessingThreshold b)
integerExceedsProcessingThreshold i =
abs i >= 2 ^ 512
integerExceedsInfinitizingThreshold i =
abs i >= 2 ^ 64
maybeDivide a b =
if b == 0
then (0, 0)
else divMod a b
iseqz [] = True
iseqz (head : tail) = head == 0 && iseqz tail
atLeastOne_iseqz [] = False
atLeastOne_iseqz (head : tail) = head == 0 || atLeastOne_iseqz tail
----------------------------------------------------------------------
zero = i2cf 0
one = i2cf 1
two = i2cf 2
three = i2cf 3
four = i2cf 4
one_fourth = r2cf 1 4
one_third = r2cf 1 3
one_half = r2cf 1 2
two_thirds = r2cf 2 3
three_fourths = r2cf 3 4
goldenRatio = repeatingTerm (IntegerTerm 1)
silverRatio = repeatingTerm (IntegerTerm 2)
sqrt2 = IntegerTerm 1 : silverRatio
sqrt5 = IntegerTerm 2 : repeatingTerm (IntegerTerm 4)
----------------------------------------------------------------------
padLeft n s
| length s < n = replicate (n - length s) ' ' ++ s
| otherwise = s
padRight n s
| length s < n = s ++ replicate (n - length s) ' '
| otherwise = s
show_cf (expression, cf, note) =
let exprStr = padLeft 19 expression in
do { putStr exprStr;
putStr " => ";
if note == "" then
putStrLn (cf2string cf)
else
let cfStr = padRight 48 (cf2string cf) in
do { putStr cfStr;
putStrLn note }
}
thirteen_elevenths = r2cf 13 11
twentytwo_sevenths = r2cf 22 7
main = do {
show_cf ("golden ratio", goldenRatio, "(1 + sqrt(5))/2");
show_cf ("silver ratio", silverRatio, "(1 + sqrt(2))");
show_cf ("sqrt(2)", sqrt2, "from the module");
show_cf ("sqrt(2)", silverRatio `sub_cf` one,
"from the silver ratio");
show_cf ("sqrt(5)", sqrt5, "from the module");
show_cf ("sqrt(5)", (two `mul_cf` goldenRatio) `sub_cf` one,
"from the golden ratio");
show_cf ("13/11", thirteen_elevenths, "");
show_cf ("22/7", twentytwo_sevenths, "approximately pi");
show_cf ("13/11 + 1/2", thirteen_elevenths `add_cf` one_half, "");
show_cf ("22/7 + 1/2", twentytwo_sevenths `add_cf` one_half, "");
show_cf ("(22/7) * 1/2", twentytwo_sevenths `mul_cf` one_half, "");
show_cf ("(22/7) / 2", twentytwo_sevenths `div_cf` two, "");
show_cf ("sqrt(2) + sqrt(2)", sqrt2 `add_cf` sqrt2, "");
show_cf ("sqrt(2) - sqrt(2)", sqrt2 `sub_cf` sqrt2, "");
show_cf ("sqrt(2) * sqrt(2)", sqrt2 `mul_cf` sqrt2, "");
show_cf ("sqrt(2) / sqrt(2)", sqrt2 `div_cf` sqrt2, "");
return ()
}
----------------------------------------------------------------------
</syntaxhighlight>
{{out}}
<pre>$ ghc continued_fraction_task.hs && ./continued_fraction_task
[1 of 1] Compiling Main ( continued_fraction_task.hs, continued_fraction_task.o )
Linking continued_fraction_task ...
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,...] from the module
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] from the silver ratio
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] from the module
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] from the golden ratio
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
13/11 + 1/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7) * 1/2 => [1;1,1,3]
(22/7) / 2 => [1;1,1,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]
</pre>
=={{header|Icon}}==
{{trans|ObjectIcon}}
<syntaxhighlight lang="icon">
$define YES 1
$define NO &null
record CF (terminated, # Are there no more terms to memoize?
memo, # Memoized terms.
term_gen) # A co-expression to generate more terms.
record ng8 (a12, a1, a2, a,
b12, b1, b2, b)
global golden_ratio
global silver_ratio
global sqrt2
global frac_13_11
global frac_22_7
global one, two, three, four
global practically_infinite_number
global much_too_big_number
procedure r2cf (n, d)
return make_cf (create r2cf_generate (n, d))
end
procedure i2cf (i)
return make_cf (create r2cf_generate (i, 1))
end
procedure cf_add (x, y)
return ng8_apply (ng8 (0, 1, 1, 0, 0, 0, 0, 1), x, y)
end
procedure cf_sub (x, y)
return ng8_apply (ng8 (0, 1, -1, 0, 0, 0, 0, 1), x, y)
end
procedure cf_mul (x, y)
return ng8_apply (ng8 (1, 0, 0, 0, 0, 0, 0, 1), x, y)
end
procedure cf_div (x, y)
return ng8_apply (ng8 (0, 1, 0, 0, 0, 0, 1, 0), x, y)
end
procedure main ()
initial {
initialize_global_CF ()
}
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2")
show ("silver ratio", silver_ratio, "(1 + sqrt(2))")
show ("sqrt(2)", sqrt2, "silver ratio minus 1")
show ("13/11", frac_13_11)
show ("22/7", frac_22_7, "approximately pi")
show ("1", one)
show ("2", two)
show ("3", three)
show ("4", four)
show ("(1 + 1/sqrt(2))/2",
ng8_apply (ng8(0, 1, 0, 0, 0, 0, 2, 0),
silver_ratio, sqrt2),
"method 1")
show ("(1 + 1/sqrt(2))/2",
ng8_apply (ng8(1, 0, 0, 1, 0, 0, 0, 8),
silver_ratio, silver_ratio),
"method 2")
show ("(1 + 1/sqrt(2))/2",
cf_mul (cf_add (one, (cf_div (one, sqrt2))), r2cf (1, 2)),
"method 3")
show ("sqrt(2) + sqrt(2)", cf_add (sqrt2, sqrt2));
show ("sqrt(2) - sqrt(2)", cf_sub (sqrt2, sqrt2));
show ("sqrt(2) * sqrt(2)", cf_mul (sqrt2, sqrt2));
show ("sqrt(2) / sqrt(2)", cf_div (sqrt2, sqrt2));
end
procedure initialize_global_CF ()
golden_ratio := make_cf (create generate_constant (1))
silver_ratio := make_cf (create generate_constant (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)
sqrt2 := cf_sub (silver_ratio, one)
return
end
procedure show (expression, cf, note)
writes (right (expression, 19), " => ")
if /note then
write (cf2string (cf))
else
write (left (cf2string (cf), 48), note)
return
end
procedure make_cf (gen)
return CF (NO, [], gen)
end
# Get a continued fraction's ith term, or fail if there is no such
# term.
procedure cf_get_term (cf, i)
local j, term
if *cf.memo <= i then {
if \cf.terminated then {
fail
} else {
every j := *cf.memo to i do {
if \ (term := @cf.term_gen) then {
put (cf.memo, term)
} else {
cf.terminated := YES
fail
}
}
}
}
return cf.memo[i + 1]
end
# Generate a continued fraction's terms, with &null as "infinity".
procedure cf_generate (cf)
local i, term
i := 0
while term := cf_get_term (cf, i) do {
suspend term
i +:= 1
}
repeat suspend &null
end
# A co-expression that generates a continued fraction's terms, with
# &null as "infinity".
procedure cf2coexpression (cf)
return create cf_generate (cf)
end
# Make a human-readable string from a continued fraction.
procedure cf2string (cf, max_terms)
local s, i, done, term
static separators
initial {
separators := ["", ";", ","]
}
/max_terms := 20
s := "["
i := 0
done := NO
while /done do {
if not (term := cf_get_term (cf, i)) then {
s ||:= "]"
done := YES
} else if i = max_terms then {
s ||:= ",...]"
done := YES
} else {
s ||:= separators[min (i + 1, 3)] || term
i +:= 1
}
}
return s
end
# Make a CF that applies the operation.
procedure ng8_apply (ng8, x, y)
return make_cf (create ng8_generate (ng8, x, y))
end
# A generator that applies a bihomographic function. A return of &null
# means "infinity".
procedure ng8_generate (ng, x, y)
local xsource, ysource
local a12, a1, a2, a
local b12, b1, b2, b
local r12, r1, r2, r
local n1, n2, n
local absorb_term_from
local term
local new_a12, new_a1, new_a2, new_a
local new_b12, new_b1, new_b2, new_b
xsource := make_source (x)
ysource := make_source (y)
a12 := ng.a12; a1 := ng.a1; a2 := ng.a2; a := ng.a
b12 := ng.b12; b1 := ng.b1; b2 := ng.b2; b := ng.b
repeat {
absorb_term_from := &null
if b12 = b1 = b2 = b = 0 then {
suspend &null # "Infinity".
} else if b2 = b = 0 then {
absorb_term_from := 'x'
} else if b2 = 0 | b = 0 then {
absorb_term_from := 'y'
} else if b1 = 0 then {
absorb_term_from := 'x'
} else if b12 ~= 0 & term := (a12/b12 = a1/b1 = a2/b2 = a/b) then {
suspend infinitized (term)
r12 := a12 % b12
r1 := a1 % b1
r2 := a2 % b2
r := a % b
a12 := b12; a1 := b1; a2 := b2; a := b
b12 := r12; b1 := r1; b2 := r2; b := r
} else {
# Put numerators over a common denominator.
n1 := a1 * b2 * b
n2 := a2 * b1 * b
n := a * b1 * b2
if abs (n1 - n) > abs (n2 - n) then
absorb_term_from := 'x'
else
absorb_term_from := 'y'
}
if \absorb_term_from === 'x' then {
term := @xsource
new_a2 := a12; new_a := a1
new_b2 := b12; new_b := b1
if \term then {
new_a12 := a2 + (a12 * term)
new_a1 := a + (a1 * term)
new_b12 := b2 + (b12 * term)
new_b1 := b + (b1 * term)
if too_big (new_a12 | new_a1 | new_a2 | new_a |
new_b12 | new_b1 | new_b2 | new_b) then
{
# All further terms are forced to "infinity".
xsource := make_source (&null)
new_a12 := a12; new_a1 := a1
new_b12 := b12; new_b1 := b1
}
}
a12 := new_a12; a1 := new_a1; a2 := new_a2; a := new_a
b12 := new_b12; b1 := new_b1; b2 := new_b2; b := new_b
}
else if \absorb_term_from === 'y' then {
term := @ysource
new_a1 := a12; new_a := a2
new_b1 := b12; new_b := b2
if \term then {
new_a12 := a1 + (a12 * term)
new_a2 := a + (a2 * term)
new_b12 := b1 + (b12 * term)
new_b2 := b + (b2 * term)
if too_big (new_a12 | new_a1 | new_a2 | new_a |
new_b12 | new_b1 | new_b2 | new_b) then
{
# All further terms are forced to "infinity".
ysource := make_source (&null)
new_a12 := a12; new_a2 := a2
new_b12 := b12; new_b2 := b2
}
}
a12 := new_a12; a1 := new_a1; a2 := new_a2; a := new_a
b12 := new_b12; b1 := new_b1; b2 := new_b2; b := new_b
}
}
end
procedure make_source (cf)
local source
if /cf then {
# Generate only "infinite" terms.
source := create generate_constant (&null)
} else if type(cf) == "co-expression" then {
# Already a co-expression.
source := cf
} else {
# Use a continued fraction as a co-expression.
source := cf2coexpression (cf)
}
return source
end
procedure infinitized (term)
initial {
/practically_infinite_number := 2^64
}
if abs (term) >= abs (practically_infinite_number) then {
term := &null
}
return term
end
procedure too_big (num)
initial {
/much_too_big_number := 2^512
}
if abs (num) < abs (much_too_big_number) then fail
return num
end
procedure generate_constant (constant)
repeat suspend constant
end
procedure r2cf_generate (n, d)
local remainder
repeat {
if d = 0 then {
suspend &null
} else {
suspend (n / d)
remainder := n % d
n := d
d := remainder
}
}
end
procedure min (x, y)
return (if x < y then x else y)
end
</syntaxhighlight>
{{out}}
<pre>$ icon bivariate_continued_fraction_task_Icon.icn
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,...] silver ratio minus 1
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
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]
</pre>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.List;
public final class ContinuedFractionArithmeticG2 {
public static void main(String[] aArgs) {
test("[3; 7] + [0; 2]", new NG( new NG8(0, 1, 1, 0, 0, 0, 0, 1), new R2cf(1, 2), new R2cf(22, 7) ),
new NG( new NG4(2, 1, 0, 2), new R2cf(22, 7) ));
test("[1; 5, 2] * [3; 7]", new NG( new NG8(1, 0, 0, 0, 0, 0, 0, 1), new R2cf(13, 11), new R2cf(22, 7) ),
new R2cf(286, 77) );
test("[1; 5, 2] - [3; 7]", new NG( new NG8(0, 1, -1, 0, 0, 0, 0, 1), new R2cf(13, 11), new R2cf(22, 7) ),
new R2cf(-151, 77) );
test("Divide [] by [3; 7]",
new NG( new NG8(0, 1, 0, 0, 0, 0, 1, 0), new R2cf(22 * 22, 7 * 7), new R2cf(22,7)) );
test("([0; 3, 2] + [1; 5, 2]) * ([0; 3, 2] - [1; 5, 2])",
new NG( new NG8(1, 0, 0, 0, 0, 0, 0, 1),
new NG( new NG8(0, 1, 1, 0, 0, 0, 0, 1),
new R2cf(2, 7), new R2cf(13, 11)),
new NG( new NG8(0, 1, -1, 0, 0, 0, 0, 1), new R2cf(2, 7), new R2cf(13, 11) ) ),
new R2cf(-7797, 5929) );
}
private static void test(String aDescription, ContinuedFraction... aFractions) {
System.out.println("Testing: " + aDescription);
for ( ContinuedFraction fraction : aFractions ) {
while ( fraction.hasMoreTerms() ) {
System.out.print(fraction.nextTerm() + " ");
}
System.out.println();
}
System.out.println();
}
private static abstract class MatrixNG {
protected abstract void consumeTerm();
protected abstract void consumeTerm(int aN);
protected abstract boolean needsTerm();
protected int configuration = 0;
protected int currentTerm = 0;
protected boolean hasTerm = false;
}
private static class NG4 extends MatrixNG {
public NG4(int aA1, int aA, int aB1, int aB) {
a1 = aA1; a = aA; b1 = aB1; b = aB;
}
public void consumeTerm() {
a = a1;
b = b1;
}
public void consumeTerm(int aN) {
int temp = a; a = a1; a1 = temp + a1 * aN;
temp = b; b = b1; b1 = temp + b1 * aN;
}
public boolean needsTerm() {
if ( b1 == 0 && b == 0 ) {
return false;
}
if ( b1 == 0 || b == 0 ) {
return true;
}
currentTerm = a / b;
if ( currentTerm == a1 / b1 ) {
int temp = a; a = b; b = temp - b * currentTerm;
temp = a1; a1 = b1; b1 = temp - b1 * currentTerm;
hasTerm = true;
return false;
}
return true;
}
private int a1, a, b1, b;
}
private static class NG8 extends MatrixNG {
public NG8(int aA12, int aA1, int aA2, int aA, int aB12, int aB1, int aB2, int aB) {
a12 = aA12; a1 = aA1; a2 = aA2; a = aA; b12 = aB12; b1 = aB1; b2 = aB2; b = aB;
}
public void consumeTerm() {
if ( configuration == 0 ) {
a = a1; a2 = a12;
b = b1; b2 = b12;
} else {
a = a2; a1 = a12;
b = b2; b1 = b12;
}
}
public void consumeTerm(int aN) {
if ( configuration == 0 ) {
int temp = a; a = a1; a1 = temp + a1 * aN;
temp = a2; a2 = a12; a12 = temp + a12 * aN;
temp = b; b = b1; b1 = temp + b1 * aN;
temp = b2; b2 = b12; b12 = temp + b12 * aN;
} else {
int temp = a; a = a2; a2 = temp + a2 * aN;
temp = a1; a1 = a12; a12 = temp + a12 * aN;
temp = b; b = b2; b2 = temp + b2 * aN;
temp = b1; b1 = b12; b12 = temp + b12 * aN;
}
}
public boolean needsTerm() {
if ( b1 == 0 && b == 0 && b2 == 0 && b12 == 0 ) {
return false;
}
if ( b == 0 ) {
configuration = ( b2 == 0 ) ? 0 : 1;
return true;
}
ab = (double) a / b;
if ( b2 == 0 ) {
configuration = 1;
return true;
}
a2b2 = (double) a2 / b2;
if ( b1 == 0 ) {
configuration = 0;
return true;
}
a1b1 = (double) a1 / b1;
if ( b12 == 0 ) {
configuration = setConfiguration();
return true;
}
a12b12 = (double) a12 / b12;
currentTerm = (int) ab;
if ( currentTerm == (int) a1b1 && currentTerm == (int) a2b2 && currentTerm == (int) a12b12 ) {
int temp = a; a = b; b = temp - b * currentTerm;
temp = a1; a1 = b1; b1 = temp - b1 * currentTerm;
temp = a2; a2 = b2; b2 = temp - b2 * currentTerm;
temp = a12; a12 = b12; b12 = temp - b12 * currentTerm;
hasTerm = true;
return false;
}
configuration = setConfiguration();
return true;
}
private int setConfiguration() {
return ( Math.abs(a1b1 - ab) > Math.abs(a2b2 - ab) ) ? 0 : 1;
}
private int a12, a1, a2, a, b12, b1, b2, b;
private double ab, a1b1, a2b2, a12b12;
}
private static interface ContinuedFraction {
public boolean hasMoreTerms();
public int nextTerm();
}
private static class R2cf implements ContinuedFraction {
public R2cf(int aN1, int aN2) {
n1 = aN1; n2 = aN2;
}
public boolean hasMoreTerms() {
return Math.abs(n2) > 0;
}
public int nextTerm() {
final int term = n1 / n2;
final int temp = n2;
n2 = n1 - term * n2;
n1 = temp;
return term;
}
private int n1, n2;
}
private static class NG implements ContinuedFraction {
public NG(NG4 aNG, ContinuedFraction aCF) {
matrixNG = aNG;
cf.add(aCF);
}
public NG(NG8 aNG, ContinuedFraction aCF1, ContinuedFraction aCF2) {
matrixNG = aNG;
cf.add(aCF1); cf.add(aCF2);
}
public boolean hasMoreTerms() {
while ( matrixNG.needsTerm() ) {
if ( cf.get(matrixNG.configuration).hasMoreTerms() ) {
matrixNG.consumeTerm(cf.get(matrixNG.configuration).nextTerm());
} else {
matrixNG.consumeTerm();
}
}
return matrixNG.hasTerm;
}
public int nextTerm() {
matrixNG.hasTerm = false;
return matrixNG.currentTerm;
}
private MatrixNG matrixNG;
private List<ContinuedFraction> cf = new ArrayList<ContinuedFraction>();
}
}
</syntaxhighlight>
{{ out }}
<pre>
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
</pre>
Line 3,262 ⟶ 7,456:
-1 -3 -5 -1 -2 -1 -26 -3
-1 -3 -5 -1 -2 -1 -26 -3
</pre>
=={{header|Mercury}}==
{{works with|Mercury|22.01.1}}
{{trans|Standard ML}}
This program was written with reference to the Standard ML, but really is a different kind of implementation: the continued fractions are represented as lazy lists.
The program is not very fast, but this might be due mainly to the <code>integer</code> type in the Mercury standard library not being very fast. I do not know. If so, an interface to the GNU Multiple Precision Arithmetic Library might speed things up quite a bit.
The program comes in two source files. The main program goes in file <code>continued_fraction_task.m</code>:
<syntaxhighlight lang="Mercury">
%%% -*- mode: mercury; prolog-indent-width: 2; -*-
%%%
%%% A program in two files:
%%% continued_fraction_task.m (this file)
%%% continued_fraction.m (the continued_fraction module)
%%%
%%% Compile with:
%%% mmc --make --use-subdirs continued_fraction_task
%%%
:- module continued_fraction_task.
:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.
:- implementation.
:- import_module continued_fraction.
:- import_module integer.
:- import_module rational.
:- import_module string.
:- pred show(string::in, continued_fraction::in, string::in,
io::di, io::uo) is det.
:- pred show(string::in, continued_fraction::in,
io::di, io::uo) is det.
show(Expression, CF, Note, !IO) :-
pad_left(Expression, (' '), 19, Expr1),
print(Expr1, !IO),
print(" => ", !IO),
(if (Note = "")
then (print(to_string(CF), !IO),
nl(!IO))
else (pad_right(to_string(CF), (' '), 48, CF1_Str),
print(CF1_Str, !IO),
print(Note, !IO),
nl(!IO))).
show(Expression, CF, !IO) :- show(Expression, CF, "", !IO).
:- func thirteen_elevenths = continued_fraction.
thirteen_elevenths = from_rational(rational(13, 11)).
:- func twentytwo_sevenths = continued_fraction.
twentytwo_sevenths = from_rational(rational(22, 7)).
main(!IO) :-
show("golden ratio", golden_ratio, "(1 + sqrt(5))/2", !IO),
show("silver ratio", silver_ratio, "(1 + sqrt(2))", !IO),
show("sqrt(2)", sqrt2, "from the module", !IO),
show("sqrt(2)", silver_ratio - one, "from the silver ratio", !IO),
show("sqrt(5)", sqrt5, "from the module", !IO),
show("sqrt(5)", (two * golden_ratio) - one, "from the golden ratio", !IO),
show("13/11", thirteen_elevenths, !IO),
show("22/7", twentytwo_sevenths, "approximately pi", !IO),
show("13/11 + 1/2", thirteen_elevenths + one_half, !IO),
show("22/7 + 1/2", twentytwo_sevenths + one_half, !IO),
show("(22/7) * 1/2", twentytwo_sevenths * one_half, !IO),
show("(22/7) / 2", twentytwo_sevenths / two, !IO),
show("sqrt(2) + sqrt(2)", sqrt2 + sqrt2, !IO),
show("sqrt(2) - sqrt(2)", sqrt2 - sqrt2, !IO),
show("sqrt(2) * sqrt(2)", sqrt2 * sqrt2, !IO),
show("sqrt(2) / sqrt(2)", sqrt2 / sqrt2, !IO),
true.
:- end_module continued_fraction_task.
</syntaxhighlight>
The <code>continued_fraction</code> module source goes in file <code>continued_fraction.m</code>:
<syntaxhighlight lang="Mercury">
%%% -*- mode: mercury; prolog-indent-width: 2; -*-
:- module continued_fraction.
:- interface.
:- import_module int.
:- import_module integer.
:- import_module lazy.
:- import_module rational.
%% A continued fraction is a kind of lazy list. The list is always
%% infinitely long. However, you need not consider terms that come
%% after an infinite term.
:- type continued_fraction
---> continued_fraction(lazy(node)).
:- type node
---> cons(term, continued_fraction).
:- type term
---> infinite_term
; integer_term(integer).
%% ng8 = {A12, A1, A2, A, B12, B1, B2, B}.
:- type ng8 == {integer, integer, integer, integer,
integer, integer, integer, integer}.
%% Get a human-readable string. The second form takes a "MaxTerms"
%% argument. The first form uses a default value for MaxTerms.
:- func to_string(continued_fraction) = string.
:- func to_string(int, continued_fraction) = string.
%% Make a term from a regular int.
:- func int_term(int) = term.
%% A "continued fraction" with only infinite terms.
:- func infinite_continued_fraction = continued_fraction.
%% A continued fraction whose term repeats infinitely.
:- func repeating_term(term) = continued_fraction.
%% A continued fraction representing an integer.
:- func from_integer(integer) = continued_fraction.
:- func from_int(int) = continued_fraction.
%% A continued fraction representing a rational number.
:- func from_rational(rational) = continued_fraction.
%% A continued fraction that is a bihomographic function of two other
%% continued fractions.
:- func apply_ng8(ng8, continued_fraction,
continued_fraction) = continued_fraction.
:- func '+'(continued_fraction,
continued_fraction) = continued_fraction.
:- func '-'(continued_fraction,
continued_fraction) = continued_fraction.
:- func '*'(continued_fraction,
continued_fraction) = continued_fraction.
:- func '/'(continued_fraction,
continued_fraction) = continued_fraction.
%% Miscellaneous continued fractions.
:- func zero = continued_fraction.
:- func one = continued_fraction.
:- func two = continued_fraction.
:- func three = continued_fraction.
:- func four = continued_fraction.
%%
:- func one_fourth = continued_fraction.
:- func one_third = continued_fraction.
:- func one_half = continued_fraction.
:- func two_thirds = continued_fraction.
:- func three_fourths = continued_fraction.
%%
:- func golden_ratio = continued_fraction. % (1 + sqrt(5))/2
:- func silver_ratio = continued_fraction. % (1 + sqrt(2))
:- func sqrt2 = continued_fraction. % The square root of two.
:- func sqrt5 = continued_fraction. % The square root of five.
:- implementation.
:- import_module string.
%%--------------------------------------------------------------------
to_string(CF) = to_string(20, CF).
to_string(MaxTerms, CF) = S :-
to_string(MaxTerms, CF, 0, "[", S).
:- pred to_string(int::in, continued_fraction::in, int::in,
string::in, string::out) is det.
to_string(MaxTerms, CF, I, S0, S) :-
CF = continued_fraction(Node),
force(Node) = cons(Term, CF1),
(if (Term = integer_term(N))
then (if (I = MaxTerms) then (S = S0 ++ ",...]")
else (TermStr = (integer.to_string(N)),
Separator = (if (I = 0) then ""
else if (I = 1) then ";"
else ","),
to_string(MaxTerms, CF1, I + 1,
S0 ++ Separator ++ TermStr, S)))
else (S = S0 ++ "]")).
%%--------------------------------------------------------------------
int_term(I) = integer_term(integer(I)).
infinite_continued_fraction = CF :-
CF = repeating_term(infinite_term).
repeating_term(T) = CF :-
CF = continued_fraction(Node),
Node = delay((func) = cons(T, repeating_term(T))).
from_integer(I) = CF :-
CF = continued_fraction(Node),
Node = delay((func) = cons(integer_term(I), infinite_continued_fraction)).
from_int(I) = from_integer(integer(I)).
from_rational(R) = CF :-
N = numer(R),
D = denom(R),
CF = from_rational_integers(N, D).
:- func from_rational_integers(integer, integer) = continued_fraction.
from_rational_integers(N, D) = CF :-
if (D = zero) then (CF = infinite_continued_fraction)
else (divide_with_rem(N, D, Q, R),
CF = continued_fraction(
delay((func) = cons(integer_term(Q),
from_rational_integers(D, R)))
)).
%%--------------------------------------------------------------------
(X : continued_fraction) + (Y : continued_fraction) = (
apply_ng8({zero, one, one, zero, zero, zero, zero, one}, X, Y)
).
(X : continued_fraction) - (Y : continued_fraction) = (
apply_ng8({zero, one, negative_one, zero, zero, zero, zero, one}, X, Y)
).
(X : continued_fraction) * (Y : continued_fraction) = (
apply_ng8({one, zero, zero, zero, zero, zero, zero, one}, X, Y)
).
(X : continued_fraction) / (Y : continued_fraction) = (
apply_ng8({zero, one, zero, zero, zero, zero, one, zero}, X, Y)
).
apply_ng8(NG, X, Y) = CF :-
NG = {A12, A1, A2, A, B12, B1, B2, B},
(if iseqz(B12, B1, B2, B) then (
%% There are no more finite terms to output.
CF = infinite_continued_fraction
)
else if iseqz(B2, B) then (
absorb_x_term(NG, NG1, X, X1, Y, Y1),
CF = apply_ng8(NG1, X1, Y1)
)
else if at_least_one_iseqz(B2, B) then (
absorb_y_term(NG, NG1, X, X1, Y, Y1),
CF = apply_ng8(NG1, X1, Y1)
)
else if iseqz(B1) then (
absorb_x_term(NG, NG1, X, X1, Y, Y1),
CF = apply_ng8(NG1, X1, Y1)
)
else (
maybe_divide(A12, B12, Q12, R12),
maybe_divide(A1, B1, Q1, R1),
maybe_divide(A2, B2, Q2, R2),
maybe_divide(A, B, Q, R),
(if (not (iseqz(B12)), Q = Q12, Q = Q1, Q = Q2)
then (
%% Output a term.
if (integer_exceeds_infinitizing_threshold(Q))
then (CF = infinite_continued_fraction)
else (NG1 = {B12, B1, B2, B, R12, R1, R2, R},
CF = continued_fraction(
delay((func) = cons(integer_term(Q),
apply_ng8(NG1, X, Y)))
))
)
else (
%% Put A1, A2, and A over a common denominator and compare some
%% magnitudes.
N1 = A1 * B2 * B,
N2 = A2 * B1 * B,
N = A * B1 * B2,
(if (abs(N1 - N) > abs(N2 - N))
then (absorb_x_term(NG, NG1, X, X1, Y, Y1),
CF = apply_ng8(NG1, X1, Y1))
else (absorb_y_term(NG, NG1, X, X1, Y, Y1),
CF = apply_ng8(NG1, X1, Y1)))
))
)).
:- pred absorb_x_term(ng8::in, ng8::out,
continued_fraction::in, continued_fraction::out,
continued_fraction::in, continued_fraction::out)
is det.
absorb_x_term(!NG, !X, !Y) :-
(!.NG) = {A12, A1, A2, A, B12, B1, B2, B},
(!.X) = continued_fraction(XNode),
force(XNode) = cons(XTerm, X1),
(if (XTerm = integer_term(N))
then (New_NG = {A2 + (A12 * N), A + (A1 * N), A12, A1,
B2 + (B12 * N), B + (B1 * N), B12, B1},
(if (ng8_exceeds_processing_threshold(New_NG))
then (
%% Pretend we have reached an infinite term.
!:NG = {A12, A1, A12, A1, B12, B1, B12, B1},
!:X = infinite_continued_fraction
)
else (!:NG = New_NG, !:X = X1)))
else (!:NG = {A12, A1, A12, A1, B12, B1, B12, B1},
!:X = infinite_continued_fraction)).
:- pred absorb_y_term(ng8::in, ng8::out,
continued_fraction::in, continued_fraction::out,
continued_fraction::in, continued_fraction::out)
is det.
absorb_y_term(!NG, !X, !Y) :-
(!.NG) = {A12, A1, A2, A, B12, B1, B2, B},
(!.Y) = continued_fraction(YNode),
force(YNode) = cons(YTerm, Y1),
(if (YTerm = integer_term(N))
then (New_NG = {A1 + (A12 * N), A12, A + (A2 * N), A2,
B1 + (B12 * N), B12, B + (B2 * N), B2},
(if (ng8_exceeds_processing_threshold(New_NG))
then (
%% Pretend we have reached an infinite term.
!:NG = {A12, A12, A2, A2, B12, B12, B2, B2},
!:Y = infinite_continued_fraction
)
else (!:NG = New_NG, !:Y = Y1)))
else (!:NG = {A12, A12, A2, A2, B12, B12, B2, B2},
!:Y = infinite_continued_fraction)).
:- pred ng8_exceeds_processing_threshold(ng8::in) is semidet.
:- pred integer_exceeds_processing_threshold(integer::in) is semidet.
ng8_exceeds_processing_threshold({A12, A1, A2, A,
B12, B1, B2, B}) :-
(integer_exceeds_processing_threshold(A12) ;
integer_exceeds_processing_threshold(A1) ;
integer_exceeds_processing_threshold(A2) ;
integer_exceeds_processing_threshold(A) ;
integer_exceeds_processing_threshold(B12) ;
integer_exceeds_processing_threshold(B1) ;
integer_exceeds_processing_threshold(B2) ;
integer_exceeds_processing_threshold(B)).
integer_exceeds_processing_threshold(Integer) :-
abs(Integer) >= pow(two, integer(512)).
:- pred integer_exceeds_infinitizing_threshold(integer::in) is semidet.
integer_exceeds_infinitizing_threshold(Integer) :-
abs(Integer) >= pow(two, integer(64)).
:- pred maybe_divide(integer::in, integer::in,
integer::out, integer::out) is det.
maybe_divide(N, D, Q, R) :-
if iseqz(D) then (Q = zero, R = zero)
else divide_with_rem(N, D, Q, R).
:- pred iseqz(integer::in) is semidet.
:- pred iseqz(integer::in, integer::in) is semidet.
:- pred iseqz(integer::in, integer::in,
integer::in, integer::in) is semidet.
iseqz(Integer) :- Integer = zero.
iseqz(A, B) :- iseqz(A), iseqz(B).
iseqz(A, B, C, D) :- iseqz(A), iseqz(B), iseqz(C), iseqz(D).
:- pred at_least_one_iseqz(integer::in, integer::in) is semidet.
at_least_one_iseqz(A, B) :- (A = zero; B = zero).
%%--------------------------------------------------------------------
:- func two_plus_sqrt5 = continued_fraction.
two_plus_sqrt5 = repeating_term(int_term(4)).
zero = from_int(0).
one = from_int(1).
two = from_int(2).
three = from_int(3).
four = from_int(4).
one_fourth = from_rational(rational(1, 4)).
one_third = from_rational(rational(1, 3)).
one_half = from_rational(rational(1, 2)).
two_thirds = from_rational(rational(2, 3)).
three_fourths = from_rational(rational(3, 4)).
golden_ratio = repeating_term(int_term(1)).
silver_ratio = repeating_term(int_term(2)).
sqrt2 = continued_fraction(delay((func) = cons(int_term(1), silver_ratio))).
sqrt5 = continued_fraction(delay((func) = cons(int_term(2), two_plus_sqrt5))).
%%--------------------------------------------------------------------
:- end_module continued_fraction.
</syntaxhighlight>
{{out}}
<pre>$ mmc --make --use-subdirs continued_fraction_task && ./continued_fraction_task
Making Mercury/int3s/continued_fraction_task.int3
Making Mercury/int3s/continued_fraction.int3
Making Mercury/ints/continued_fraction.int
Making Mercury/ints/continued_fraction_task.int
Making Mercury/cs/continued_fraction.c
Making Mercury/cs/continued_fraction_task.c
Making Mercury/os/continued_fraction.o
Making Mercury/os/continued_fraction_task.o
Making continued_fraction_task
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,...] from the module
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] from the silver ratio
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] from the module
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] from the golden ratio
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
13/11 + 1/2 => [1;1,2,7]
22/7 + 1/2 => [3;1,1,1,4]
(22/7) * 1/2 => [1;1,1,3]
(22/7) / 2 => [1;1,1,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]
</pre>
Line 3,503 ⟶ 8,112:
-1 -3 -5 -1 -2 -1 -26 -3
-1 -3 -5 -1 -2 -1 -26 -3 </pre>
=={{header|ObjectIcon}}==
{{trans|Scheme}}
Where the Scheme uses thunks, the Object Icon uses co-expressions.
<syntaxhighlight lang="ObjectIcon">
# -*- ObjectIcon -*-
import io
global golden_ratio
global silver_ratio
global sqrt2
global frac_13_11
global frac_22_7
global one, two, three, four
procedure r2cf (n, d)
return CF (create r2cf_generate (n, d))
end
procedure i2cf (i)
return CF (create r2cf_generate (i, 1))
end
procedure cf_add (x, y)
return NG8(0, 1, 1, 0, 0, 0, 0, 1).apply(x, y)
end
procedure cf_sub (x, y)
return NG8(0, 1, -1, 0, 0, 0, 0, 1).apply(x, y)
end
procedure cf_mul (x, y)
return NG8(1, 0, 0, 0, 0, 0, 0, 1).apply(x, y)
end
procedure cf_div (x, y)
return NG8(0, 1, 0, 0, 0, 0, 1, 0).apply(x, y)
end
procedure main ()
initial initialize_global_CF ()
show ("golden ratio", golden_ratio, "(1 + sqrt(5))/2")
show ("silver ratio", silver_ratio, "(1 + sqrt(2))")
show ("sqrt(2)", sqrt2, "silver ratio minus 1")
show ("13/11", frac_13_11)
show ("22/7", frac_22_7, "approximately pi")
show ("1", one)
show ("2", two)
show ("3", three)
show ("4", four)
show ("(1 + 1/sqrt(2))/2",
NG8(0, 1, 0, 0, 0, 0, 2, 0).apply(silver_ratio, sqrt2),
"method 1")
show ("(1 + 1/sqrt(2))/2",
NG8(1, 0, 0, 1, 0, 0, 0, 8).apply(silver_ratio, silver_ratio),
"method 2")
show ("(1 + 1/sqrt(2))/2",
cf_mul (cf_add (one, (cf_div (one, sqrt2))), r2cf (1, 2)),
"method 3")
show ("sqrt(2) + sqrt(2)", cf_add (sqrt2, sqrt2));
show ("sqrt(2) - sqrt(2)", cf_sub (sqrt2, sqrt2));
show ("sqrt(2) * sqrt(2)", cf_mul (sqrt2, sqrt2));
show ("sqrt(2) / sqrt(2)", cf_div (sqrt2, sqrt2));
end
procedure initialize_global_CF ()
/golden_ratio := CF (create generate_constant (1))
/silver_ratio := CF (create generate_constant (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)
/sqrt2 := cf_sub (silver_ratio, one)
return
end
procedure show (expression, cf, note)
io.writes (right (expression, 19), " => ")
if /note then
io.write (cf.to_string())
else
io.write (left (cf.to_string(), 48), note)
return
end
class CF () # Continued fraction.
private terminated # Are there no more terms to memoize?
private memo # Memoized terms.
private term_gen # A co-expression to generate more terms.
public static default_max_terms
public new (gen)
terminated := &no
memo := []
term_gen := gen
return
end
public get_term (i) # Get the ith term (or fail if there is none).
local j, term
if *memo <= i then
{
if \terminated then
fail
else
{
every j := *memo to i do
{
if \ (term := @term_gen) then
put (memo, term)
else {
terminated := &yes
fail
}
}
}
}
return memo[i + 1]
end
public generate () # Generate the terms, with &null as "infinity".
local i, term
i := 0
while term := get_term (i) do
{
suspend term
i +:= 1
}
repeat suspend &null
end
public to_coexpression () # Co-expression that generates terms.
return create generate ()
end
public to_string (max_terms) # Make a human-readable string.
local s, i, done, term
static separators
initial separators := ["", ";", ","]
/max_terms := (\default_max_terms | 20)
s := "["
i := 0
done := &no
while /done do
{
if not (term := get_term (i)) then
{
s ||:= "]"
done := &yes
}
else if i = max_terms then
{
s ||:= ",...]"
done := &yes
}
else
{
s ||:= separators[min (i + 1, 3)] || term
i +:= 1
}
}
return s
end
end # end class CF
class NG8 () # A bihomographic function.
private na12, na1, na2, na
private nb12, nb1, nb2, nb
public static practically_infinite_number
public static much_too_big_number
public new (a12, a1, a2, a, b12, b1, b2, b)
na12 := a12
na1 := a1
na2 := a2
na := a
nb12 := b12
nb1 := b1
nb2 := b2
nb := b
return
end
public apply (x, y) # Make a CF that applies the operation.
return CF (create generate (x, y))
end
public generate (x, y) # Apply the operation to generate terms.
local xsource, ysource
local a12, a1, a2, a
local b12, b1, b2, b
local r12, r1, r2, r
local n1, n2, n
local absorb_term_from
local term
local new_a12, new_a1, new_a2, new_a
local new_b12, new_b1, new_b2, new_b
xsource := make_source (x)
ysource := make_source (y)
a12 := na12; a1 := na1; a2 := na2; a := na
b12 := nb12; b1 := nb1; b2 := nb2; b := nb
repeat
{
absorb_term_from := &null
if b12 = b1 = b2 = b = 0 then
suspend &null # "Infinity".
else if b2 = b = 0 then
absorb_term_from := 'x'
else if b2 = 0 | b = 0 then
absorb_term_from := 'y'
else if b1 = 0 then
absorb_term_from := 'x'
else if b12 ~= 0 & term := (a12/b12 = a1/b1 = a2/b2 = a/b) then
{
suspend infinitized (term)
r12 := a12 % b12
r1 := a1 % b1
r2 := a2 % b2
r := a % b
a12 := b12; a1 := b1; a2 := b2; a := b
b12 := r12; b1 := r1; b2 := r2; b := r
}
else
{
# Put numerators over a common denominator.
n1 := a1 * b2 * b
n2 := a2 * b1 * b
n := a * b1 * b2
if abs (n1 - n) > abs (n2 - n) then
absorb_term_from := 'x'
else
absorb_term_from := 'y'
}
if \absorb_term_from === 'x' then
{
term := @xsource
new_a2 := a12; new_a := a1
new_b2 := b12; new_b := b1
if \term then
{
new_a12 := a2 + (a12 * term)
new_a1 := a + (a1 * term)
new_b12 := b2 + (b12 * term)
new_b1 := b + (b1 * term)
if too_big (new_a12 | new_a1 | new_a2 | new_a |
new_b12 | new_b1 | new_b2 | new_b) then
{
# All further terms are forced to "infinity".
xsource := make_source (&null)
new_a12 := a12; new_a1 := a1
new_b12 := b12; new_b1 := b1
}
}
a12 := new_a12; a1 := new_a1; a2 := new_a2; a := new_a
b12 := new_b12; b1 := new_b1; b2 := new_b2; b := new_b
}
else if \absorb_term_from === 'y' then
{
term := @ysource
new_a1 := a12; new_a := a2
new_b1 := b12; new_b := b2
if \term then
{
new_a12 := a1 + (a12 * term)
new_a2 := a + (a2 * term)
new_b12 := b1 + (b12 * term)
new_b2 := b + (b2 * term)
if too_big (new_a12 | new_a1 | new_a2 | new_a |
new_b12 | new_b1 | new_b2 | new_b) then
{
# All further terms are forced to "infinity".
ysource := make_source (&null)
new_a12 := a12; new_a2 := a2
new_b12 := b12; new_b2 := b2
}
}
a12 := new_a12; a1 := new_a1; a2 := new_a2; a := new_a
b12 := new_b12; b1 := new_b1; b2 := new_b2; b := new_b
}
}
end
private make_source (cf)
local source
if /cf then
# Generate only "infinite" terms.
source := create generate_constant (&null)
else if type(cf) == "co-expression" then
# Already a co-expression.
source := cf
else
# Use a continued fraction as a co-expression.
source := cf.to_coexpression ()
return source
end
private infinitized (term)
initial /practically_infinite_number := 2^64
if abs (term) >= abs (practically_infinite_number) then
term := &null
return term
end
private too_big (num)
initial /much_too_big_number := 2^512
if abs (num) < abs (much_too_big_number) then fail
return num
end
end # end class NG8
procedure generate_constant (constant)
repeat suspend constant
end
procedure r2cf_generate (n, d)
local remainder
repeat
{
if d = 0 then
suspend &null
else
{
suspend (n / d)
remainder := n % d
n := d
d := remainder
}
}
end
</syntaxhighlight>
{{out}}
<pre>$ oiscript bivariate_continued_fraction_task_OI.icn
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,...] silver ratio minus 1
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
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]
</pre>
=={{header|OCaml}}==
{{trans|Standard ML}}
{{libheader|Zarith}}
To facilitate comparison of OCaml and Standard ML, I follow the Standard ML implementation closely.
You will need the Zarith module, which is a commonly used interface to GNU Multiple Precision.
<syntaxhighlight lang="ocaml">
(* Compile and run with (for instance):
ocamlfind ocamlc -linkpkg -package zarith bivariate_continued_fraction_task.ml && ./a.out
*)
module type CONTINUED_FRACTION =
sig
(* A term_generator thunk generates terms, which a t structure
memoizes. *)
type t
type term_generator = unit -> Z.t option
type ng8 = Z.t * Z.t * Z.t * Z.t * Z.t * Z.t * Z.t * Z.t
(* Create a continued fraction. *)
val make : term_generator -> t
(* Does the indexed term exist? *)
val exists : t -> int -> bool
(* Retrieve the indexed term. *)
val get : t -> int -> Z.t
(* Use a t as a term_generator thunk. *)
val to_term_generator : t -> term_generator
(* Get a human-readable string. *)
val default_max_terms : int ref
val to_string : ?max_terms:int -> t -> string
(* Create a continued fraction representing an integer. *)
val of_bigint : Z.t -> t
val of_int : int -> t
(* Create a continued fraction representing a rational number. *)
val of_bigrat : Q.t -> t
val of_bigints : Z.t -> Z.t -> t
val of_ints : int -> int -> t
(* Create a continued fraction that has one term repeated
forever. *)
val constant_term_of_bigint : Z.t -> t
val constant_term_of_int : int -> t
(* Create a continued fraction by arithmetic. (I have not bothered
here to implement ng4, although one likely would wish to have
ng4 as well.) *)
val ng8_of_ints : (int * int * int * int
* int * int * int * int) -> ng8
val ng8_stopping_processing_threshold : Z.t ref
val ng8_infinitization_threshold : Z.t ref
val apply_ng8 : ng8 -> t -> t -> t
val ( + ) : t -> t -> t
val ( - ) : t -> t -> t
val ( * ) : t -> t -> t
val ( / ) : t -> t -> t
val ( ~- ) : t -> t
val pow : t -> int -> t
(* Miscellaneous continued fractions. *)
val zero : t
val one : t
val two : t
val three : t
val four : t
val one_fourth : t
val one_third : t
val one_half : t
val two_thirds : t
val three_fourths : t
val golden_ratio : t
val silver_ratio : t
val sqrt2 : t
val sqrt5 : t
end
module Continued_fraction : CONTINUED_FRACTION =
struct
type term_generator = unit -> Z.t option
type record_t =
{
terminated : bool; (* Is the generator exhausted? *)
memo_count : int; (* How many terms are memoized? *)
memo : Z.t Array.t; (* Memoized terms. *)
generate : term_generator; (* The source of terms. *)
}
type t = record_t ref
type ng8 = Z.t * Z.t * Z.t * Z.t * Z.t * Z.t * Z.t * Z.t
let make generator =
ref { terminated = false;
memo_count = 0;
memo = Array.make 32 Z.zero;
generate = generator }
let resize_if_necessary (cf : t) i =
let record = !cf in
if record.terminated then
()
else if i < Array.length record.memo then
()
else
let new_size = 2 * (i + 1) in
let new_memo = Array.make new_size Z.zero in
let rec copy_terms i =
if i = record.memo_count then
()
else
let term = record.memo.(i) in
new_memo.(i) <- term;
copy_terms (i + 1)
in
let new_record = { record with memo = new_memo } in
copy_terms 0;
cf := new_record
let rec update_terms (cf : t) i =
let record = !cf in
if record.terminated then
()
else if i < record.memo_count then
()
else
match record.generate () with
| None ->
let new_record = { record with terminated = true} in
cf := new_record
| Some term ->
let () = record.memo.(record.memo_count) <- term in
let new_record = { record with memo_count =
succ record.memo_count } in
cf := new_record;
update_terms cf i
let exists (cf : t) i =
resize_if_necessary cf i;
update_terms cf i;
i < (!cf).memo_count
let get (cf : t) i =
let record = !cf in
if record.memo_count <= i then
raise (Invalid_argument
"Continued_fraction.get:out_of_bounds")
else
record.memo.(i)
let to_term_generator (cf : t) =
let i : int ref = ref 0 in
fun () -> let j = !i in
if exists cf j then
begin
i := succ j;
Some (get cf j)
end
else
None
let default_max_terms = ref 20
let to_string ?(max_terms = !default_max_terms) (cf : t) =
if max_terms < 1 then
raise (Invalid_argument
"Continued_fraction.to_string:max_terms_out_of_bounds")
else
let rec loop i accum =
if not (exists cf i) then
accum ^ "]"
else if i = max_terms then
accum ^ ",...]"
else
let separator = if i = 0 then
""
else if i = 1 then
";"
else
"," in
let term_string = Z.to_string (get cf i) in
loop (i + 1) (accum ^ separator ^ term_string)
in
loop 0 "["
let of_bigint i =
let finished = ref false in
make (fun () -> (if !finished then
None
else
begin
finished := true;
Some i
end))
let of_int i = of_bigint (Z.of_int i)
let i2cf = of_int
let constant_term_of_bigint i = make (fun () -> Some i)
let constant_term_of_int i = constant_term_of_bigint (Z.of_int i)
let constant_term_cf = constant_term_of_int
let of_bigrat ratnum =
let ratnum = ref ratnum in
make (fun () -> (if Q.is_real !ratnum then
let (n, d) = (Q.num !ratnum,
Q.den !ratnum) in
let (q, r) = Z.ediv_rem n d in
begin
ratnum := { num = d; den = r };
Some q
end
else
None))
let of_bigints n d = of_bigrat { num = n; den = d }
let of_ints n d = of_bigints (Z.of_int n) (Z.of_int d)
let ng8_of_ints (a12, a1, a2, a, b12, b1, b2, b) =
let f = Z.of_int in
(f a12, f a1, f a2, f a, f b12, f b1, f b2, f b)
let ng8_stopping_processing_threshold = ref Z.(one lsl 512)
let ng8_infinitization_threshold = ref Z.(one lsl 64)
let too_big term =
Z.(abs (term) >= abs (!ng8_stopping_processing_threshold))
let any_too_big (a, b, c, d, e, f, g, h) =
too_big (a) || too_big (b)
|| too_big (c) || too_big (d)
|| too_big (e) || too_big (f)
|| too_big (g) || too_big (h)
let infinitize term =
if Z.(abs (term) >= abs (!ng8_infinitization_threshold)) then
None
else
Some term
let no_terms_source () = None
let equal_zero number = (Z.sign number = 0)
let divide a b =
if equal_zero b then
(Z.zero, Z.zero)
else
Z.ediv_rem a b
let apply_ng8 (ng : ng8) =
fun (x : t) (y : t) ->
begin
let ng = ref ng
and xsource = ref (to_term_generator x)
and ysource = ref (to_term_generator y) in
let all_b_are_zero () =
let (_, _, _, _, b12, b1, b2, b) = !ng in
equal_zero b && equal_zero b2
&& equal_zero b1 && equal_zero b12
in
let all_four_equal (a, b, c, d) =
a = b && a = c && a = d
in
let absorb_x_term () =
let (a12, a1, a2, a, b12, b1, b2, b) = !ng in
match (!xsource) () with
| Some term ->
let new_ng = Z.((a2 + (a12 * term),
a + (a1 * term), a12, a1,
b2 + (b12 * term),
b + (b1 * term), b12, b1)) in
if any_too_big new_ng then
(* Pretend all further x terms are infinite. *)
begin
ng := (a12, a1, a12, a1, b12, b1, b12, b1);
xsource := no_terms_source
end
else
ng := new_ng
| None -> ng := (a12, a1, a12, a1, b12, b1, b12, b1)
in
let absorb_y_term () =
let (a12, a1, a2, a, b12, b1, b2, b) = !ng in
match (!ysource) () with
| Some term ->
let new_ng = Z.((a1 + (a12 * term), a12,
a + (a2 * term), a2,
b1 + (b12 * term), b12,
b + (b2 * term), b2)) in
if any_too_big new_ng then
(* Pretend all further y terms are infinite. *)
begin
ng := (a12, a12, a2, a2, b12, b12, b2, b2);
ysource := no_terms_source
end
else
ng := new_ng
| None -> ng := (a12, a12, a2, a2, b12, b12, b2, b2)
in
let rec loop () =
if all_b_are_zero () then
None (* There are no more terms to output. *)
else
let (_, _, _, _, b12, b1, b2, b) = !ng in
if equal_zero b && equal_zero b2 then
(absorb_x_term (); loop ())
else if equal_zero b || equal_zero b2 then
(absorb_y_term (); loop ())
else if equal_zero b1 then
(absorb_x_term (); loop ())
else
let (a12, a1, a2, a, _, _, _, _) = !ng in
let (q12, r12) = divide a12 b12
and (q1, r1) = divide a1 b1
and (q2, r2) = divide a2 b2
and (q, r) = divide a b in
if Z.sign b12 <> 0 && all_four_equal (q12, q1, q2, q) then
begin
ng := (b12, b1, b2, b, r12, r1, r2, r);
(* Return a term--or, if a magnitude threshold is
reached, return no more terms . *)
infinitize q
end
else
begin
(* Put a1, a2, and a over a common denominator and
compare some magnitudes. *)
let n1 = Z.(a1 * b2 * b)
and n2 = Z.(a2 * b1 * b)
and n = Z.(a * b1 * b2) in
if Z.(abs (n1 - n) > abs (n2 - n)) then
(absorb_x_term (); loop ())
else
(absorb_y_term (); loop ())
end
in
make (fun () -> loop ())
end
let ( + ) = apply_ng8 (ng8_of_ints (0, 1, 1, 0, 0, 0, 0, 1))
let ( - ) = apply_ng8 (ng8_of_ints (0, 1, (-1), 0, 0, 0, 0, 1))
let ( * ) = apply_ng8 (ng8_of_ints (1, 0, 0, 0, 0, 0, 0, 1))
let ( / ) = apply_ng8 (ng8_of_ints (0, 1, 0, 0, 0, 0, 1, 0))
let ( ~- ) =
let neg = apply_ng8 (ng8_of_ints (0, 0, (-1), 0, 0, 0, 0, 1)) in
fun x -> neg x x
let pow =
let one = of_int 1 in
let reciprocal =
apply_ng8 (ng8_of_ints (0, 0, 0, 1, 0, 1, 0, 0)) in
let reciprocal x = reciprocal x x in
fun cf i -> let rec loop (x : t) (n : int) (accum : t) =
if Int.(1 < n) then
let nhalf = Int.(div n 2)
and xsquare = x * x in
if Int.(add nhalf nhalf <> n) then
loop xsquare nhalf (accum * x)
else
loop xsquare nhalf accum
else if Int.(n = 1) then
(accum * x)
else
accum
in
if 0 <= i then
loop cf i one
else
reciprocal (loop cf Int.(neg i) one)
let zero = of_int 0
let one = of_int 1
let two = of_int 2
let three = of_int 3
let four = of_int 4
let one_fourth = of_ints 1 4
let one_third = of_ints 1 3
let one_half = of_ints 1 2
let two_thirds = of_ints 2 3
let three_fourths = of_ints 3 4
let golden_ratio = constant_term_of_int 1
let silver_ratio = constant_term_of_int 2
let sqrt2 = silver_ratio - one
let sqrt5 = (two * golden_ratio) - one
end
module CF = Continued_fraction
;;
let i2cf = CF.of_int
and r2cf = CF.of_ints
and constant_term_cf = CF.constant_term_of_int
and cf2string = CF.to_string
;;
let make_pad n = String.make n ' '
;;
let show (expression, cf, note) =
let cf_string = cf2string cf in
let expr_sz = String.length expression
and cf_sz = String.length cf_string
and note_sz = String.length note in
let expr_pad_sz = max (19 - expr_sz) 0
and cf_pad_sz = if note_sz = 0 then 0 else max (48 - cf_sz) 0 in
let expr_pad = make_pad expr_pad_sz
and cf_pad = make_pad cf_pad_sz in
print_string expr_pad;
print_string expression;
print_string " => ";
print_string cf_string;
print_string cf_pad;
print_string note;
print_endline ""
;;
show ("golden ratio", CF.golden_ratio, "(1 + sqrt(5))/2");
show ("silver ratio", CF.silver_ratio, "(1 + sqrt(2))");
show ("sqrt2", CF.sqrt2, "");
show ("sqrt5", CF.sqrt5, "");
show ("1/4", CF.one_fourth, "");
show ("1/3", CF.one_third, "");
show ("1/2", CF.one_half, "");
show ("2/3", CF.two_thirds, "");
show ("3/4", CF.three_fourths, "");
show ("13/11", r2cf 13 11, "");
show ("22/7", r2cf 22 7, "approximately pi");
show ("0", CF.zero, "");
show ("1", CF.one, "");
show ("2", CF.two, "");
show ("3", CF.three, "");
show ("4", CF.four, "");
show ("4 + 3", CF.(four + three), "");
show ("4 - 3", CF.(four - three), "");
show ("4 * 3", CF.(four * three), "");
show ("4 / 3", CF.(four / three), "");
show ("4 ** 3", CF.(pow four 3), "");
show ("4 ** (-3)", CF.(pow four (-3)), "");
show ("negative 4", CF.(-four), "");
CF.(show ("(1 + 1/sqrt(2))/2",
(one + (one / sqrt2)) / two, "method 1"));
CF.(show ("(1 + 1/sqrt(2))/2",
silver_ratio * pow sqrt2 (-3), "method 2"));
CF.(show ("(1 + 1/sqrt(2))/2",
(pow silver_ratio 2 + one) / (four * two), "method 3"));
show ("sqrt2 + sqrt2", CF.(sqrt2 + sqrt2), "");
show ("sqrt2 - sqrt2", CF.(sqrt2 - sqrt2), "");
show ("sqrt2 * sqrt2", CF.(sqrt2 * sqrt2), "");
show ("sqrt2 / sqrt2", CF.(sqrt2 / sqrt2), "");
</syntaxhighlight>
{{out}}
<pre>$ ocamlfind ocamlc -linkpkg -package zarith bivariate_continued_fraction_task.ml && ./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))
sqrt2 => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt5 => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...]
1/4 => [0;4]
1/3 => [0;3]
1/2 => [0;2]
2/3 => [0;1,2]
3/4 => [0;1,3]
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
0 => [0]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
4 + 3 => [7]
4 - 3 => [1]
4 * 3 => [12]
4 / 3 => [1;3]
4 ** 3 => [64]
4 ** (-3) => [0;64]
negative 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
sqrt2 + sqrt2 => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt2 - sqrt2 => [0]
sqrt2 * sqrt2 => [2]
sqrt2 / sqrt2 => [1]
</pre>
=={{header|Owl Lisp}}==
{{trans|Mercury}}
{{works with|Owl Lisp|0.2.1}}
<syntaxhighlight lang="scheme">
;;--------------------------------------------------------------------
;;
;; A continued fraction will be represented as an infinite lazy list
;; of terms, where a term is either an integer or the symbol 'infinity
;;
;;--------------------------------------------------------------------
(define (cf2maxterms-string maxterms cf)
(let repeat ((p cf)
(i 0)
(s "["))
(let ((term (lcar p))
(rest (lcdr p)))
(if (eq? term 'infinity)
(string-append s "]")
(if (>= i maxterms)
(string-append s ",...]")
(let ((separator (case i
((0) "")
((1) ";")
(else ",")))
(term-str (number->string term)))
(repeat rest (+ i 1)
(string-append s separator term-str))))))))
(define (cf2string cf)
(cf2maxterms-string 20 cf))
;;--------------------------------------------------------------------
(define (repeated-term-cf term)
(lunfold (lambda (x) (values x x))
term
(lambda (x) #false)))
(define cf-end (repeated-term-cf 'infinity))
(define (i2cf term) (pair term cf-end))
(define (r2cf fraction)
;; The entire finite-length continued fraction is constructed in
;; reverse order. The list is then rebuilt in the correct order, and
;; given an infinite number of 'infinity terms as its tail.
(let repeat ((n (numerator fraction))
(d (denominator fraction))
(revlst #null))
(if (zero? d)
(lappend (reverse revlst) cf-end)
(let-values (((q r) (truncate/ n d)))
(repeat d r (cons q revlst))))))
(define (->cf x)
(cond ((integer? x) (i2cf x))
((rational? x) (r2cf x))
(else x)))
;;--------------------------------------------------------------------
(define (maybe-divide a b)
(if (zero? b)
(values #false #false)
(truncate/ a b)))
(define integer-exceeds-processing-threshold?
(let ((threshold+1 (expt 2 512)))
(lambda (i) (>= (abs i) threshold+1))))
(define (any-exceed-processing-threshold? lst)
(any integer-exceeds-processing-threshold? lst))
(define integer-exceeds-infinitizing-threshold?
(let ((threshold+1 (expt 2 64)))
(lambda (i) (>= (abs i) threshold+1))))
(define (ng8-values ng)
(values (list-ref ng 0)
(list-ref ng 1)
(list-ref ng 2)
(list-ref ng 3)
(list-ref ng 4)
(list-ref ng 5)
(list-ref ng 6)
(list-ref ng 7)))
(define (absorb-x-term ng x y)
(let-values (((a12 a1 a2 a b12 b1 b2 b) (ng8-values ng)))
(define (no-more-x)
(values (list a12 a1 a12 a1 b12 b1 b12 b1) cf-end y))
(let ((term (lcar x)))
(if (eq? term 'infinity)
(no-more-x)
(let ((new-ng (list (+ a2 (* a12 term))
(+ a (* a1 term)) a12 a1
(+ b2 (* b12 term))
(+ b (* b1 term)) b12 b1)))
(cond ((any-exceed-processing-threshold? new-ng)
;; Pretend we have reached the end of x.
(no-more-x))
(else (values new-ng (lcdr x) y))))))))
(define (absorb-y-term ng x y)
(let-values (((a12 a1 a2 a b12 b1 b2 b) (ng8-values ng)))
(define (no-more-y)
(values (list a12 a12 a2 a2 b12 b12 b2 b2) x cf-end))
(let ((term (lcar y)))
(if (eq? term 'infinity)
(no-more-y)
(let ((new-ng (list (+ a1 (* a12 term)) a12
(+ a (* a2 term)) a2
(+ b1 (* b12 term)) b12
(+ b (* b2 term)) b2)))
(cond ((any-exceed-processing-threshold? new-ng)
;; Pretend we have reached the end of y.
(no-more-y))
(else (values new-ng x (lcdr y)))))))))
(define (apply-ng8 ng x y)
(let repeat ((ng ng)
(x (->cf x))
(y (->cf y)))
(let-values (((a12 a1 a2 a b12 b1 b2 b) (ng8-values ng)))
(cond ((every zero? (list b12 b1 b2 b)) cf-end)
((every zero? (list b2 b))
(let-values (((ng x y) (absorb-x-term ng x y)))
(repeat ng x y)))
((any zero? (list b2 b))
(let-values (((ng x y) (absorb-y-term ng x y)))
(repeat ng x y)))
((zero? b1)
(let-values (((ng x y) (absorb-x-term ng x y)))
(repeat ng x y)))
(else
(let-values (((q12 r12) (maybe-divide a12 b12))
((q1 r1) (maybe-divide a1 b1))
((q2 r2) (maybe-divide a2 b2))
((q r) (maybe-divide a b)))
(cond
((and (not (zero? b12))
(= q q12) (= q q1) (= q q2))
;; Output a term.
(if (integer-exceeds-infinitizing-threshold? q)
cf-end ; Pretend the term is infinite.
(let ((new-ng (list b12 b1 b2 b r12 r1 r2 r)))
(pair q (repeat new-ng x y)))))
(else
;; Put a1, a2, and a over a common denominator
;; and compare some magnitudes.
(let ((n1 (* a1 b2 b))
(n2 (* a2 b1 b))
(n (* a b1 b2)))
(let ((absorb-term
(if (> (abs (- n1 n)) (abs (- n2 n)))
absorb-x-term
absorb-y-term)))
(let-values (((ng x y) (absorb-term ng x y)))
(repeat ng x y))))))))))))
;;--------------------------------------------------------------------
(define (make-ng8-operation ng) (lambda (x y) (apply-ng8 ng x y)))
(define cf+ (make-ng8-operation '(0 1 1 0 0 0 0 1)))
(define cf- (make-ng8-operation '(0 1 -1 0 0 0 0 1)))
(define cf* (make-ng8-operation '(1 0 0 0 0 0 0 1)))
(define cf/ (make-ng8-operation '(0 1 0 0 0 0 1 0)))
;;--------------------------------------------------------------------
(define golden-ratio (repeated-term-cf 1))
(define silver-ratio (repeated-term-cf 2))
(define sqrt2 (pair 1 silver-ratio))
(define sqrt5 (pair 2 (repeated-term-cf 4)))
;;--------------------------------------------------------------------
(define (show expression cf note)
(let* ((cf (cf2string cf))
(expr-len (string-length expression))
(expr-pad-len (max 0 (- 18 expr-len)))
(expr-pad (make-string expr-pad-len #\space)))
(display expr-pad)
(display expression)
(display " => ")
(display cf)
(unless (string=? note "")
(let* ((cf-len (string-length cf))
(cf-pad-len (max 0 (- 48 cf-len)))
(cf-pad (make-string cf-pad-len #\space)))
(display cf-pad)
(display note)))
(newline)))
(show "13/11 + 1/2" (cf+ 13/11 1/2) "(cf+ 13/11 1/2)")
(show "22/7 + 1/2" (cf+ 22/7 1/2) "(cf+ 22/7 1/2)")
(show "22/7 * 1/2" (cf* 22/7 1/2) "(cf* 22/7 1/2)")
(show "golden ratio" golden-ratio "(repeated-term-cf 1)")
(show "silver ratio" silver-ratio "(repeated-term-cf 2)")
(show "sqrt(2)" sqrt2 "(pair 1 silver-ratio)")
(show "sqrt(2)" (cf- silver-ratio 1) "(cf- silver-ratio 1)")
(show "sqrt(5)" sqrt5 "(pair 2 (repeated-term-cf 4)")
(show "sqrt(5)" (cf- (cf* 2 golden-ratio) 1)
"(cf- (cf* 2 golden-ratio) 1)")
(show "sqrt(2) + sqrt(2)" (cf+ sqrt2 sqrt2) "(cf+ sqrt2 sqrt2)")
(show "sqrt(2) - sqrt(2)" (cf- sqrt2 sqrt2) "(cf- sqrt2 sqrt2)")
(show "sqrt(2) * sqrt(2)" (cf* sqrt2 sqrt2) "(cf* sqrt2 sqrt2)")
(show "sqrt(2) / sqrt(2)" (cf/ sqrt2 sqrt2) "(cf/ sqrt2 sqrt2)")
(show "(1 + 1/sqrt(2))/2" (cf/ (cf+ 1 (cf/ 1 sqrt2)) 2)
"(cf/ (cf+ 1 (cf/ 1 sqrt2)) 2)")
(show "(1 + 1/sqrt(2))/2" (apply-ng8 '(0 1 0 0 0 0 2 0)
silver-ratio sqrt2)
"(apply-ng8 '(0 1 0 0 0 0 2 0) sqrt2 sqrt2)")
(show "(1 + 1/sqrt(2))/2" (apply-ng8 '(1 0 0 1 0 0 0 8)
silver-ratio silver-ratio)
"(apply-ng8 '(1 0 0 1 0 0 0 8) silver-ratio silver-ratio)")
;;--------------------------------------------------------------------
</syntaxhighlight>
{{out}}
<pre>$ ol continued-fraction-task-Owl.scm
13/11 + 1/2 => [1;1,2,7] (cf+ 13/11 1/2)
22/7 + 1/2 => [3;1,1,1,4] (cf+ 22/7 1/2)
22/7 * 1/2 => [1;1,1,3] (cf* 22/7 1/2)
golden ratio => [1;1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...] (repeated-term-cf 1)
silver ratio => [2;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] (repeated-term-cf 2)
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] (pair 1 silver-ratio)
sqrt(2) => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...] (cf- silver-ratio 1)
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] (pair 2 (repeated-term-cf 4)
sqrt(5) => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...] (cf- (cf* 2 golden-ratio) 1)
sqrt(2) + sqrt(2) => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...] (cf+ sqrt2 sqrt2)
sqrt(2) - sqrt(2) => [0] (cf- sqrt2 sqrt2)
sqrt(2) * sqrt(2) => [2] (cf* sqrt2 sqrt2)
sqrt(2) / sqrt(2) => [1] (cf/ sqrt2 sqrt2)
(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,...] (cf/ (cf+ 1 (cf/ 1 sqrt2)) 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,...] (apply-ng8 '(0 1 0 0 0 0 2 0) sqrt2 sqrt2)
(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,...] (apply-ng8 '(1 0 0 1 0 0 0 8) silver-ratio silver-ratio)</pre>
=={{header|Phix}}==
Line 3,777 ⟶ 9,509:
ix = 1
iy = 2
yoverflow = 4
def too_big (values):
Line 3,791 ⟶ 9,524:
def absorb_x_term ():
(a12, a1, a2, a, b12, b1, b2, b) = env[ng]
if env[
term = None
else:
Line 3,803 ⟶ 9,536:
else:
env[ng] = (a12, a1, a12, a1, b12, b1, b12, b1)
env[
else:
env[ng] = (a12, a1, a12, a1, b12, b1, b12, b1)
Line 3,810 ⟶ 9,543:
def absorb_y_term ():
(a12, a1, a2, a, b12, b1, b2, b) = env[ng]
if env[
term = None
else:
Line 3,822 ⟶ 9,555:
else:
env[ng] = (a12, a12, a2, a2, b12, b12, b2, b2)
env[
else:
env[ng] = (a12, a12, a2, a2, b12, b12, b2, b2)
Line 3,868 ⟶ 9,601:
return retval
return ContinuedFraction (func, [ng8_tuple, 0, 0, False, False])
#---------------------------------------------------------------------
Line 4,378 ⟶ 10,111:
sqrt(2) * sqrt(2) => [2]
sqrt(2) / sqrt(2) => [1]
</pre>
=={{header|Standard ML}}==
{{trans|Scheme}}
{{works with|Poly/ML|5.9}}
{{works with|MLton|20210117}}
There are a few extras in here.
<syntaxhighlight lang="sml">
(*------------------------------------------------------------------*)
signature CONTINUED_FRACTION =
sig
(* A termGenerator thunk generates terms, which a continuedFraction
structure memoizes. *)
type termGenerator = unit -> IntInf.int option
type continuedFraction
(* Create a continued fraction. *)
val make : termGenerator -> continuedFraction
(* Does the indexed term exist? *)
val exists : continuedFraction * int -> bool
(* Retrieving the indexed term. *)
val sub : continuedFraction * int -> IntInf.int
(* Using a continuedFraction as a termGenerator thunk. *)
val toTermGenerator : continuedFraction -> termGenerator
(* Producing a human-readable string. *)
val toMaxTermsString : int -> continuedFraction -> string
val defaultMaxTerms : int ref
val toString : continuedFraction -> string
(* - - - - - - - - - - - - - - - - - - - *)
(* Creating some specific kinds of continued fractions: *)
(* Representing an integer. *)
val fromIntInf : IntInf.int -> continuedFraction
val fromInt : int -> continuedFraction
val i2cf : int -> continuedFraction (* Synonym for fromInt. *)
(* With one term repeated forever. *)
val withConstantIntInfTerm : IntInf.int -> continuedFraction
val withConstantIntTerm : int -> continuedFraction
(* Representing a rational number. *)
val fromIntInfNumerDenom : IntInf.int * IntInf.int ->
continuedFraction
val fromIntNumerDenom : int * int -> continuedFraction
val r2cf : int * int ->
continuedFraction (* Synonym for fromIntNumerDenom. *)
(* Representing arithmetic. (I have not bothered here to implement
ng4, although one likely would wish to have ng4 as well.) *)
type ng8 = IntInf.int * IntInf.int * IntInf.int * IntInf.int *
IntInf.int * IntInf.int * IntInf.int * IntInf.int
val ng8_fromInt : int * int * int * int *
int * int * int * int -> ng8
val ng8StoppingProcessingThreshold : IntInf.int ref
val ng8InfinitizationThreshold : IntInf.int ref
val apply_ng8 : ng8 ->
continuedFraction * continuedFraction ->
continuedFraction
val + : continuedFraction * continuedFraction -> continuedFraction
val - : continuedFraction * continuedFraction -> continuedFraction
val * : continuedFraction * continuedFraction -> continuedFraction
val / : continuedFraction * continuedFraction -> continuedFraction
val ~ : continuedFraction -> continuedFraction
val pow : continuedFraction * int -> continuedFraction
(* - - - - - - - - - - - - - - - - - - - *)
(* Miscellanous continued fractions: *)
val zero : continuedFraction
val one : continuedFraction
val two : continuedFraction
val three : continuedFraction
val four : continuedFraction
val one_fourth : continuedFraction
val one_third : continuedFraction
val one_half : continuedFraction
val two_thirds : continuedFraction
val three_fourths : continuedFraction
val goldenRatio : continuedFraction
val silverRatio : continuedFraction
val sqrt2 : continuedFraction
val sqrt5 : continuedFraction
end
structure ContinuedFraction : CONTINUED_FRACTION =
struct
type termGenerator = unit -> IntInf.int option
type cfRecord = {
terminated : bool, (* Is the generator exhausted? *)
memoCount : int, (* How many terms are memoized? *)
memo : IntInf.int array, (* Memoized terms. *)
generate : termGenerator (* The source of terms. *)
}
type continuedFraction = cfRecord ref
fun make generator =
ref {
terminated = false,
memoCount = 0,
memo = Array.array (32, IntInf.fromInt 0),
generate = generator
}
fun resizeIfNecessary (cf : continuedFraction, i) =
let
val record = !cf
in
if #terminated record then
()
else if i < Array.length (#memo record) then
()
else
let
val newSize = 2 * (i + 1)
val newMemo = Array.array (newSize, IntInf.fromInt 0)
fun copyTerms i =
if i = #memoCount record then
()
else
let
val term = Array.sub (#memo record, i)
in
Array.update (newMemo, i, term);
copyTerms (i + 1)
end
val newRecord : cfRecord = {
terminated = false,
memoCount = #memoCount record,
memo = newMemo,
generate = #generate record
}
in
copyTerms 0;
cf := newRecord
end
end
fun updateTerms (cf : continuedFraction, i) =
let
val record = !cf
in
if #terminated record then
()
else if i < #memoCount record then
()
else
case (#generate record) () of
Option.NONE =>
let
val newRecord : cfRecord = {
terminated = true,
memoCount = #memoCount record,
memo = #memo record,
generate = #generate record
}
in
cf := newRecord
end
| Option.SOME term =>
let
val () = Array.update (#memo record, #memoCount record,
term)
val newRecord : cfRecord = {
terminated = false,
memoCount = (#memoCount record) + 1,
memo = #memo record,
generate = #generate record
}
in
cf := newRecord;
updateTerms (cf, i)
end
end
fun exists (cf : continuedFraction, i) =
(resizeIfNecessary (cf, i);
updateTerms (cf, i);
i < #memoCount (!cf))
fun sub (cf : continuedFraction, i) =
let
val record = !cf
in
if #memoCount record <= i then
raise Domain
else
Array.sub (#memo record, i)
end
fun toTermGenerator (cf : continuedFraction) =
let
val i : int ref = ref 0
in
fn () =>
let
val j = !i
in
if exists (cf, j) then
(i := j + 1;
Option.SOME (sub (cf, j)))
else
Option.NONE
end
end
fun toMaxTermsString maxTerms =
if maxTerms < 1 then
raise Domain
else
fn (cf : continuedFraction) =>
let
fun loop (i, accum) =
if not (exists (cf, i)) then
accum ^ "]"
else if i = maxTerms then
accum ^ ",...]"
else
let
val separator =
if i = 0 then
""
else if i = 1 then
";"
else
","
val termString = IntInf.toString (sub (cf, i))
in
loop (i + 1, accum ^ separator ^ termString)
end
in
loop (0, "[")
end
val defaultMaxTerms : int ref = ref 20
val toString = toMaxTermsString (!defaultMaxTerms)
fun fromIntInf i =
let
val done : bool ref = ref false
in
make (fn () =>
if !done then
Option.NONE
else
(done := true;
Option.SOME i))
end
fun fromInt i =
fromIntInf (IntInf.fromInt i)
val i2cf = fromInt
fun withConstantIntInfTerm i =
make (fn () => Option.SOME i)
fun withConstantIntTerm i =
withConstantIntInfTerm (IntInf.fromInt i)
fun fromIntInfNumerDenom (n, d) =
let
val zero = IntInf.fromInt 0
val state = ref (n, d)
in
make (fn () =>
let
val (n, d) = !state
in
if d = zero then
Option.NONE
else
let
val (q, r) = IntInf.divMod (n, d)
in
state := (d, r);
Option.SOME q
end
end)
end
fun fromIntNumerDenom (n, d) =
fromIntInfNumerDenom (IntInf.fromInt n, IntInf.fromInt d)
val r2cf = fromIntNumerDenom
type ng8 = IntInf.int * IntInf.int * IntInf.int * IntInf.int *
IntInf.int * IntInf.int * IntInf.int * IntInf.int
fun ng8_fromInt (a12, a1, a2, a, b12, b1, b2, b) =
let
val f = IntInf.fromInt
in
(f a12, f a1, f a2, f a, f b12, f b1, f b2, f b)
end
val ng8StoppingProcessingThreshold =
ref (IntInf.pow (IntInf.fromInt 2, 512))
val ng8InfinitizationThreshold =
ref (IntInf.pow (IntInf.fromInt 2, 64))
fun tooBig term =
abs (term) >= abs (!ng8StoppingProcessingThreshold)
fun anyTooBig (a, b, c, d, e, f, g, h) =
tooBig (a) orelse tooBig (b) orelse
tooBig (c) orelse tooBig (d) orelse
tooBig (e) orelse tooBig (f) orelse
tooBig (g) orelse tooBig (h)
fun infinitize term =
if abs (term) >= abs (!ng8InfinitizationThreshold) then
Option.NONE
else
Option.SOME term
fun noTermsSource () =
Option.NONE
val equalZero =
let
val zero = IntInf.fromInt 0
in
fn b => (b = zero)
end
val divide =
let
val zero = IntInf.fromInt 0
in
fn (a, b) =>
if equalZero b then
(zero, zero)
else
IntInf.divMod (a, b)
end
fun apply_ng8 ng =
fn (x, y) =>
let
val ng = ref ng
and xsource = ref (toTermGenerator x)
and ysource = ref (toTermGenerator y)
fun all_b_areZero () =
let
val (_, _, _, _, b12, b1, b2, b) = !ng
in
equalZero b andalso
equalZero b2 andalso
equalZero b1 andalso
equalZero b12
end
fun allFourEqual (a, b, c, d) =
a = b andalso a = c andalso a = d
fun absorb_x_term () =
let
val (a12, a1, a2, a, b12, b1, b2, b) = !ng
in
case (!xsource) () of
Option.SOME term =>
let
val new_ng = (a2 + (a12 * term),
a + (a1 * term), a12, a1,
b2 + (b12 * term),
b + (b1 * term), b12, b1)
in
if anyTooBig new_ng then
(* Pretend all further x terms are infinite. *)
(ng := (a12, a1, a12, a1, b12, b1, b12, b1);
xsource := noTermsSource)
else
ng := new_ng
end
| Option.NONE =>
ng := (a12, a1, a12, a1, b12, b1, b12, b1)
end
fun absorb_y_term () =
let
val (a12, a1, a2, a, b12, b1, b2, b) = !ng
in
case (!ysource) () of
Option.SOME term =>
let
val new_ng = (a1 + (a12 * term), a12,
a + (a2 * term), a2,
b1 + (b12 * term), b12,
b + (b2 * term), b2)
in
if anyTooBig new_ng then
(* Pretend all further y terms are infinite. *)
(ng := (a12, a12, a2, a2, b12, b12, b2, b2);
ysource := noTermsSource)
else
ng := new_ng
end
| Option.NONE =>
ng := (a12, a12, a2, a2, b12, b12, b2, b2)
end
fun loop () =
(* Although my Scheme version of this program used mutual
tail recursion, here there is only single tail
recursion. The difference is that in SML we cannot
rely on properness of mutual tail calls, the way we
can in standard Scheme. *)
if all_b_areZero () then
Option.NONE (* There are no more terms to output. *)
else
let
val (_, _, _, _, b12, b1, b2, b) = !ng
in
if equalZero b andalso equalZero b2 then
(absorb_x_term (); loop ())
else if equalZero b orelse equalZero b2 then
(absorb_y_term (); loop ())
else if equalZero b1 then
(absorb_x_term (); loop ())
else
let
val (a12, a1, a2, a, _, _, _, _) = !ng
val (q12, r12) = divide (a12, b12)
and (q1, r1) = divide (a1, b1)
and (q2, r2) = divide (a2, b2)
and (q, r) = divide (a, b)
in
if b12 <> 0 andalso
allFourEqual (q12, q1, q2, q) then
(ng := (b12, b1, b2, b, r12, r1, r2, r);
(* Return a term--or, if a magnitude threshold
is reached, return no more terms . *)
infinitize q)
else
let
(* Put a1, a2, and a over a common
denominator and compare some
magnitudes. *)
val n1 = a1 * b2 * b
and n2 = a2 * b1 * b
and n = a * b1 * b2
in
if abs (n1 - n) > abs (n2 - n) then
(absorb_x_term (); loop ())
else
(absorb_y_term (); loop ())
end
end
end
in
make (fn () => loop ())
end
val op+ = apply_ng8 (ng8_fromInt (0, 1, 1, 0, 0, 0, 0, 1))
val op- = apply_ng8 (ng8_fromInt (0, 1, ~1, 0, 0, 0, 0, 1))
val op* = apply_ng8 (ng8_fromInt (1, 0, 0, 0, 0, 0, 0, 1))
val op/ = apply_ng8 (ng8_fromInt (0, 1, 0, 0, 0, 0, 1, 0))
val op~ =
let
val neg = apply_ng8 (ng8_fromInt (0, 0, ~1, 0, 0, 0, 0, 1))
in
fn (x) => neg (x, x)
end
val pow =
let
val one = i2cf 1
val reciprocal =
fn (x) =>
apply_ng8 (ng8_fromInt (0, 0, 0, 1, 0, 1, 0, 0))
(x, x)
in
fn (cf, i) =>
let
fun loop (x, n, accum) =
if 1 < n then
let
val nhalf = n div 2
and xsquare = x * x
in
if Int.+ (nhalf, nhalf) <> n then
loop (xsquare, nhalf, accum * x)
else
loop (xsquare, nhalf, accum)
end
else if n = 1 then
accum * x
else
accum
in
if 0 <= i then
loop (cf, i, one)
else
reciprocal (loop (cf, Int.~ i, one))
end
end
val zero = i2cf 0
val one = i2cf 1
val two = i2cf 2
val three = i2cf 3
val four = i2cf 4
val one_fourth = r2cf (1, 4)
val one_third = r2cf (1, 3)
val one_half = r2cf (1, 2)
val two_thirds = r2cf (2, 3)
val three_fourths = r2cf (3, 4)
val goldenRatio = withConstantIntTerm 1
val silverRatio = withConstantIntTerm 2
val sqrt2 = silverRatio - one
val sqrt5 = (two * goldenRatio) - one
end
(*------------------------------------------------------------------*)
fun makePad n =
String.implode (List.tabulate (n, fn i => #" "))
fun show (expression : string,
cf : ContinuedFraction.continuedFraction,
note : string) =
let
val cfString = ContinuedFraction.toString cf
val exprSz = size expression
and cfSz = size cfString
and noteSz = size note
val exprPadSize = Int.max (19 - exprSz, 0)
and cfPadSize = if noteSz = 0 then 0 else Int.max (48 - cfSz, 0)
val exprPad = makePad exprPadSize
and cfPad = makePad cfPadSize
in
print exprPad;
print expression;
print " => ";
print cfString;
print cfPad;
print note;
print "\n"
end;
let
open ContinuedFraction
in
show ("golden ratio", goldenRatio, "(1 + sqrt(5))/2");
show ("silver ratio", silverRatio, "(1 + sqrt(2))");
show ("sqrt2", sqrt2, "");
show ("sqrt5", sqrt5, "");
show ("1/4", one_fourth, "");
show ("1/3", one_third, "");
show ("1/2", one_half, "");
show ("2/3", two_thirds, "");
show ("3/4", three_fourths, "");
show ("13/11", r2cf (13, 11), "");
show ("22/7", r2cf (22, 7), "approximately pi");
show ("0", zero, "");
show ("1", one, "");
show ("2", two, "");
show ("3", three, "");
show ("4", four, "");
show ("4 + 3", four + three, "");
show ("4 - 3", four - three, "");
show ("4 * 3", four * three, "");
show ("4 / 3", four / three, "");
show ("4 ** 3", pow (four, 3), "");
show ("4 ** (-3)", pow (four, ~3), "");
show ("negative 4", ~four, "");
show ("(1 + 1/sqrt(2))/2",
(one + (one / sqrt2)) / two, "method 1");
show ("(1 + 1/sqrt(2))/2",
silverRatio * pow (sqrt2, ~3), "method 2");
show ("(1 + 1/sqrt(2))/2",
(pow (silverRatio, 2) + one) / (four * two), "method 3");
show ("sqrt2 + sqrt2", sqrt2 + sqrt2, "");
show ("sqrt2 - sqrt2", sqrt2 - sqrt2, "");
show ("sqrt2 * sqrt2", sqrt2 * sqrt2, "");
show ("sqrt2 / sqrt2", sqrt2 / sqrt2, "");
()
end;
(*------------------------------------------------------------------*)
(* local variables: *)
(* mode: sml *)
(* sml-indent-level: 2 *)
(* sml-indent-args: 2 *)
(* end: *)
</syntaxhighlight>
{{out}}
<pre>$ poly --script bivariate_continued_fraction_task.sml
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))
sqrt2 => [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]
sqrt5 => [2;4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,...]
1/4 => [0;4]
1/3 => [0;3]
1/2 => [0;2]
2/3 => [0;1,2]
3/4 => [0;1,3]
13/11 => [1;5,2]
22/7 => [3;7] approximately pi
0 => [0]
1 => [1]
2 => [2]
3 => [3]
4 => [4]
4 + 3 => [7]
4 - 3 => [1]
4 * 3 => [12]
4 / 3 => [1;3]
4 ** 3 => [64]
4 ** (-3) => [0;64]
negative 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
sqrt2 + sqrt2 => [2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,...]
sqrt2 - sqrt2 => [0]
sqrt2 * sqrt2 => [2]
sqrt2 / sqrt2 => [1]
</pre>
Line 4,512 ⟶ 10,885:
=={{header|Wren}}==
{{trans|Kotlin}}
<syntaxhighlight lang="
construct new() {
_cfn = 0
| |||