Bilinear Fractional Transformations

A bilinear fractional transformation (BiLFT) is a two-argument function of the form

(lambda (x y)
  (/ (+ (* A x y) (* B x) (* C y) D)
     (+ (* E x y) (* F x) (* G y) H)))

William Gosper figured out how to use BiLFTs to perform linear fractional operations on infinite compositions of LFTs. It isn’t as hard as it might seem.

BiLFTs do not form a group under functional composition with LFTs, but they have this interesting property: if you hold either input to the BiLFT constant, the BiLFT becomes a LFT on the other input. So if we consider just one of the inputs at a time (or the output), we can apply some group-like operations.

You can compose LFTs with BiLFTs in three ways. First, you can run the output of a BiLFT into the input of a LFT. Second and third, you can run the output of a LFT into the either the x or y input of a BiLFT. Composing a LFT with a BiLFT in any of these ways produces another BiLFT, so there is some notion of closure under composition. Composing the identity LFT with a BiLFT in any of these ways does not change anything. For any composition of a LFT with a BiLFT, you can compose the inverse of the LFT to undo the composition.

BiLFTs can also form infinite compositions, but since there are two inputs, each input can be infinite. We cannot use a Scheme stream, but we can create a two-tailed stream-like object called a binary-expression. A binary-expression is a composition of two infinite compositions. We represent a binary-expression as an object with two delayed cdrs.

(defclass binary-expression ()
  ((generation :initarg :generation
               :initform 0)
   (bilft :initarg :bilft
          :initform (error "Required initarg :bilft was omitted."))
   (delayed-left :initarg :delayed-left
                 :initform (error "Required initarg :delayed-left was omitted."))
   (delayed-right :initarg :delayed-right
                  :initform (error "Required initarg :delayed-right was omitted."))))

Like a LFT, we use binary-expressions by trying to operate on the BiLFT, but forcing one of the tails and refining the binary-expression when necessary.

