NetNews Usenet Archive 1992 #23/NN_1992_23.iso/spool/sci/math/13321

Text File | 1992-10-16 | 5.0 KB | 112 lines

Newsgroups: sci.math
Path: sparky!uunet!morrow.stanford.edu!CSD-NewsHost.Stanford.EDU!Sunburn.Stanford.EDU!pratt
From: [e-mail] (Vaughan R. Pratt)
Subject: Re: Real -> fraction algorithm needed (also: FAREY)
Message-ID: <[e-mail]>
Sender: [e-mail]
Organization: Computer Science Department,  Stanford University.
References: <[e-mail]>
Date: Fri, 16 Oct 1992 20:00:45 GMT
Lines: 100

Boy, talk about serendipity.  First we get someone who wants to see
Farey sequences used somewhere:

In article <[e-mail]> [e-mail] (Naved Rehman) writes:
>are there any applications of FAREY SEQUENCE?

Then along comes a perfect customer for Farey sequences:

In article <[e-mail]> [e-mail] (Ted Hadley) writes:
>Greetings. I need an algorithm which I can code in 'C' which performs the
>following task. This, I feel, is the best place to ask. If this is a common
>request, please skip the details and cite a reference. Much appreciated. Here
>goes:
>
>Part 1: General case:
>
>Given an real number X, compute to integers A and B such that A/B most
>closely approximate X

So to make this post even more multithreaded let's suppose we want a
rational approximation to pi.  Start out with the approximations 0/1
and 1/0, guaranteed to be respectively lower and upper bounds on any
nonnegative real.  Add their numerators and denominators to get 1/1,
which being below pi becomes the new lower bound.  And so on until
happy.  Formally (i.e. in C :-)

  main(){int a=0,b=1,c=1,d=0,e,f;
  while (b < 1000) {
    e = a+c, f = b+d;
    if (e < f*3.14159265358979323846) a = e, b = f; else c = e, d = f;
    printf("%d/%d ", e, f);}}

which yields 1/1 2/1 3/1 4/1 7/2 10/3 13/4 16/5 19/6 22/7 25/8 47/15
69/22 91/29 113/36 135/43 157/50 179/57 201/64 223/71 245/78 267/85
289/92 311/99 333/106 355/113 688/219 1043/332 1398/445 1753/558
2108/671 2463/784 2818/897 3173/1010.

The last of these is about 3.141584, good to 3 parts in a million (see
the bit near the end of this post about expected errors and q^2).

>Part 2: First constraint:
>
>Repeat (1) such that the error, defined as the fractional part of A/(BX), is
>within a certain bound (say +/- E) AND A and B are minimized (i.e., the
>smallest possible values which produce a quotient with acceptible error).

Evidently e and f both grow monotonically, so by stopping as soon as we
have a value lying within [pi-E,pi+E] we have automatically met your
second condition that A and B are minimized (sorry I renamed them e and
f).

>Part 3: Second contraint: (this is the preferred form)
>
>Given X, E, and Amax and Bmax, where Amax and Bmax are the maximum allowable
>integers for A and B, respectively, compute all values of A and B for which
>the error E is acceptible. For example: |A| < Amax and |B| < Bmax and some
>value for E, list all possible {A, B} for which |fract(A/(BX))| < E.

What is really neat about the Farey sequence is that it generates every
rational exactly once, in its reduced form, i.e. the complete Farey
sequence is just the set Q of nonegative rationals, together with
infinity, standardly ordered.  So to solve part 3, compute that part of
the complete Farey sequence consisting of those rationals with
numerators and denominators bounded as above and lying between the
nearest rationals on each side of but not in [pi-E,pi+E].  (These two
"outside" rationals are repeatedly replaced by progressively better
approximations to pi-E below and pi+E above.) This gets you exactly
those rationals that you ordered, namely those lying within your
interval---the ones just outside the interval were needed only for the
process, not the answer.

In thinking about all this, bear in mind that *any* consecutive pair
a/b,c/d in the Farey sequence, however formed, satisfies

        |a  c|
        |    |  =  -1               (*)
        |b  d|

This is certainly true for the initial pair 0/1,1/0, and follows
whenever e/f is interpolated since we are now looking at the
determinants formed by adding the left column to the right and the
right column to the left respectively.  So you can interpolate new
entries in the Farey sequence in any order you like and nothing breaks
and moreover no rational can ever become "locked out" by adding
interpolants in an unfortunate order.

It also follows from (*) that (a,b) = (c,d) = 1, since otherwise the
common factor would divide the determinant.  Moreover c/d-a/b =
(bc-ad)/bd = 1/bd, telling you that consecutive finite elements of the
Farey sequence up to a given bound q on the denominators are spaced at
least 1/q^2 apart and at most 1/q apart (for the latter think about
numbers like .000001, .499999, and .500001), giving you an idea of your
prospects for rationally approximate a given real given such a bound q
on the denominator.

For a *really* practical application of Farey sequences see the bottom
of p.156, right column, of my paper on fast rendering of conic splines
in Siggraph'85, which uses some of the above ideas.
--
Vaughan Pratt   [e-mail]     exp(n pi sqrt(163)/3) for n=1..30