<commit_before>
		Continued fractions

	    Approximation of numbers by continued fractions

The function {ContFrac} converts a (rational) number $r$ into a regular
continued fraction,
$$ r = n[0] + 1/(n[1] + 1/(n[2]+...)) $$.
Here all numbers $n[i]$ ("terms" of a continued fraction) are integers and all except $n[0]$ must be positive.
(Continued fractions may not converge unless their terms are positive and
bounded from below.)

The algorithm for converting a rational number $r=n/m$ into a continued
fraction is simple. First, we determine the integer part of $r$, which is
$Div(n,m)$. If it is negative, we need to subtract one, so that $r=n[0]+x$ and
the remainder $x$ is nonnegative and less than $1$. The remainder
$x=Mod(n,m)/m$ is then inverted, $r[1] := 1/x = m/Mod(n,m)$ and so we have
completed the first step in the decomposition, $r = n[0] + 1/r[1]$; now $n[0]$
is integer but $r[1]$ is perhaps not integer. We repeat the same procedure on
$r[1]$, obtain the next integer term $n[1]$ and the remainder $r[2]$ and so on,
until such $n$ that $r[n]$ is an integer and there is no more work to do. This
process will always terminate because all floating-point values are actually
rationals in disguise.

Continued fractions are useful in many ways. For example, if we know that a
certain number $x$ is rational but have only a floating-point representation of
$x$ with a limited precision, say, 1.5662650602409638, we can try to guess its
rational form (in this example $x=130/83$). The function {GuessRational} uses
continued fractions to find a rational number with "optimal" (small) numerator
and denominator that is approximately equal to a given floating-point number.

Consider the following example. The number $17/3$ has a continued fraction
expansion {{5,1,2}}. Evaluated as a floating point number with limited
precision, it may become something like $17/3+0.00001$, where the small number
represents a roundoff error. The continued fraction expansion of this number is
{{5, 1, 2, 11110, 1, 5, 1, 3, 2777, 2}}. The presence of an unnaturally large
term $11110$ clearly signifies the place where the floating-point error was
introduced; all terms following it should be discarded to recover the continued fraction {{5,1,2}} and from it the initial number $17/3$.

If a continued fraction for a number $x$ is
cut right before an unusually large term, and
evaluated, the resulting rational number is very close to close to $x$ but has an unusually small
denominator. This works because partial continued fractions provide "optimal"
rational approximations for the final (irrational) number, and because the
magnitude of the terms of the partial fraction is related to the magnitude of
the denominator of the resulting rational approximation.

{GuessRational(x, prec)} needs to choose the place where it should cut the
continued fraction. The algorithm for this is somewhat less precise than it
could be but it works well enough. The idea is to cut the continued fraction
when adding one more term would change the result by less than the specified
precision. To realize this in practice, we need an estimate of how much a
continued fraction changes when we add one term.

The routine {GuessRational} uses an easy upper bound for the difference of continued fractions that differ only by an additional last term:
$$ Abs(delta) := Abs( 1/(a[1]+1/(...+1/a[n])) - 1/(a[1]+1/(...+1/a[n+1])) ) < 1/ ((a[1]*...*a[n])^2 * a[n+1]) $$.
Thus we should compute the product of successive terms $a[i]$ of the continued fraction and stop at $a[n]$ at which this product exceeds the maximum number of digits. The routine {GuessRational} has a second parameter {prec} which is by default 1/2 times the number of decimal digits of current precision; it stops at $a[n]$ at which the product $a[1]*...*a[n]$ exceeds $10^prec$.

The above estimate for $delta$ hinges on the inequality
$$ 1/(a+1/(b+...)) < 1/a $$
and is suboptimal if some terms $a[i]=1$, because the product of $a[i]$ does not increase when one of the terms is equal to 1, whereas in fact these terms do make $delta$ smaller. A somewhat better estimate would be obtained if we use the inequality
$$ 1/(a+1/(b+1/(c+...))) < 1/(a+1/(b+1/c)) $$.
This does not lead to a significant improvement if $a>1$ but makes a difference when $a=1$. In the product $a[1]*...*a[n]$, the terms $a[i]$ which are equal to 1 should be replaced by
$$ a[i]+1/(a[i+1]+1/a[i+2])$$.
Since the comparison of $a[1]*...*a[n]$ with $10^prec$ is qualitative, it it enough to do calculations for it with limited precision.

