Tuesday, September 1, 2015

Linear fractional transforms and permutations


So what does group theory have to do with linear fractional transforms? I'm glad you asked. The answer is pretty complicated, though.

There's no good place to start here, so let's just dive in. This function computes the cross-ratio of four numbers:

(define (cross-ratio A B C D)
  (/ (* (- a b) (- c d))
     (* (- a d) (- c b))))

1 ]=> (cross-ratio 1 9 5 13)

;Value: 4/3

There's a geometric interpretation of the cross ratio, but for the moment, just think of it as a function that takes any four numbers and produces a value.

1 ]=> (cross-ratio 3 2 1 9)

;Value: -4/3

1 ]=> (cross-ratio 5 1 8 4)

;Value: 16/7

The cross ratio is preserved under linear fractional transforms:

1 ]=> (define lft2 (make-linear-fractional-transform 3 -1 2 7))

;Value: lft2

1 ]=> (cross-ratio (lft2 3) (lft2 2) (lft2 1) (lft2 9))

;Value: -4/3

1 ]=> (cross-ratio (lft 3) (lft 2) (lft 1) (lft 9))

;Value: -4/3

1 ]=> (cross-ratio 3 2 1 9)

;Value: -4/3

The cross ratio is also preserved under some permutations of its arguments.

1 ]=> (cross-ratio 3 2 1 9)

;Value: -4/3

1 ]=> (cross-ratio 1 9 3 2)

;Value: -4/3

1 ]=> (cross-ratio 2 3 9 1)

;Value: -4/3
but not all
1 ]=> (cross-ratio 3 2 9 1)

;Value: 4/7

1 ]=> (cross-ratio 2 9 3 1)

;Value: 7/3

Now here's the interesting part. There's a linear fractional transform that will "undo" the permutation of the arguments to cross-ration.

1 ]=> ((make-linear-fractional-transform -1 1 0 1) 7/3)

;Value: -4/3

1 ]=> ((make-linear-fractional-transform 1 0 1 -1) 4/7)

;Value: -4/3
There are 24 permutations of the argument list to cross-ratio, and here are the equivalence classes:
((-4/3 (1 9 3 2) (2 3 9 1) (3 2 1 9) (9 1 2 3))
 (-3/4 (1 2 3 9) (2 1 9 3) (3 9 1 2) (9 3 2 1))
  (3/7 (1 2 9 3) (2 1 3 9) (3 9 2 1) (9 3 1 2))
  (4/7 (1 9 2 3) (2 3 1 9) (3 2 9 1) (9 1 3 2))
  (7/4 (1 3 2 9) (2 9 1 3) (3 1 9 2) (9 2 3 1))
  (7/3 (1 3 9 2) (2 9 3 1) (3 1 2 9) (9 2 1 3)))
And these are the linear fractional transforms that permute among these:
(#[linear-fractional-transform 89333 x]
#[linear-fractional-transform 89334 1/x]
#[linear-fractional-transform 89335 (1 - x)]
#[linear-fractional-transform 89403 1/(1 - x)]
#[linear-fractional-transform 89370 x/(x - 1)]
#[linear-fractional-transform 89393 (x - 1)/x])

This is cool. Instead of thinking of a linear fractional transform as a continuous function that traces out a hyperbola, or as an interval on the real number line, we can think of a linear fractional transform as a function that computes a permutation.


Here is a function called d3 that is defined on the symbols a, b, c, d, e, and f.

1 ]=> (d3 'a 'b)

;Value: d
With six symbols, there's only 36 possible outcomes, so we can enumerate them:
    a b c d e f
a ((e d f b a c) 
b  (f e d c b a) 
c  (d f e a c b) 
d  (c a b f d e) 
e  (a b c d e f) 
f  (b c a e f d))

Here's a hack. The identity element can be seen to be 'e, so list that row and column first:

    e a b c d f
e ((e a b c d f) 
a  (a e d f b c) 
b  (b f e d c a) 
c  (c d f e a b) 
d  (d c a b f e) 
f  (f b c a e d))

Now the table headers are exactly the same as the first entries in the respective rows or columns, so you can do without them.

