Contlog - representing numbers as continued logarithms

Introduction

Consider this version of the binary GCD algorithm, for finding the greatest common
divisor of two positive numbers using only shifting and subtraction.

int gcd(int pair[2])
{
	int shift = ffs(pair[0] | pair[1]) - 1;
	pair[0] >>= shift;
	pair[1] >>= shift;
	for (int i = pair[0] >= pair[1]; pair[0] != pair[1]; i ^= 1) {
		while (pair[i] << 1 < pair[i^1]) {
			pair[i] <<= 1;
			printf("%d", i);
		}
		pair[i^1] -= pair[i];
		printf("%d", i);
	}
	pair[0] >>= ffs(pair[0]) - 1;
	printf("1\n");
	return (pair[0] << shift);
}

The print statements in the function log the transformations made to
the values that ultimately made them match, to end the iteration.
They also write a binary representation of the rational number that is
the ratio of the two inputs.

In his notes on continued fraction arithmetic, Gosper introduced
continued logarithms as a binary alternative to continued fractions
for representing real numbers.  The representation here is only
slightly different from the one he proposed.  He observed that a
string of bits like

110110001

can represent the fraction

2**1 + 2**1 / (
       2**0 + 2**0 / (
              2**1 + 2**1 / (
	             2**2 + 2**2 / (
		            2**0))))

which is 44/13.  The first substring of n+1 1-bits shows that the
value is between 2**n and 2**(n+1).  After subtracting 2**n, and
dividing by 2**n, the value left is expressed in the bits that follow,
beginning with a string of 0-bits.  Each set of matching bits
represents the integer part of a binary logarithm, thus the name
'continued logarithms'.

The representation here has the virtue that the mapping from a set of
bits, interpreted as a binary integer, to a rational number, is
monotonic.  That is, if integer i maps to rational number p, and
integer j maps to rational number q, then i < j if and only if p < q.
This is not a characteristic of Gosper's initially proposed
representation.

Strictly speaking, Gosper's encoding only represented numbers greater
than or equal to one.  The representation presented here can represent
proper fractions by beginning the bitstring with 0-bits instead of
1-bits.  It can also include a twos-complement representation for
negative values.


Extracting ratios from this representation

While this representation is great for representing exactly rational
numbers with small numerator and denominator in just a few bits, it
sacrifices the representation of large numbers, especially ones near
powers of two.  For example 127/1 is represented by

1111111011111101111101111011101101

so the signed contlog representation of 127 requires 36 bits.
Computing a 32-bit approximation to 127/1 requires rounding; this
implementation rounds to nearest, or to even integer representations
in case of a tie.

Translating that truncated bit string back to a rational number, leads
to the result 37722176/297025, which is not likely to be an expected
result.  To provide a more expected result, this implementation
considers a fixed size bit string to represent the interval that
includes all the ratios that would be converted to that bit string,
and picks the one to represent the interval that has the smallest
numerator and denominator.  For a fixed-size bit string that ends in a
zero, that includes the values at the endpoints of the interval, and
for the others, it does not.  Otherwise, two consecutive bit strings
could map to the same ratio.  For the case of 127/1, this produces the
expected result 127/1.  For 1000/999, the result looks like 1000/999,
and not 470999/470528.

Arithmetic

Where Gosper seeks to compute with unbounded continued fractions, this
implementation aspires only to fixed-width operands and results.
Thus, it is feasible to completely decode one operand into a numerator
and denominator before considering the other operand at all, and this
implementation does that.  Those two values are scaled up, doubling
both until one or both of them would overflow if doubled again.  The
values are used to form a homographic function that is evaluated at
the value of the other operand to compute the composition of the two
operands.  So, for example, to compute x + y, first find that y=n/d, and then form

(n + dx) / (d + 0x)

and update the four coefficients as the program consumes the leading
bits of x.  When such an update leads to overflow, scale the
coefficents down, to fit within the fixed number of bits available.

After every incremental update, examine the coefficients to see if
some of the output bits can be extracted.


Test program

The 'contlog' program produced from this source code is a simple
calculator that uses the contlog representation, and allows both
integer ratios and hexadecimal bitstrings to be used to represent
operands, and displays results in both forms, and also in a familiar
floating point form.  Sample usage:

$ contlog 4/7
x: 4/7 (26000000) =     0.571428571429
$ contlog 4/7 - 5/9
x: 4/7 (26000000) =     0.571428571429
y: 5/9 (24000000) =     0.555555555556
x-y: 1/63 (01042260) =     0.015873015873
$ contlog 55555555
x: 2178309/1346269 (55555555) =     1.618033988750
$ contlog 2/1 /
x: 2/1 (60000000) =     2.000000000000
sqrt(x): 8119/5741 (4e38e38e) =     1.414213551646
x/sqrt(x): 8119/5741 (4e38e38e) =     1.414213551646