| Here's a BASIC program to investigate your problem. You can change the
error calculation at line 10 to get the results you are interested in.
The program prints a line each time it finds a better approximation than
the last. It does between 2M and 3M iterations, instead of M*M needed by an
exhaustive search. Enjoy!
1 ! RATIONAL APPROXIMATION TO A REAL NUMBER
! AUTHOR: LYNN YARBROUGH
INPUT "REAL"; R
PRINT "P", "Q", "P/Q", "ERROR"
P = 0
Q = 1
MN = R
10 E = ABS(P/Q-R)
IF E<MN THEN
MN = E
PRINT P,Q,P/Q,E
END IF
IF P/Q < R THEN
P = P + 1
ELSE
Q = Q + 1
END IF
GOTO 10 IF P<256 AND Q<256
END
|
| re.0 Each of the two minimization problems is solved by one of two
closest BRAs (bounded rational approximations). Thus, simply
find the two closest, and determine which one minimizes the
expression -- easier said than done.
I've seen a method that generates a sequence of BRAs, in order
of increasing value. Perhaps this would help (if someone else
could find the reference).
The greatest absolute error over the representable range of BPAs
is 1/2, and occurs when A (the number being approximated) is X+1/2,
where X is an integer, and B/2 < X < B, where B is the bound. The
greatest relative error occurs near zero, or when A is roughly 1/(2B).
I assume 'better' answers to these questions are wanted (for example,
when 1/B <= A <= 1).
re.1 Is this guaranteed to find the closest BRA? Yes, as it turns out.
It's possible to find a good initial approximation, and then continue
with the algorithm until the best BRA is found. The denominator Q
can be initially chosen to be floor(B/2) (where B is the bound), and
the numerator P as floor(A*Q), since, if the BRA equals p'/q', with
q' < floor(B/2), the algorithm will also find the solution (2p')/(2q').
A better algorithm should be possible. MACSYMA, for example, computes
BRAs with a fairly large bound.
re.2 The 'pattern' may be due to the continued fraction representation of
the number. For example, the sequence:
0/1, ..., P/Q, (P+Q)/(P+2Q), ...
rapidly converges to a well-known constant.
In fact, 'truncating' the continued fraction expansion of the number,
and replacing the truncated terms with either 0 or 1 may give the BRA.
|
| An interesting related theorem, known to Euler:
Let a be any irrational. Then there are an infinitly many rationals
p/q such that |a - p/q| < 1/q^2.
For a RATIONAL, on the other hand, it's easy to show that there are only
finitely many such rationals. (Obviously, we only consider the case (p,q) = 1.)
This distinction provides one of the fundamental techniques for proving numbers
irrational.
There is a sort of converse to this theorem, due to Liouville (generalizing a
result for quadratics known for many years):
Suppose a is the root of an nth degree polynomial with integer
coeficients, n>1. (Technical term: a is algebraic of order n
(or less, as written)). Then there exists a c = c(a) > 0 such
that |a - p/q| > c/q^n for all rationals p/q.
This turns out to be the basis of THE fundamental technique used for proving
numbers transcental - you find "too many" p/q's that approximate the alleged
transcental "too closely". A simple example:
Let z = Sigma 10^-n!. Let pj=10^j!*Sigma(1,j)10^-n!, and
qj = 10^j!. Then:
|z - pj/qj|=Sigma(j+1,infinity)10^n
< 10^-(j+1)!*(1+.1+.01+...)
= (10/9)qj^(-j-1)
< qj^-j
This is taken from "Transcendental Number Theory", by Alan Baker - a book I
do NOT recommend unless you like spending (literally) hours per page.
-- Jerry
|
| A topic in number theory that's directly related to your problem is that
of Farey sequences. A Farey sequence of order N on the interval [0,1] is
the set of rational numbers a/b, (a,b <= N) reduced to lowest terms and
sorted in ascending order. There are about 3*(N**2)/(PI**2) distinct
fractions for a given N... (for N = 255 there are 19821 fractions).
(The exact number is the SUM(phi(k)), k=1,...,N where phi(k) is Euler's
totient function: the number of integers < k, relatively prime to k.)
For a given real number x in [0,1] there will be a fraction a/b such that
ABS(x - a/b) <= 1/(b*(N+1)).
Here's a recursion that will generate all fractions Xi/Yi of order N:
X(0) = 0; Y(0) = 1; X(1) = 1; Y(1) = N;
LOOP: T = INT((Y(k)+N)/Y(k+1));
X(k+2) = T*X(k+1) - X(k);
Y(k+2) = T*Y(k+1) - Y(k);
k = k+1;
REPEAT;
The "mediant" of two fractions a/b and e/f is defined as (a+e)/(b+f) and lies
between a/b and e/f. In fact, each fraction in the series is the mediant
of its two neighbors.
Farey sequences are discussed in a recently published book that many readers
of this file will probably find interesting, "Number Theory in Science
and Communication", by M. R. Schroeder (Springer Verlag). Surprisingly
Knuth doesn't mention them at all in "The Art of Computer Programming".
- Jim
|
| The following program produces a plot of the Farey fractions Y/X of
order 240/M-1 (M selects a submultiple of the REGIS screen size,
2 is a good value) on a two dimensional lattice. This demonstrates
an interesting geometric interpretation of these fractions: they are
exactly those lattice points that would become 'visible' to an observer at
the origin as he rotates his look angle counterclockwise. As can be
seen, about 1/3 of all lattice points are visible, and the distribution
is quite uniform with only the slightst amount of lacunarity. To add some
additional information, points with a prime numerator or denominator are
plotted in green, others are plotted in red.
The two dimensional Fourier or Hadamard transform of this image also makes
an intriguing display.
- Jim
PROGRAM FPLOT
C
C Plot points in two dimensional lattice that are
C visible from the origin using Farey sequences
C
C Points with either X or Y a prime are plotted in green,
C Points with both X and Y composite are plotted in red.
C
C Loosely adapted from Stan's spiral prime number plotting program
C
IMPLICIT INTEGER*4 (A-Z)
C
DIMENSION PRIMES(100)
C
C Create a short list of prime numbers
C
PRIMES(1) = 2
NP = 1
DO 90 K=3,256,2
DO 80 J=1,NP
IF (MOD(K,PRIMES(J)).EQ.0) GOTO 90
80 CONTINUE
NP = NP+1
PRIMES(NP) = K
90 CONTINUE
C
M = IASK('Enter submultiple, M: ')
FLAG = IASK('Enter variant: (1 = primes, 2 = nonprimes, 3 = both')
CALL CLEAR
C
C Start the routine to calculate the successive Farey fractions
C
N = (240/M)-1
X0 = 1
Y0 = 0
X1 = N
Y1 = 1
300 CONTINUE
T = (X0+N)/X1
X2 = T*X1-X0
Y2 = T*Y1-Y0
X0 = X1
Y0 = Y1
X1 = X2
Y1 = Y2
DO 200 K=1,NP
IF (X0.EQ.PRIMES(K).OR.Y0.EQ.PRIMES(K)) GOTO 250
200 CONTINUE
C
C Plot in red if both numerator and denominator are composite
C
IF ((FLAG.AND.2).NE.0) CALL PLOT(X0,Y0,M,'R')
GOTO 270
C
C Plot in green if either numerator or denominator is prime
C
250 CONTINUE
IF ((FLAG.AND.1).NE.0) CALL PLOT(X0,Y0,M,'G')
270 CONTINUE
IF (X2.NE.Y2) GOTO 300
C
50 TYPE 100, 27
100 FORMAT(1X,A1,'\')
C
CALL EXIT
END
SUBROUTINE CLEAR
IMPLICIT INTEGER*4(A-Z)
ESC=27
TYPE 100, ESC,ESC,ESC
100 FORMAT(1X,A1,'[H',A1,'[2J',A1,'P1ps(M0(L0)(AL0))S(I0)S(E)S[0,0]')
RETURN
END
SUBROUTINE PLOT(X,Y,M,COLOR)
IMPLICIT INTEGER*4(A-Z)
PARAMETER (XOFFSET=350)
PARAMETER (YOFFSET=240)
CHARACTER*1 PREV_COLOR/'X'/
CHARACTER*1 COLOR
IF (COLOR.NE.PREV_COLOR) THEN
TYPE 101, COLOR
101 FORMAT('+W(I(',A1,$,'))')
PREV_COLOR=COLOR
END IF
TYPE 100, X*2*M,480-Y*2*M
100 FORMAT(' P[',I3,',',I3,']V[]')
RETURN
END
FUNCTION IASK(C)
C
C Ask user for an integer parameter
C
IMPLICIT INTEGER*4 (I-N)
CHARACTER *(*) C
L=LEN(C)
100 TYPE 10,C
10 FORMAT (' 'A<L>' '$)
READ (5,*,ERR=200) IV
IASK = IV
RETURN
200 TYPE 20,7
20 FORMAT ('+'A1,$)
GOTO 100
END
|