((e a b c d f) 
 (a e d f b c) 
 (b f e d c a) 
 (c d f e a b) 
 (d c a b f e) 
 (f b c a e d))

Another weird thing is that you can permute any rows but the first, or permute any columns but the first, and still have an equivalent table.

Anyway, d3 is the "symmetry group" of a triangle. Here the operations are

e = leave it alone
f = rotate clockwise 120 degrees
d = rotate counter-clockwise 120 degrees
a = hold vertex A in place, and flip the triangle over swapping vertex B and C
b = hold vertex B in place, and flip the triangle over swapping vertex A and C
c = hold vertex C in place, and flip the triangle over swapping vertex A and B
Now consider this,
e = #[linear-fractional-transform 89333 x]
f = #[linear-fractional-transform 89393 (x - 1)/x]
d = #[linear-fractional-transform 89403 1/(1 - x)]
a = #[linear-fractional-transform 89370 x/(x - 1)]
b = #[linear-fractional-transform 89335 (1 - x)]
c = #[linear-fractional-transform 89334 1/x]

Just as d3 permutes a triangle, these linear fractional transforms permute among the equivalence classes of cross ratios.



Thursday, August 27, 2015

n-ary functions and argument lists

Suppose we have a binary operation Op2 = (lambda (left right) ...). If it is closed over some type, you can make expression trees.
(Op2 (Op2 <arg1> <arg2>)
     (Op2
          (Op2 <arg3> <arg4>)
          <arg5>))
If Op2 is associative as well, these are equal:
(Op2 (Op2 <arg1> <arg2>)
     (Op2 <arg3>
         (Op2 <arg4> <arg5>)))

(Op2 <arg1> 
     (Op2 (Op2 <arg2> <arg3>)
          (Op2 <arg4> <arg5>)))

This makes the tree structure irrelevant so long as the fringe of the tree stays in order, so we can flatten the tree by making an N-ary version of the binary operation:
(define (binary->nary Op2)
  (lambda (a1 a2 . an)
    (fold-left Op2 (Op2 a1 a2) an)))

((binary->n-ary Op2) <arg1> <arg2> <arg3> <arg4> <arg5> <arg6> etc.)
The value of this expression naturally depends upon the values of the arguments. Changing the argument list is highly likely to change the value of the entire expression. However, we can make certain changes to the argument list without changing the value of the entire expression. If we know what those changes are, we can manipulate the argument list before invoking the operator, and still get the same answer. Naturally, most of the changes we can make to the argument list depend on the specifics of Op2, but it turns out that some interesting changes are possible without knowing any specifics, only knowing a few high-level properties of Op2.

Obviously, the first thing we can do is reduce the argument list through evaluation. Simply replace the first two arguments with the value of (Op2 <arg1> <arg2>)
((binary->n-ary Op2) <arg1> <arg2> <arg3> <arg4> <arg5> <arg6> etc.) =

((binary->n-ary Op2) (Op2 <arg1> <arg2>) <arg3> <arg4> <arg5> <arg6> etc.) =

((binary->n-ary Op2) <result> <arg3> <arg4> <arg5> <arg6> etc.)
Since Op2 is associative, we can replace any 2 adjacent arguments with their combination.

Now suppose there is an identity element among the arguments we can give to Op2.
(Op2 <arg> id) = <arg>  and
(Op2 id <arg>) = <arg>
We can do this:
(define (binary->nary Op2)
  (lambda an
    (fold-left Op2 id an)))

(define Op (binary->nary Op2))
Which is cleaner than the original. We also get a new way to manipulate the argument list to Op. We can add the identity element anywhere we wish, or we can delete the identity element wherever we find one.
(Op <arg1> <arg2> <arg3> Id <arg4> <arg5> <arg6> etc.) =

(Op <arg1> <arg2> <arg3> <arg4> <arg5> <arg6> etc.) =

(Op <arg1> Id <arg2> <arg3> <arg4> <arg5> Id <arg6> etc.)

One more restriction. We want Op2 to be invertible. Suppose (Op2 <arg1> <arg2>) = <result>. Op2 is invertible if, given any two of <arg1>, <arg2>, and <result>, the third can be uniquely determined. If you have one <arg> and a <result>, you can run things backwards and get the other <arg>.