This algorithm works well if $x$ is computed with enough precision; namely, it
must be computed to at least as many digits as there are in the numerator and
the denominator of the fraction combined. Also, the parameter {prec} should not
be too large (or else the algorithm will find another rational number with a larger
denominator that approximates $x$ "better" than the precision to which you know $x$).

The related function {NearRational(x, prec)} works somewhat differently. The
goal is to find an "optimal" rational number, i.e. with smallest numerator and
denominator, that is within the distance $10^(-prec)$ of a given value $x$. The
algorithm for this comes from the 1972 HAKMEM document, Item 101C. Their
description is terse but clear:

	Problem: Given an interval, find in it the
	rational number with the smallest numerator and
	denominator.
	Solution: Express the endpoints as continued
	fractions.  Find the first term where they differ
	and add 1 to the lesser term, unless it's last. 
	Discard the terms to the right.  What's left is
	the continued fraction for the "smallest"
	rational in the interval.  (If one fraction
	terminates but matches the other as far as it
	goes, append an infinity and proceed as above.)

The HAKMEM text (<NAME>, <NAME>, and <NAME>: Memo No. 239, MIT AI Lab, 1972, available as HTML online from various places) contains several interesting insights relevant to continued fractions and other numerical algorithms.

	    Precision of approximation by continued fractions

Sometimes an analytic function $f(x)$ can be approximated using a continued
fraction that contains $x$ in its terms. Examples include: the inverse tangent
$ArcTan(x)$, the error function $Erf(x)$ and the incomplete gamma function
$Gamma(a,x)$ (see below for details). For these functions, continued fractions
provide a method of numerical calculation that works when the Taylor series
converges slowly or not at all. However, continued fractions usually converge
quickly for one value of $x$ but slowly for another. Also, it is not as easy to
obtain an analytic error bound for a continued fraction approximation as it is
for power series.

In this section we describe two methods that compute a continued fraction: the
simple "bottom-up" and the more complicated "top-down" method. The "bottom-up"
method is faster but requires to know the number of terms in advance. The
"top-down" method provides an automatic error estimate and can be used to
evaluate a continued fraction with more and more terms until the desired
precision is achieved. The formula for the precision of the continued fraction
approximation used in the "top-down" method sometimes allows to estimate the
number of terms in advance.

*REM FIXME

		Newton's method and its improvements

The Newton-Raphson method of numerical solution of algebraic equations can be used to obtain multiple-precision values of several elementary functions.

The basic formula is widely known: If $f(x)=0$ must be solved, one starts with a value of $x$ that is close to some root and iterates $$x'=x-f(x)*(D(x)f(x))^(-1)$$.
This formula is based on the approximation of the function $f(x)$ by a tangent line at some point $x$. A Taylor expansion in the neighborhood of the root shows that (for an initial value $x[0]$ sufficiently close to the root) each iteration gives at least twice as many correct digits of the root as the previous one ("quadratic convergence"). Therefore the complexity of this algorithm is proportional to a logarithm of the required precision and to the time it takes to evaluate the function and its derivative. Generalizations of this method require computation of higher derivatives of the function $f(x)$ but successive approximations to the root converge several times faster (the complexity is still logarithmic).<commit_msg>some more docs about the continued fractions, before I forget this.<commit_after>
		Continued fractions

	    Approximation of numbers by continued fractions

The function {ContFrac} converts a (rational) number $r$ into a regular
continued fraction,
$$ r = n[0] + 1/(n[1] + 1/(n[2]+...)) $$.
Here all numbers $n[i]$ ("terms" of a continued fraction) are integers and all except $n[0]$ must be positive.
(Continued fractions may not converge unless their terms are positive and
bounded from below.)

The algorithm for converting a rational number $r=n/m$ into a continued
fraction is simple. First, we determine the integer part of $r$, which is
$Div(n,m)$. If it is negative, we need to subtract one, so that $r=n[0]+x$ and
the remainder $x$ is nonnegative and less than $1$. The remainder
$x=Mod(n,m)/m$ is then inverted, $r[1] := 1/x = m/Mod(n,m)$ and so we have
completed the first step in the decomposition, $r = n[0] + 1/r[1]$; now $n[0]$
is integer but $r[1]$ is perhaps not integer. We repeat the same procedure on
$r[1]$, obtain the next integer term $n[1]$ and the remainder $r[2]$ and so on,
until such $n$ that $r[n]$ is an integer and there is no more work to do. This
process will always terminate because all floating-point values are actually
rationals in disguise.

