# Continued Logarithm
A variant of continued fractions with simpler calculations.
# 1 Concept
Introduced at the end of:
Bill Gosper, “Continued Fraction Arithmetic”, unpublished, circa 1978 – https://perl.plover.com/yak/cftalk/INFO/gosper.txt
There is a nice explanation in:
T Brabec, “Hardware Implementation of Continued Logarithm Arithmetic,” 12th GAMM - IMACS International Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics (SCAN 2006), Duisburg, Germany, 2006, pp. 9-9, doi: 10.1109/SCAN.2006.23.
(Binary) continued logarithms are simpler to implement than continued fractions because the only operations needed are
- addition
- subtraction
- multiplication by 2 (a simple bit shift in binary)
- comparison
# 1.1 Digits
The first digit of \(x\) is:
-
\(-\) when \(x \in [-\infty, 0)\), followed by digits of \(- x\).
-
\(/\) when \(x \in [0, 1)\), followed by digits of \(1 / x\).
-
\(0\) when \(x \in [1, 2)\), followed by digits of \(1 / (x - 1)\).
-
\(1\) when \(x \in [2, \infty)\), followed by digits of \(x / 2\).
-
\(\infty\) when \(x = \infty\), terminating the sequence.
\(-\) and \(/\) occur at most once (and in that order), at the start of the expansion.
\(\infty\) occurs at most once, at the end of the expansion of a rational number.
Infinity can also be represented by an infinite sequence of \(111\dots\), but (hopefully) rational operations on rational numbers will always give the terminating form.
Another ambiguity is that \(00\infty = 10\infty\), the second form should be preferred.
# 2 Implementation
# 2.1 JavaScript
My first implementation in JavaScript was doomed because values were represented by stateful coroutines, and I had no way to clone/duplicate them to allow evaluating expressions that have variables occuring more than once, never mind doing that while sharing computational effort.
# 2.2 Haskell
Haskell’s lazy evaluation model proved perfect, providing streaming computation with sharing. Computation is driven by demand for output digits.
My implementation is at code.mathr.co.uk/continued-logarithms:
git clone https://code.mathr.co.uk/continued-logarithms.git
Tested with ghc
and hugs
(Hugs is 50x slower but works on smaller systems).
A Hackage release is pending further testing.
# 3 Algorithms
Algorithms involve homographic functions, that is, functions with form invariant under forward / backward composition with Möbius transformations.
Let Hom a b c d
represent the function:
\[\frac{a x + b}{c x + d}.\]
Let BiHom a b c d e f g h
represent the two variable function:
\[\frac{a x y + b y + c x + d}{e x y + f y + g x + h}.\]
Let SqrHom a b c d e f
represent the function:
\[\frac{a x^2 + b x + c}{d x^2 + e x + f}.\]
The digits have actions:
\(/\)
Hom 0 1 1 0
\(-\)
Hom (-1) 0 0 1
\(0\)
Hom 1 1 1 0
\(1\)
Hom 2 0 0 1
# 3.1 Hom
Let Hom a b c d
represent the function:
\[\frac{a x + b}{c x + d}\]
Input of a digit whose action is Hom x y z w
is by matrix multiplication on the right:
input (Hom x y z w) (Hom a b c d) = Hom
(a * x + b * z) (a * y + b * w)
(c * x + d * z) (c * y + d * w)
If the digit is \(\infty\), output a constant (likely \(\frac{a}{c}\)).
Output of a digit whose action is Hom x y z w
is by matrix multiplication by the inverse of the action on the left:
output (Hom x y z w) (Hom a b c d) = Hom
( w * a - y * c) ( w * b - y * d)
(-z * a + x * c) (-z * b + x * d)
If the digit is \(\infty\), the computation halts.
Evaluation of \[\frac{a x + b}{c x + d}\] where \(a, b, c, d\) are integers and \(x\) is a continued logarithm proceeds as follows:
If the first digit of \(x\) is \(-\), input it.
If the first digit of \(x\) is \(/\), input it.
Now it can be assumed that \(x \ge 1\).
Symbolic differentiation shows that there are stationary points only when \(ad - bc = 0\), in which case the function is constant, so evaluating at endpoints is fine.
There is a vertical asymptote when \(x = -\frac{d}{c}\), which will cause problems if it is bigger than \(1\), which happens if \(c\) and \(c + d\) have differing signs.
If the asymptote is bigger than \(1\), input the next digit of \(x\).
Otherwise we know the function is monotonic on \([1, \infty]\), so we can find the interval of its image by computing the interval of the image of the endpoints, which are \(\frac{a + b}{c + d}\) and \(\frac{a}{c}\).
If both of these are in \([-\infty, 0)\), output \(-\).
If both of these are in \([0,1)\), output \(/\).
If both of these are \([1,2)\), output \(0\).
If both of these are \([2,\infty)\), output \(1\).
If both of these are \(\infty\), output \(\infty\) and halt.
Otherwise we cannot reach any conclusion yet, so input the next digit of \(x\).
# 3.2 BiHom
Assume \(x,y \ge 1\), otherwise input their leading digits.
By the argument above for Hom
,
if stationary points exist for any argument,
then the function is constant with respect to that argument,
so evaluation at endpoints of the input interval suffices.
There are vertical asymptotes when \(e x y + f y + g x + h = 0\).
Fix \(x\), then \((e x + f) y + (g x + h)\) is linear in \(y\) and thus crosses \(0\) unless \(e x + f = 0\). The zero crossing is at \(y = -(g x + h) / (e x + f)\), which is bigger than \(1\) when \(e + g\) and \(e + f + g + h\) have differing signs.
Similarly fixing \(y\), then the zero crossing is bigger than \(1\) when \(e + f\) and \(e + f + g + h\) have differing signs.
That is, all of \(e + f, e + g, e + f + g + h\) must have the same sign for the function to be monotonic. Furthermore, \(e\) must have the same sign too, considering there must be no sign changes when both \(x,y \to \infty\) and \(x, y \to 1\).
Now the monotonic function can be evaluated at its endpoints giving:
\[\frac{a}{e}, \frac{a + b}{e + f}, \frac{a + c}{e + g}, \frac{a + b + c + d}{e + f + g + h}.\]
Compare this interval with the respective digit intervals and output if possible, otherwise input the next digit of \(x\) and/or \(y\).
# 3.3 SqrHom
Evaluation of \[\frac{a x^2 + b x + c}{d x^2 + e x + f}\] where \(a, b, c, d, e, f\) are integers and \(x\) is a continued logarithm proceeds as follows:
If the first digit of \(x\) is \(-\), input it.
If the first digit of \(x\) is \(/\), input it.
Now it can be assumed that \(x \ge 1\).
Symbolic differentiation shows that there are stationary points when \[(ae-bd) x^2 + (2(af-cd)) x + (bf-ce) = 0 = A x^2 + B x + C.\]
If \(A = B = C = 0\) then the function is constant.
If \(A = B = 0 \ne C\) then there are no stationary points.
If \(A = 0 \ne B\) then the stationary point is \(x = -C / B\) which is only a problem if it’s bigger than \(1\), that is when \(B\) and \(B + C\) have differing signs.
If \(A \ne 0\) then there are stationary points if \(D = B^2 - 4 A C \ge 0\). The largest stationary point is bigger than \(1\) when \((2 A + B)^2 - D\) and \(A\) have differing signs.
There are vertical asymptotes when \(d x^2 + e x + f = 0\).
If \(d = e = f = 0\) then the function is infinite.
If \(d = e = 0 \ne f\) then there are no asymptotes.
If \(d = 0 \ne e\) then the asymptote is at \(x = -f / e\) which is bigger than \(1\) when \(e\) and \(e + f\) have differing signs.
If \(d \ne 0\) then there are asymptotes if \(E = e^2 - 4 d f \ge 0\). The largest asymptote is bigger than \(1\) when \((2 d + e)^2 - E\) and \(d\) have differing signs.
If there are stationary points or asymptotes bigger than \(1\), input the next digit of \(x\).
Otherwise, the function is monotonic on \([1, \infty]\), and evaluating at the endpoints gives \(\frac{a + b + c}{d + e + f}\) and \(\frac{a}{d}\).
Compare this interval with the respective digit intervals and output if possible, otherwise input the next digit of \(x\).
Possibly avoiding SqrHom
and using BiHom
expressions
with duplicated input variables might be more efficient,
avoiding the integer multiplications,
or I am doing something wrong in my analysis
and a version without integer multiplication is possible.
# 3.4 Arithmetic
\(+\), \(-\), \(*\), \(/\) can be performed
with suitable BiHom
evaluations.
# 3.5 Merge
Suppose one can construct a list of pairs of lower and upper bounds (expressed as continued logarithms), with the lower bounds increasing and the upper bounds decreasing, such that if any pair of lower and upper bounds agree to N digits, then all the subsequent pairs of bounds agree to at least N digits too.
If this list is finite assume it ends with an equal pair, otherwise assume that the bounds always eventually get closer together (otherwise the algorithm may not terminate).
Then from it one can construct a single continued logarithm such that it is between the bounds.
# 3.6 Inverse
Given a monotonic strictly increasing function \(f\) from \(I \subset [0, \infty]\) to \(J = f(I)\), its inverse \(f^{-1}\) can be found by binary search to emit a sequence of bounds to be merged.
# 3.7 Logarithm
Base \(n\) logarithm can be performed by binary search and merge using Farey addition: \(\log_n x = \frac{a}{b}\) means \(\frac{x^b}{n^a} = 1\). Set up the search bounds like \[\frac{x^b}{n^a} > 1 > \frac{x^d}{n^c}\] with \(a = 0, b = 1, c = 1, d = 0\), that is, \(x > 1 > \frac{1}{n}\). Proceed by \(e = a + c, f = b + d\), and consider \[\frac{x^f}{n^e} = \frac{x^b}{n^a} \times \frac{x^d}{n^c}.\]
# 3.8 Exponential
\(x^y\) can be performed by inverting \(\log_x\).
# 3.9 e
Natural logarithm and exponential can be found using \(e\), computed using a merge of bounds from its continued fraction.
# 3.10 Hyperbolic
Can be implemented in terms of exp.
Inverse hyperbolic via logs or function inversion.
# 3.11 Trigonometric
Inverse tangent can be computed using a merge of bounds from the generalized continued fraction.
Other inverse trigonometry in terms of inverse tangent.
Tangent via inverting inverse tangent after range reduction using double-angle identity.
Sine and cosine in terms of tangent.
# 3.12 pi
\(\pi\) can be computed using inverse tangent.
# 3.13 Conversions
Conversion to a tightening sequence of lower/upper bounds of generalized rationals (allowing 0 denominator) is simple: after the sign, each tail is in \([1, \infty]\) so bounds can be calculated.
Conversion to a rational (when it terminates) is by merging the above.
Conversion to integer (e.g. floor, ceiling) is by merging bounds: if floor lower = floor upper then output the common value, otherwise continue (possibly not terminating).
# 4 Evaluation
# 4.1 Performance
Benchmarked using Haskell compiled with ghc -O2
of an unzoomed Mandelbrot at 60x60 pixels with 100 iterations limit.
-
Double is more or less instant.
-
Continued Logarithm takes about 4m10s.
-
Rational did not finish; it would have exhausted memory, had I not interrupted it after 20mins consuming 1GB. See exact Mandelbrot.
# 4.2 Accuracy
Continued Logarithm should be as good as exact rational, with an implicit rational interval representation for other numbers that can get as narrow as needed.
A True
result of the comparison norm z > 4
means that \(|z|^2\) is certainly bigger than \(4\),
(calculating only as many “digits” as needed
for the lower bound on \(|z|^2\) to be bigger than \(4\)).
# 4.3 Termination
Rationals have terminating expansion.
If all input is rational, and only rational arithmetic is used (i.e. no square roots or similar operations), then the result will be rational with a terminating representation.
If irrational arithmetic is used, all bets are off: determining even the first digit of \[(\sqrt{2})\times(\sqrt{2}) = {2}_{10} = {10\infty}_{CL(2)}\] is probably impossible.
So, there remains a gap for a representation that can always output useful information, for example telling that \(\sqrt{2}^2 \approx 2\), if not \(\sqrt{2}^2 = 2\).