Requiring Op2 to be invertible has many consequences, some of them quite non-obvious. An obvious consequence, though, is that we can define inverse elements. If (Op2 <argA> <argB>) = Id, then we say that <argB> is the inverse of <argA> (and vice versa). We find the inverse of an argument by fixing the output as the identity element and running Op2 backwards to find the other argument.

This gives us the final way to manipulate the argument list. If an element appears next to its inverse, both can be removed:
(Op <arg1> <arg2> <arg3> (inverse <arg3>) <arg5> <arg6> etc.) =
(Op <arg1> <arg2> (Op2 <arg3> (inverse <arg3>)) <arg5> <arg6> etc.) =
(Op <arg1> <arg2> Id <arg5> <arg6> etc.) =
(Op <arg1> <arg2> <arg5> <arg6> etc.)

So here are all the restrictions on Op2:
  • Closed over an set of arguments
  • Associative
  • Has an identity argument
  • Invertible
If Op2 has these properties (and a lot of binary operations do), then we can define an n-ary Op and play with its argument list. If you do this, you might notice that it looks kind of familiar:
(op f g (inverse g) j id h) = (op f j id h) = (op f j h)

The argument list sort of looks like a function pipeline. The allowed manipulations of the argument list are compatible with a function pipeline, too. In fact, it could be a function pipeline if Op is the compose operator, and f, g, and h are appropriate invertible unary functions. But whatever it is, the point is that it looks enough like a function pipeline that we can pretend that it is.

Sunday, August 23, 2015

Playing with linear fractional transforms

I wanted to play with continued fractions and linear fractional transforms, so I wrote some code to make it easier. A linear fractional transform (also called a homographic function or Mobius transform) is a function like this:
In MIT/GNU Scheme:
;; x => (3x + 1)/(x + 4)

1 ]=> (define foo (make-linear-fractional-transform 3 1 1 4))

;Value: foo

1 ]=> (foo 2)

;Value: 7/6
I used an entity object so, in addition to invoking it on a number, there are two more ways to manipulate a linear fractional transform:
;; A predicate
1 ]=> (linear-fractional-transform? foo)

;Value: #t

;; And a CPS accessor
1 ]=> (lft/spread-coefficients foo (lambda (A B C D) (list A B C D)))

;Value 307: (3 1 1 4)
I also added a print method:
1 ]=> foo

;Value 308: #[linear-fractional-transform 308 (3x + 1)/(x + 4)]

As I mentioned in a prior post, you can partly apply a linear fractional transform:
1 ]=> foo

;Value 308: #[linear-fractional-transform 308 (3x + 1)/(x + 4)]

1 ]=> (lft/partly-apply foo 2)

;Value 315: #[linear-fractional-transform 315 (7x + 3)/(6x + 1)]
Since I want to reason about applying a linear fractional transform to an argument, I wrote an abstraction for that:
;; Apply LFT foo to continued fraction phi.
1 ]=> (make-lft-application foo phi)

;Value 311: #[lft-application 311 (3x + 1)/(x + 4) {1 ...}]
So now we can write a procedure that takes an application, peels off the first term in the continued fraction, feeds it to the linear fractional transform, and creates a new application:
(define (application/step lft-application)
  (let ((lft (application-function lft-application))
 (cf  (application-continued-fraction lft-application)))
    (make-lft-application
     (lft/partly-apply lft (head cf))
     (tail cf))))

1 ]=> (define appl (make-lft-application lft/identity sqrt-two))

;Value: appl

1 ]=> appl

;Value 317: #[lft-application 317 x {1 2 2 2 ...}]

1 ]=> (application/step appl)

;Value 318: #[lft-application 318 (x + 1)/x {2 2 2 ...}]

1 ]=> (application/step (application/step appl))

;Value 319: #[lft-application 319 (3x + 1)/(2x + 1) {2 2 ...}]

1 ]=> (application/step (application/step (application/step appl)))

;Value 320: #[lft-application 320 (7x + 3)/(5x + 2) {2 ...}]
All these lft-application objects should be (numerically) equal.

In an earlier post I showed how a linear fractional transform can be partly evaluated by determining the integer-part of the transform. The integer-part of an application is the integer-part of the application-function. You can get the fractional part by subtracting the integer-part.

