3

I'm trying to implement basic arithmetic on Bill Gosper's continued logarithms, which are a 'mutation' of continued fractions allowing the term co-routines to emit and consume very small messages even on very large or very small numbers.

Reversible arithmetic, such as {+,-,*,/} are fairly straightforwardly described by Gosper at least in a unary representation, but I'm having difficulty implementing the modulo operator which effectively truncates information from the division operation.

I've realized the modulo operator can be mostly defined with operations I already have:

a mod b == a - b * floor(a / b)

leaving my only problem with floor.

I've also read that the run-length encoded format for continued logarithms effectively describes

'... the integer part of the log base 2 of the number remaining to be described.'

So yielding the first term right away (pass through) produces the correct output so far, but leaves a significant portion to be determined which I assume requires some sort of carry mechanism.

I've written the following code to test input terms and the expected output terms, but I'm mainly looking for high level algorithm ideas behind implementing floor.

An example input (1234 / 5) to output pair is

Input: [7, 0, 3, 0, 0, 0, 0, 1, 3, 3, 1]

Output: [7, 0, 3, 1, 4, 2, 1, 1]

from fractions import Fraction

def const(frac):
        """ CL bistream from a fraction >= 1 or 0. """
        while frac:
                if frac >= 2:
                        yield 1
                        frac = Fraction(frac, 2)
                else:
                        yield 0
                        frac -= 1
                        frac = Fraction(1, frac) if frac else 0

def rle(bit_seq):
        """ Run-length encoded CL bitstream. """
        s = 0
        for bit in bit_seq:
                s += bit
                if not bit:
                        yield s
                        s = 0

def floor(rle_seq):
        """ RLE CL terms of the greatest integer less than rle_seq. """
        #pass
        yield from output

""" Sample input/output pairs for floor(). """
num = Fraction(1234)
for den in range(1, int(num)+1):
        input  = list(rle(const(num  / den)))
        output = list(rle(const(num // den))) # Integer division!
        print(">  ", input)
        print(">> ", output) 
        print(">>*", list(floor(input))) 
        print()
        assert(list(floor(input)) == output)

How can I implement the floor operator in the spirit of continued fraction arithmetic by consuming terms only when necessary and emitting terms right away, and especially only using the run-length encoded format (in binary) rather than the unary expansion Gosper tends to describe.

  • I am not flagging this because I lack the expertise to feel certain one way or another, but do you think this might make more sense as a question in the mathematics stack exchange? – Max von Hippel Nov 11 '17 at 22:57
  • 1
    @MaxvonHippel I'm not sure, either. I was originally going to post there but I felt because of the code sample I gave it would be more appropriate here. Feel free to move it if it's off topic here. –  Nov 11 '17 at 22:58
  • I'm not sure that the context is enormously important here (that is, the code block & explanation of how you're working on modulo). I think the question might be clearer if you just specify that you are trying to implement `floor` using this specific philosophy of mathematical implementation, link to the book (as you did), and then give the input / output (as you did) and maybe your working hypothesis. Because your question is really about `floor`, not `modulo`. – Max von Hippel Nov 11 '17 at 23:03
  • That said, as far as the site is concerned, I'd say just leave it here for a while and see how people who know more than I do about the topic respond. If the consensus is that it should be on the math site then move it; it's definitely not a bad question for SO, it just might be a better question for mathematics. – Max von Hippel Nov 11 '17 at 23:04

1 Answers1

0

By assuming that the next coefficient in the run-length encoding is infinite, you can get a lower bound. By assuming that the next term is 1, you can get an upper bound.

You can simply process as many run-length encoded coefficients until you know that both the lower and the upper bound are in the half-open interval [N, N + 1). In this case you know that the floor of the continued logarithm is N. This is similar to what Bill Gosper does at the start of the linked document.

Note, however, that this process doesn't necessarily terminate. For example, when you multiply sqrt(2) by sqrt(2), you get, of course, the number 2. However, the continued logarithm for sqrt(2) is infinite. To evaluate the product sqrt(2) * sqrt(2) you will need all the coefficients to know that you will end up with 2. With any finite number of terms, you can't decide if the product is less than 2 or at least equal to it.

Note that this problem is not specific to continued logarithms, but it is a fundamental problem that occurs in any system in which you can have two numbers for which the representation is infinite but the product can be represented with a finite number of coefficients.

To illustrate this, suppose that these coroutines don't spit out run-length encoded values, but decimal digits, and we want to calculate floor(sqrt(2) * sqrt(2)). After how many steps can we be sure that the product will be at least 2? Let's take 11 digits, just to see what happens: 1.41421356237 * 1.41421356237 = 1.9999999999912458800169

As you might guess, we get arbitrarily close to 2, but will never 'reach' 2. Indeed, without knowing that the source of the digits is sqrt(2), it might just happen that the digits terminate after that point and that the product ends up below 2. Similarly, all following digits might be 9's, which would result in a product slightly above 2.

(A simpler example would be to take the floor of a routine that produces 0.9999...)

So in these kind of arbitrary-precision numerical systems you can end up in situations where you can only calculate some interval (N - epsilon, N + epsilon), where you can make epsilon arbitrarily small, but never equal to zero. It is not possible to take the floor of this expression, as -- by the numerical methods employed -- it is not possible to decide if the real value will end up below or above N.

Ruben
  • 180
  • 1
  • 11