From 1733f1da1c479e6392e73c180933b816aa2104de Mon Sep 17 00:00:00 2001
From: Glenn Rice <grice1@missouriwestern.edu>
Date: Wed, 29 Mar 2023 17:29:14 -0500
Subject: [PATCH] Improved non-recursive version of contFrac in
 contextFraction.pl.

This version uses a loop instead of a recursive function.  Benchmarking
shows it is at least twice as fast.  It also should use less memory
overhead as it does not store the steps in constructing the continued
fraction in arrays.  It only keeps what is needed to continue.  That is
the last two values of the convergent numerator and denominator.

A basic unit test is added to check that this does the right thing.
---
 macros/contexts/contextFraction.pl | 73 ++++++++++--------------------
 t/contexts/fraction.t              | 23 ++++++++--
 2 files changed, 43 insertions(+), 53 deletions(-)

diff --git a/macros/contexts/contextFraction.pl b/macros/contexts/contextFraction.pl
index 3cb7850e1b..a7e37fa62a 100644
--- a/macros/contexts/contextFraction.pl
+++ b/macros/contexts/contextFraction.pl
@@ -296,59 +296,34 @@ sub Init {
 	main::PG_restricted_eval('sub Fraction {Value->Package("Fraction()")->new(@_)};');
 }
 
-#
 # contFrac($x, $maxdenominator)
-#
-# Recursive subroutine that takes positive real input $x and outputs
-# an array (a,b) where a/b is a very good fraction approximation with
-# b no larger than maxdenominator
-#
-
+# Subroutine that takes a positive real input $x and outputs an array
+# (a,b) where a/b is a very good fraction approximation with b no
+# larger than maxdenominator.
 sub contFrac {
-	my $x              = shift;
-	my $maxdenominator = shift;
-	my %sequences      = @_;      # an => continued fraction sequence (reference)
-								  # hn => sequence of numerators (reference)
-								  # kn => sequence of denominators (reference)
-
-	# dereference sequences
-	my @an = (int($x));
-	@an = @{ $sequences{"an"} } if defined($sequences{"an"});
-	my @hn = (int($x));
-	@hn = @{ $sequences{"hn"} } if defined($sequences{"hn"});
-	my @kn = (1);
-	@kn = @{ $sequences{"kn"} } if defined($sequences{"kn"});
-
-	# calculate what real the continued fraciton process leaves at this level
+	my ($x, $maxdenominator) = @_;
+
 	my $step = $x;
-	for my $i (0 .. $#an - 1) { $step = ($step - $an[$i])**(-1); }
-	# if this is an integer, stop
-	if ($step == int($step)) { return ($hn[-1], $kn[-1]); }
-
-	$step = ($step - $an[-1])**(-1);
-
-	# next integer from continued fraction sequence
-	# next numerator and denominator, according to continued fraction formulas
-	my $newa = int($step);
-	my $newh;
-	my $newk;
-	if   ($#an > 0) { $newh = $newa * $hn[-1] + $hn[-2]; }
-	else            { $newh = $newa * $an[0] + 1; }
-	if   ($#an > 0) { $newk = $newa * $kn[-1] + $kn[-2]; }
-	else            { $newk = $newa; }
-
-	# machine rounding error may begin to make denominators skyrocket out of control
-	if ($newk > $maxdenominator) { return ($hn[-1], $kn[-1]); }
-
-	#otherwise, create sequence references and pass one level deeper
-	@an = (@an, $newa);
-	@hn = (@hn, $newh);
-	@kn = (@kn, $newk);
-	my $anref = \@an;
-	my $hnref = \@hn;
-	my $knref = \@kn;
-	return contFrac($x, $maxdenominator, an => $anref, hn => $hnref, kn => $knref);
+	my $n    = int($step);
+	my ($h0, $h1, $k0, $k1) = (1, $n, 0, 1);
+
+	# End when $step is an integer.
+	while ($step != $n) {
+		$step = 1 / ($step - $n);
+
+		# Compute the next integer from the continued fraction sequence.
+		$n = int($step);
+
+		# Compute the next numerator and denominator according to the continued fraction formulas.
+		my ($newh, $newk) = ($n * $h1 + $h0, $n * $k1 + $k0);
+
+		# Machine rounding error may begin to make denominators skyrocket out of control
+		last if ($newk > $maxdenominator);
+
+		($h0, $h1, $k0, $k1) = ($h1, $newh, $k1, $newk);
+	}
 
+	return ($h1, $k1);
 }
 
 #
diff --git a/t/contexts/fraction.t b/t/contexts/fraction.t
index d873bba78e..6019df0ce0 100644
--- a/t/contexts/fraction.t
+++ b/t/contexts/fraction.t
@@ -21,9 +21,24 @@ import Parser::Legacy;
 
 Context('Fraction');
 
-ok my $a1 = Compute('1/2');
-ok my $a2 = Compute('2/4');
-
-is $a1->value, $a2->value, 'contextFraction: reduce fractions';
+subtest 'contextFraction: Basic computation and reduction' => sub {
+	ok my $a1 = Compute('1/2'), 'compute 1/2';
+	ok my $a2 = Compute('2/4'), 'compute 2/4';
+
+	is $a1->value, $a2->value, 'comparison (1/2 = 2/4)';
+};
+
+subtest 'contextFraction: Conversion of real to fraction' => sub {
+	my ($result, $direct);
+	for my $num (1 .. 100) {
+		for my $den (1 .. 100) {
+			my $real = Real($num / $den);
+			push(@$result, Fraction($real)->value);
+			push(@$direct, Fraction($num, $den)->value);
+		}
+	}
+
+	is $result, $direct, 'converted real gives correct fraction';
+};
 
 done_testing();