A digression

If you apply a linear fractional transform to zero, it's obvious the answer is B/D. On the other hand, if you apply a transform to a sufficiently large x, you'll get as close as you want to A/C.

If the denominator of a linear fractional transform is zero for some value of x, there should be a vertical asymptote at that point. That's the pole of the transform. The pole is at (- D)/C. The pole will be at zero if D is zero. It will be at a negative number if D and C are the same sign and at a positive number if D and C differ in sign.

If you take a linear fractional transform with a pole at a negative number, and you sweep the input from 0 up on to infinity, the output will vary smoothly and monotonically from B/D approaching A/C and staying between both values at all times.
1 ]=> lft1

;Value 675: #[linear-fractional-transform 675 (3x + 1)/(4x + 2)]

1 ]=> (lft1 0)

;Value: 1/2

1 ]=> (lft1 1000000)

;Value: 3000001/4000002

1 ]=> (exact->inexact 3000001/4000002)

;Value: .7499998750000625

(On the other hand, if the pole is at a positive number, as you sweep the input from 0 up to infinity, the output starts at B/D, but flees away from A/C until the input gets to the pole. Then the output approaches A/C, but from the opposite direction. In any case, if the pole is positive, then the output will vary from B/D and eventually approach A/C, but never being between them.)

We can represent intervals as linear fractional transforms. The endpoints of the interval are A/C and B/D.

To get the width of the interval, just subtract the endpoints: A/C - B/D = (A*D - B*C)/(C * D)


Imagine you are performing some calculation with continued fractions. Since there may be an infinite number of terms, the calculation will proceed incrementally, using up terms as needed and generating other terms. So you can think of a more complex calculation as a tree, where a node in the tree is a linear fractional transform and the continued fraction terms flow between the nodes.

When we do an application/step, we move a term from the continued fraction into the linear fractional transform. Now consider a term as an element of information. We've moved this information out of the continued fraction and into the linear fractional transform. The information is apparently "stored" in the linear fractional transform until it is extracted as an output term for the next stage in the computation. But if you think about it, the "format" of the information is different depending upon whether it is flowing between nodes, where it is a series of continued fraction terms, or if it is stored in a linear fractional transform, where it is encoded somehow in the values of the coefficients. The act of partly-evaluating a linear fractional transform is in effect "encoding" some information as a continued fraction term. Partly applying a linear fractional transform is in effect "decoding" the continued fraction term (presumably generated by an earlier computation). Why not change to a more efficient communication channel?

When a node sends information to another node, instead of converting the information to several continued fraction terms to be assembled at the other end, we'll send the information as a single linear fractional transform. It contains the desired information already in the right "format". (See Peter Potts's work.)

A digression


What happens if we compose two linear fractional transforms?
(compose (lambda (x)
           (/ (+ (* A x) B)
              (+ (* C x) D)))
         (lambda (y)
           (/ (+ (* p y) q)
              (+ (* r y) s))))
We get
(lambda (x)
   (/ (+ (* A (/ (+ (* p x) q)
                 (+ (* r x) s))) B)
      (+ (* C (/ (+ (* p x) q)
                 (+ (* r x) s))) D)))
Which, after some algebra, turns into this:
(lambda (x)
   (/ (+ (* (+ (* A p) (* B r)) x) (+ (* A q) (* B s)))
      (+ (* (+ (* C p) (* D r)) x) (+ (* C q) (* D s)))))
Which is equivalent to this:
(lambda (x)
  (let ((E (+ (* A p) (* B r)))
        (F (+ (* A q) (* B s)))
 (G (+ (* C p) (* D r)))
 (H (+ (* C q) (* D s))))

   (/ (+ (* E x) F)
      (+ (* G x) H))))
Which you can see is another linear fractional transform.

If we have a linear fractional transform
(lambda (x)
  (/ (+ (* A x) B)
     (+ (* C x) D)))
It's inverse (if it has one) is:
(lambda (x)
  (/ (+ (* D x) (- B))
     (+ (* (- C) x) A))))
Which is yet another linear fractional transform. These things are everywhere.

