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