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) |
|||
| (16 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 1,823:
===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].
Line 2,296 ⟶ 2,298:
@(exrat_make (0, 1), exrat_make (0, 1))
else
(* Do integer division
particular function does floor division if the divisor is
positive, ceiling division if the divisor is negative. Thus the
Line 2,623 ⟶ 2,625:
{{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.
Line 6,124 ⟶ 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>
Line 6,476 ⟶ 6,752:
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 6,922 ⟶ 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 7,537 ⟶ 8,486:
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 9,187 ⟶ 10,885:
=={{header|Wren}}==
{{trans|Kotlin}}
<syntaxhighlight lang="
construct new() {
_cfn = 0
| |||