Let's see, if we have a binary operation binop that is
  1. Closed over some set, i.e. given any two elements of the set, the operation applied to the elements produces another element of the set. In other words, binop takes two arguments, returns one value, and the type of both arguments and return value are the same.
  2. Associative, i.e. (binop a (binop b c)) = (binop (binop a b) c)
  3. Has an identity argument. A "left identity" is an argument such that (binop left-identity x) = x. A "right identity" is an argument such that (binop x right-identity) = x. An "identity" argument works as a left or right identity.
  4. Is invertible, i.e. for any objects a and b, there is a unique object x such that (binop a x) = b and a unique object y such that (binop y b) = a

then we have a group.

The compose function is a binary operation. When you compose a linear fractional transform with another, you get a third linear fractional transform.
1 ]=> (define lft1 (make-linear-fractional-transform 3 1 4 2))

;Value: lft1

1 ]=> (define lft2 (make-linear-fractional-transform 5 1 1 0))

;Value: lft2

1 ]=> (lft/compose lft1 lft2)

;Value 662: #[linear-fractional-transform 662 (16x + 3)/(22x + 4)]
Linear fractional transforms are associative.
1 ]=> (define lft3 (make-linear-fractional-transform 7 2 1 3))

;Value: lft3

1 ]=> (lft/compose lft1 (lft/compose lft2 lft3))

;Value 663: #[linear-fractional-transform 663 (115x + 41)/(158x + 56)]

1 ]=> (lft/compose (lft/compose lft1 lft2) lft3)

;Value 664: #[linear-fractional-transform 664 (115x + 41)/(158x + 56)]

The linear fractional transform (make-linear-fractional-transform 1 0 0 1) is both a left and right identity when composed with another linear fractional transform.
1 ]=> (define lft/identity (make-linear-fractional-transform 1 0 0 1))

;Value: lft/identity

1 ]=> (lft/compose lft/identity lft1)

;Value 665: #[linear-fractional-transform 665 (3x + 1)/(4x + 2)]

1 ]=> (lft/compose lft1 lft/identity)

;Value 666: #[linear-fractional-transform 666 (3x + 1)/(4x + 2)]
Given lft1 and lft2, there is a unique linear fractional transform x such that (compose lft1 x) = lft2, and a unique linear fractional transform y such that (compose y lft1) = lft2. x should be (compose (inverse lft1) lft2), and y should be (compose lft2 (inverse lft1))
1 ]=> lft1

;Value 675: #[linear-fractional-transform 675 (3x + 1)/(4x + 2)]

1 ]=> lft2

;Value 687: #[linear-fractional-transform 687 (5x + 1)/x]

1 ]=> (define x (lft/compose (lft/inverse lft1) lft2)))

;Value: x

1 ]=> (lft/compose lft1 x)

;Value 690: #[linear-fractional-transform 690 (5x + 1)/x]

1 ]=> (define y (lft/compose lft2 (lft/inverse lft1)))

;Value: y

1 ]=> (lft/compose y lft1)

;Value 691: #[linear-fractional-transform 691 (5x + 1)/x]
It sure looks like linear fractional transforms form a group under function composition.
I guess it's time to learn a little group theory.

Monday, September 22, 2014

A couple more homographic function tricks

A generalized continued fraction is an expression of the form:
You can partly apply a homographic function to a generalized continued fraction if you have a stream of the ai and bi
(define (partly-apply-general hf nums dens)
  (let ((a (first  hf))
        (b (second hf))
        (c (third  hf))
        (d (fourth hf)))
    (if (empty-stream? nums)
        (values (list a a
                      c c)
                nums
                dens)
        (let ((num (head nums))
              (den (head dens)))
          (values (list (+ (* a den) (* b num)) a
                        (+ (* c den) (* d num)) c)
                  (tail nums)
                  (tail dens))))))

(define (print-hf-general hf nums dens)
  (call-with-values (lambda () (partly-evaluate hf))
    (lambda (term hf*)
      (if (not term)
          (call-with-values (lambda () (partly-apply-general hf nums dens))
            print-hf-general)
          (begin
            (display term) 
            ;; take reciprocal and multiply by 10
            (let ((a (first hf*))
                  (b (second hf*))
                  (c (third hf*))
                  (d (fourth hf*)))
              (print-hf-general (list (* c 10) (* d 10)
                                      a        b)
                           nums dens)))))))