Continued fractions are useful in many ways. For example, if we know that a
certain number $x$ is rational but have only a floating-point representation of
$x$ with a limited precision, say, 1.5662650602409638, we can try to guess its
rational form (in this example $x=130/83$). The function {GuessRational} uses
continued fractions to find a rational number with "optimal" (small) numerator
and denominator that is approximately equal to a given floating-point number.

Consider the following example. The number $17/3$ has a continued fraction
expansion {{5,1,2}}. Evaluated as a floating point number with limited
precision, it may become something like $17/3+0.00001$, where the small number
represents a roundoff error. The continued fraction expansion of this number is
{{5, 1, 2, 11110, 1, 5, 1, 3, 2777, 2}}. The presence of an unnaturally large
term $11110$ clearly signifies the place where the floating-point error was
introduced; all terms following it should be discarded to recover the continued fraction {{5,1,2}} and from it the initial number $17/3$.

If a continued fraction for a number $x$ is
cut right before an unusually large term, and
evaluated, the resulting rational number is very close to close to $x$ but has an unusually small
denominator. This works because partial continued fractions provide "optimal"
rational approximations for the final (irrational) number, and because the
magnitude of the terms of the partial fraction is related to the magnitude of
the denominator of the resulting rational approximation.

{GuessRational(x, prec)} needs to choose the place where it should cut the
continued fraction. The algorithm for this is somewhat less precise than it
could be but it works well enough. The idea is to cut the continued fraction
when adding one more term would change the result by less than the specified
precision. To realize this in practice, we need an estimate of how much a
continued fraction changes when we add one term.

The routine {GuessRational} uses an easy upper bound for the difference of continued fractions that differ only by an additional last term:
$$ Abs(delta) := Abs( 1/(a[1]+1/(...+1/a[n])) - 1/(a[1]+1/(...+1/a[n+1])) ) < 1/ ((a[1]*...*a[n])^2 * a[n+1]) $$.
Thus we should compute the product of successive terms $a[i]$ of the continued fraction and stop at $a[n]$ at which this product exceeds the maximum number of digits. The routine {GuessRational} has a second parameter {prec} which is by default 1/2 times the number of decimal digits of current precision; it stops at $a[n]$ at which the product $a[1]*...*a[n]$ exceeds $10^prec$.

The above estimate for $delta$ hinges on the inequality
$$ 1/(a+1/(b+...)) < 1/a $$
and is suboptimal if some terms $a[i]=1$, because the product of $a[i]$ does not increase when one of the terms is equal to 1, whereas in fact these terms do make $delta$ smaller. A somewhat better estimate would be obtained if we use the inequality
$$ 1/(a+1/(b+1/(c+...))) < 1/(a+1/(b+1/c)) $$.
This does not lead to a significant improvement if $a>1$ but makes a difference when $a=1$. In the product $a[1]*...*a[n]$, the terms $a[i]$ which are equal to 1 should be replaced by
$$ a[i]+1/(a[i+1]+1/a[i+2])$$.
Since the comparison of $a[1]*...*a[n]$ with $10^prec$ is qualitative, it it enough to do calculations for it with limited precision.

This algorithm works well if $x$ is computed with enough precision; namely, it
must be computed to at least as many digits as there are in the numerator and
the denominator of the fraction combined. Also, the parameter {prec} should not
be too large (or else the algorithm will find another rational number with a larger
denominator that approximates $x$ "better" than the precision to which you know $x$).

The related function {NearRational(x, prec)} works somewhat differently. The
goal is to find an "optimal" rational number, i.e. with smallest numerator and
denominator, that is within the distance $10^(-prec)$ of a given value $x$. The
algorithm for this comes from the 1972 HAKMEM document, Item 101C. Their
description is terse but clear:

	Problem: Given an interval, find in it the
	rational number with the smallest numerator and
	denominator.
	Solution: Express the endpoints as continued
	fractions.  Find the first term where they differ
	and add 1 to the lesser term, unless it's last. 
	Discard the terms to the right.  What's left is
	the continued fraction for the "smallest"
	rational in the interval.  (If one fraction
	terminates but matches the other as far as it
	goes, append an infinity and proceed as above.)

The HAKMEM text (<NAME>, <NAME>, and <NAME>: Memo No. 239, MIT AI Lab, 1972, available as HTML online from various places) contains several interesting insights relevant to continued fractions and other numerical algorithms.

	    Accurate computation of continued fractions