(defun refine-left (binary-expression)
  (let ((delayed-left (slot-value binary-expression 'delayed-left)))
    (if (null delayed-left)
        (error "Attempt to refine-left on empty stream.")
        (let ((left (force delayed-left)))
          (make-instance 'binary-expression
                         :generation (+ (slot-value binary-expression 'generation) 1)
                         :bilft (compose-bilft-lft-x
                                 (slot-value binary-expression 'bilft)
                                 (if (empty-stream? left)
                                     (make-lft 1 1
                                               0 0)
                                     (stream-car left)))
                         :delayed-left (if (empty-stream? left)
                                           nil
                                           (stream-delayed-cdr left))
                         :delayed-right (slot-value binary-expression 'delayed-right))))))

(defun refine-right (binary-expression)
  (let ((delayed-right (slot-value binary-expression 'delayed-right)))
    (if (null delayed-right)
        (error "Attempt to refine-right on empty stream.")
        (let ((right (force delayed-right)))
          (make-instance 'binary-expression
                         :generation (+ (slot-value binary-expression 'generation) 1)
                         :bilft (compose-bilft-lft-y
                                 (slot-value binary-expression 'bilft)
                                 (if (empty-stream? right)
                                     (make-lft 1 1
                                               0 0)
                                     (stream-car right)))
                         :delayed-left (slot-value binary-expression 'delayed-left)
                         :delayed-right (if (empty-stream? right)
                                            nil
                                            (stream-delayed-cdr right)))))))

refine-fairly alternates forcing the delayed-left or delayed-right tail while refine-disjoint looks at the overlap in the ranges of the BiLFT.

(defun refine-fairly (binary-expression)
  (if (zerop (mod (slot-value binary-expression 'generation) 2))
      (if (null (slot-value binary-expression 'delayed-left))
          (refine-right binary-expression)
          (refine-left binary-expression))
      (if (null (slot-value binary-expression 'delayed-right))
          (refine-left binary-expression)
          (refine-right binary-expression))))

(defun refine-disjoint (binary-expression)
  (if (bilft-disjoint? (slot-value binary-expression 'bilft))
      (refine-right binary-expression)
      (refine-left binary-expression)))

You can do arithmetic on infinite compositions by using the appropriate BiLFT:

(defparameter bilft-add 
  (make-bilft 0 1 1 0
              0 0 0 1))

(defparameter bilft-subtract 
  (make-bilft 0 1 -1 0
              0 0 0 1))

(defparameter bilft-multiply 
  (make-bilft 1 0 0 0
              0 0 0 1))

(defparameter bilft-divide 
  (make-bilft 0 1 0 0
              0 0 1 0))

Converting a binary-expression to a LFT stream

binary-expressions can be operated on in an analagous way to an infinite composition of LFTs, but in order to nest binary expressions, we need to be able to turn one into an infinite composition of LFTs so it can become the left or right input of the other. The way to do this is to generate a stream of LFTs and their inverses. The inverses are composed into the output of the binary-expression while the LFTs are composed downstream into one of the inputs of the next binary-expression.

The math works such that we can select any LFT that has an inverse, but we want our binary-expression to represent a narrowing transformation, so we try composing a few different LFTs and their inverses and see if we still have a narrowing transformation. If we find one, we proceed with the computation, but if we cannot find a LFT and its inverse that preserves the narrowing, we refine the binary-expression by forcing one of the two delayed inputs and composing it.

(defun decompose (left composed)
  (values left (funcall (inverse-lft left) composed)))

(defun decompose-range? (left composed if-success if-failure)
  (multiple-value-bind (left right) (decompose left composed)
    (if (range? right)  ;; does it still narrow?
        (funcall if-success left right)
        (funcall if-failure))))

(defun try-decompose-digit (lft if-success if-failure)
  (decompose-range?
   (make-lft 2 1 0 1) lft
   if-success
   (lambda ()
     (decompose-range?
      (make-lft 1 0 1 2) lft
      if-success
      (lambda ()
        (decompose-range?
         (make-lft 3 1 1 3) lft
         if-success
         if-failure))))))
    
(defun binary-expression-decompose-digit (binary-expression)
  (try-decompose-digit
   (slot-value binary-expression 'bilft)
   (lambda (digit bilft)
     (values digit (make-instance 'binary-expression
                                  :generation (slot-value binary-expression 'generation)
                                  :bilft bilft
                                  :delayed-left (slot-value binary-expression 'delayed-left)
                                  :delayed-right (slot-value binary-expression 'delayed-right))))
   (lambda ()
     (binary-expression-decompose-digit (refine-disjoint binary-expression)))))

(defun binary-expression->lft-stream (binary-expression)
  (multiple-value-bind (digit remainder) (binary-expression-decompose-digit binary-expression)
    (cons-lft-stream digit 
                     (binary-expression->lft-stream remainder))))

Now we have the machinery we need to do arithmetic on pairs of LFT streams.

Infinite Expression Trees

You can only go so far with a finite number of arithmetic operations, but you can go much further with an infinite number. For example, you can create converging infinite series. We can recursively generate an infinite tree of binary expressions. We unfold the tree and call a generating function at each recursive level.

(defun unfold-expression-tree-1 (left generate counter)
  (funcall (funcall generate counter)
           left
           (delay-lft-stream (unfold-expression-tree-1 left generate (+ counter 1)))))

Often the root of the tree is special, so the usual way we call this is

(defun unfold-expression-tree (root-bilft left generate)
  (funcall root-bilft
           left
           (delay-lft-stream (unfold-expression-tree-1 left generate 1))))

Square Root of Infinite Composition

To compute the square root of an infinite composition, we create an infinite tree. Each level of the tree refines the estimate of the square root of the estimate from the next level down in the tree.

(defun sqrt-lft-stream (lft-stream)
  (unfold-expression-tree-1
    lft-stream
    (constantly (make-bilft 1 2 1 0
                            0 1 2 1))
    0))

but we need not construct an infinite tree if each level is identical. We just re-use the level:

(defun sqrt-lft-stream (lft-stream)
  (let ((stream nil))
    (setq stream
          (funcall (make-bilft 1 2 1 0
                               0 1 2 1)
                   lft-stream
                   (delay-lft-stream stream)))
    stream))

This feeds back the square root into its own computation. (These ideas are from Peter Potts and Reinhold Heckmann.)

e^x and log(x)

Peter Potts gives these formulas for e^x and log(x) where x is an infinite composition of LFTs:

(defun %exp-lft-stream (lft-stream)
  (unfold-expression-tree-1
   (lambda (n)
       (make-bilft (+ (* 2 n) 2) (+ (* 2 n) 1) (* 2 n) (+ (* 2 n) 1)
                   (+ (* 2 n) 1) (* 2 n) (+ (* 2 n) 1) (+ (* 2 n) 2)))
   0
   (funcall (inverse-lft (make-lft 1 -1 1 1)) lft-stream)))

(defun %log-lft-stream (lft-stream)
  (unfold-expression-tree
   (make-bilft 1 1 -1 -1
               0 1  1  0)
   (lambda (n)
     (make-bilft n (+ (* n 2) 1)       (+ n 1) 0
                 0       (+ n 1) (+ (* n 2) 1) n))
   lft-stream))

These infinite expression trees converge only on a limited range, so we need to use identities on the arguments and results to extend the range to cover a wide set of inputs.