For example, we can compute pi from this generalized continued fraction:
(print-hf-general (list 0 4 1 0)
              ;; [1 1 4 9 16 25 ...]
       (cons-stream 1
      (let squares ((i 1))
    (cons-stream (* i i) (squares (1+ i)))))
              ;; [1 3 5 7 9 11 ...]
       (let odds ((j 1)) 
     (cons-stream j (odds (+ j 2)))))

314159265358979323846264338327950^G
; Quit!
A bi-homographic function is a function of the form:
(define (bi-homographic a b c d e f g h)
  (lambda (x y)
    (/ (+ (* a x y) (* b x) (* c y) d)
       (+ (* e x y) (* f x) (* g y) h))))
Like a homographic function, you can partly evaluate a bi-homographic function and generate a continued fraction. You can also partly apply a bi-homographic function to a pair of continued fractions. When you do this, you have a choice of which continued fraction to be the object of the partial application. There's about twice as much nasty math involved, but the gist is that a bi-homographic function takes two continued fractions as arguments and produces one continued fraction as a result.

It turns out that addition, subtraction, multiplication and division are bi-homographic functions, so one can incrementally compute sums and products of continued fractions.

Thursday, September 18, 2014

A useful, if somewhat pointless, trick with homographic functions

In my previous posts I showed that if you are applying a homographic function to a continued fraction, you can partly evaluate the answer before you completely apply the function. Instead of representing homographic functions as lambda expressions, I'll represent them as a list of the coefficients a, b, c, and d in (lambda (t) (/ (+ (* a t) b) (+ (* c t) d))). I'll represent a simple continued fraction as a stream of the integer terms in the denominators.
Here is how you partly apply a homographic function to a continued fraction:
(define (partly-apply hf cf)
  (let ((a (first  hf))
        (b (second hf))
        (c (third  hf))
        (d (fourth hf)))
    (if (empty-stream? cf)
        (values (list a a
                      c c)
                cf)
        (let ((term (head cf)))
          (values (list (+ (* a term) b) a
                        (+ (* c term) d) c)
                  (tail cf))))))