Sometimes an analytic function $f(x)$ can be approximated using a continued
fraction that contains $x$ in its terms. Examples include: the inverse tangent
$ArcTan(x)$, the error function $Erf(x)$ and the incomplete gamma function
$Gamma(a,x)$ (see below for details). For these functions, continued fractions
provide a method of numerical calculation that works when the Taylor series
converges slowly or not at all. However, continued fractions usually converge
quickly for one value of $x$ but slowly for another. Also, it is not as easy to
obtain an analytic error bound for a continued fraction approximation as it is
for power series.

In this section we describe two methods that compute a continued fraction: the
simple "bottom-up" and the more complicated "top-down" method. The "bottom-up"
method is faster but requires to know the number of terms in advance. The
"top-down" method provides an automatic error estimate and can be used to
evaluate a continued fraction with more and more terms until the desired
precision is achieved. The formula for the precision of the continued fraction
approximation used in the "top-down" method sometimes allows to estimate the
number of terms in advance.

A continued fraction
$$ a[0]+b[0]/(a[1]+b[1]/(a[2]+...))$$
is specified by a set of <i>terms</i> ($a[i]$, $b[i]$). [If continued fractions are used to approximate analytic functions such as $ArcTan(x)$, then ($a[i]$, $b[i]$) will be some infinite sequences that depend on $x$.]
Denote the truncated fraction containing only terms from $m$ to $n$, i.e.
$$ a[m]+b[m]/(a[m+1]+b[m+1]/(...+b[n]/a[n]))$$,
by $F[m][n]$.
In this notation, the continued fraction that we need to compute is $F[0][n]$. Our problem is first, to select a large enough $n$ so that $F[0][n]$ gives enough precision, and second, to compute the value of $F[0][n]$.

The "bottom-up" method is to select some value of $n$ and start evaluating the fraction from the bottom. First we take $F[n][n]=a[n]$ and then we use the obvious relation of backward recurrence,
$$ F[m][n] = a[m]+b[m]/F[m+1][n] $$,
to obtain successively $F[n-1][n]$, ..., $F[0][n]$. This method requires one long division at each step.

An alternative implementation may be to imagine that all $a[i]$, $b[i]$ are
integers, and to obtain the numerator and the denominator of $F[0][n]$
separately as two simultaneous backward recurrences. If 
$F[m+1][n] = p[m+1]/q[m+1]$, then 
$p[m]=a[m]*p[m+1]+b[m]*q[m+1]$ and $q[m]=p[m+1]$. This method would start with
$p[n]=a[n]$, $q[n]=1$ and require two long multiplications at each step; the only
division will be performed at the very end. Sometimes this method allows to
reduce the roundoff error.

There is an improvement to this method that can sometimes increase the achieved
precision without computing more terms. Namely, we choose for $F[0][n]$ not the
obvious $a[n]$ but another number that more closely approximates the remainder
of the fraction which we are not computing.
*REM FIXME

The "bottom-up" method obviously requires to know the number of terms $n$ in
advance; calculations have to be started all over again if more terms are
needed. Also, it provides no error estimate.

The "top-down" method is slower and is based on the error estimate which we shall now derive. We need to estimate the difference between successive approximations $F[0][n]$ and $F[0][n+1]$.
The basic recurrence relation is
$$ F[m][n+1] - F[m][n] = -(b[m]*(F[m+1][n+1]-F[m+1][n]))/(F[m+1][n+1]*F[m+1][n]) $$.

*REM FIXME

		Newton's method and its improvements

The Newton-Raphson method of numerical solution of algebraic equations can be used to obtain multiple-precision values of several elementary functions.

The basic formula is widely known: If $f(x)=0$ must be solved, one starts with a value of $x$ that is close to some root and iterates $$x'=x-f(x)*(D(x)f(x))^(-1)$$.
This formula is based on the approximation of the function $f(x)$ by a tangent line at some point $x$. A Taylor expansion in the neighborhood of the root shows that (for an initial value $x[0]$ sufficiently close to the root) each iteration gives at least twice as many correct digits of the root as the previous one ("quadratic convergence"). Therefore the complexity of this algorithm is proportional to a logarithm of the required precision and to the time it takes to evaluate the function and its derivative. Generalizations of this method require computation of higher derivatives of the function $f(x)$ but successive approximations to the root converge several times faster (the complexity is still logarithmic).
{'fraction_non_alphanumeric': 0.08687228360475467, 'fraction_numerical': 0.01665183013038697, 'mean_word_length': 4.608810572687225, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4027827', 'n_tokens_mistral': 5492, 'n_tokens_neox': 5084, 'n_words': 2925}