Partly evaluating a homographic function involves looking at the limits of the function as t starts at 1 and goes to infinity:
(define (partly-evaluate hf)
  (let ((a (first hf))
        (b (second hf))
        (c (third hf))
        (d (fourth hf)))

    (if (and (same-sign? c (+ c d))
             (let ((limit1 (quotient      a       c))
                   (limit2 (quotient (+ a b) (+ c d))))
               (= limit1 limit2)))
        (let ((term (quotient a c)))
          (let ((new-c (- a (* c term)))
                (new-d (- b (* d term))))
            (values term (list c d new-c new-d))))
        (values #f #f))))
We can combine these two steps and make something useful. For example, we can print the value of applying a homographic function to a continued fraction incrementally, printing the most significant digits before computing further digits.
(define (print-hf-cf hf cf)
  (call-with-values (lambda () (partly-evaluate hf))
    (lambda (term hf*)
      (if (not term)
          (call-with-values (lambda () (partly-apply hf cf))
            print-hf-cf)
          (begin
            (display term) 
            ;; take reciprocal and multiply by 10
            (let ((a (first hf*))
                  (b (second hf*))
                  (c (third hf*))
                  (d (fourth hf*)))
              (print-hf-cf (list (* c 10) (* d 10)
                                 a        b)
                           cf)))))))
But how often are you going to apply a homographic function to a continued fraction? Fortunately, the identity function is homographic (coefficients are 1 0 0 1), so applying it to a continued fraction doesn't change the value. The square root of 2 is a simple continued fraction with coefficients [1 2 2 2 ...] where the 2s repeat forever. We apply the identity homographic function to it and print the results:
(printcf (list 1 0 0 1) sqrt-two)
14142135623730950488016887242096980785696718^G
; Quit!
As you can see, we start printing the square root of two and we don't stop printing digits until the user interrupts.

A fancier version could truncate the output and print a decimal point after the first iteration.

Friday, September 5, 2014

Another stupid homographic function trick

In my last post I showed that if you take a homographic function and apply it to a fraction, you can partly apply the function to the integer part of the fraction and get a new homographic function. The new function can be applied to the non-integer part of the fraction to generate an answer equivalent to the original function applied to the original fraction.

It turns out that you can go in the other direction as well. You can partly evaluate a homographic function. For example, consider this homographic function:
((lambda (t)
   (/ (+ (* 70 t) 29)
      (+ (* 12 t)  5))) n)
Which we intend to apply to some positive number n. Even if all we know is that n is positive, we can deduce that the value of the homographic function is between 29/5 (when t is 0) and 70/12 (as t goes to infinity). The integer part of those values are the same, so we can factor that out:
(+ 5 (/ 1 ((lambda (t)
               (/ (+ (* 12 t) 5)
                  (+ (* 10 t) 4))) n)))
The partial answer has an integer value of 5 and a fractional part that contains a new homographic function applied to our original n. We can do it again:
(+ 5 (/ 1
        (+ 1 (/ 1
                ((lambda (t)
                   (/ (+ (* 10 t) 4)
                      (+ (* 2 t) 1))) n)))))
The fractional part of the answer can itself be factored into another integer and a new homographic function applied to our original n.

Wednesday, September 3, 2014

Stupid homographic function trick

A function of the form
is called a homographic function.  Here is one in Scheme:
(lambda (t)
   (/ (+ (* 3 t) 4)
      (+ (* 5 t) 2)))
And here is what it's graph looks like:
If you multiply all the coefficients (a, b, c, and d) by the same number, it doesn't change the function. For instance, this homographic function:
(lambda (t)
   (/ (+ (* 21 t) 28)
      (+ (* 35 t) 14)))
is the same as the one above. If one of your coefficients isn't an integer, don't worry, you can multiply everything by the denominator and get an equivalent homographic function. On the other hand, you can divide all your coefficients by their greatest common divisor and get an equivalent homographic function with smaller coefficients. We'll keep our homographic functions in smallest integer form.

A rational number can be written as the sum of an integer and a fraction less than one. For example, 23/5 = 4 + 3/5.

Let's apply a homographic function to a rational number:
((lambda (t)
   (/ (+ (* a t) b)
      (+ (* c t) d))) (+ x y/z))

;; substitute
(/ (+ (* a (+ x y/z)) b)
   (+ (* c (+ x y/z)) d))

;; distribute the multiplication
(/ (+ (* a x) (* a y/z) b)
   (+ (* c x) (* c y/z) d))

;; multiply top and bottom by z/y
(/ (* z/y (+ (* a x) (* a y/z) b))
   (* z/y (+ (* c x) (* c y/z) d)))

;; distribute the multiplication
(/ (+ (* a x z/y) (* a y/z z/y) (* b z/y))
   (+ (* c x z/y) (* c y/z z/y) (* d z/y)))

;; simplify
(/ (+ (* a x z/y) a (* b z/y))
   (+ (* c x z/y) c (* d z/y)))

;; rearrange terms
(/ (+ (* a x z/y) (* b z/y) a)
   (+ (* c x z/y) (* d z/y) c))

;; factor out z/y
(/ (+ (* (+ (* a x) b) z/y) a)
   (+ (* (+ (* c x) d) z/y) c))

Now we do something tricky. We abstract out the z/y term:
((lambda (t)
   (/ (+ (* (+ (* a x) b) t) a)
      (+ (* (+ (* c x) d) t) c))) (/ z y))
If introduce a let form, we can see something interesting:
((lambda (t)
   (let ((a1 (+ (* a x) b)) 
         (b1 a)
         (c1 (+ (* c x) d))
         (d1 c))
     (/ (+ (* a1 t) b1)
        (+ (* c1 t) d1)))) (/ z y))
We find a new homographic function being applied to a new rational number. The new homographic function has coefficients related to the original one, and the new rational number is the reciprocal of the fractional part of the original rational number. So if we have a homographic function hf applied to a fraction of the form x + y/z, we can easily find a new homographic function hf' that when applied to z/y will produce the same answer as the original. Applying a homographic function to a fraction has the effect of "eating" the integer part of the fraction, which generates a new homographic function that is applied to the reciprocal of the fractional part.