Find the best rational approximation to π
Setup: The problem is to find the best ratio of two integers
that approximates π and output the numerator and denominator of
that ratio. We will be restricted to numbers less than 10000. One
way to solve this problem is to build a 10000x10000 matrix where
the i,jth number is i/j-π and then look for the
smallest number in the matrix - if that number is at row i
and column j then the numerator is i and the
denominator is j. This method is so inefficient that matlab
cannot solve it in a reasonable amount of time. The reason is many
many more numbers than needed are considered. For example, in
row i, if for some j, i/j-π is positive,
then smaller values for j will result in even a more positive
difference and they do not need to be considered. Also,
if i/j-π is negative, larger values of j need not
be considered. This leads to this approach for
finding the best ratio.
Express what a computer is doing in English: The computer
performs the following operations numerous times:
- Looks at the next numerator or denominator in sequence. Say it calls
them numer and denom.
- Compares the absolute value of the difference of
numer/denom and π with the value of a variable
called best_diff.
- Replaces best_diff's value with that of the absolute
value of the new difference if it is less than the value of best_diff.
- Saves numer and denom to best_numer
and best_denom when best_diff is replaced.
The computer also initializes some variables. Namely, the computer
sets initial values for numer, denom, best_numer,
best_denom, best_diff.
A description of an algorithm for approximating pi is as follows:
- Initialize all variables: numer and best_numer
get 3, denom and best_denom get 1, best_diff gets
3-π.
- Repeat the following while numer < 0
- If abs(numer/denom-pi) < best_diff then
- Set best_diff to abs(numer/denom-pi)
- Set best_numer and best_denom
- If the difference is positive, increment denom
Otherwise increment numer
Translation to Matlab code: by line number
- numer = 3;
denom = 1;
best_numer = 3;
best_denom = 1;
best_diff = 3-pi;
- while numer < 10000
- if abs(numer/denom - pi) < abs(best_diff)
- best_diff = numer/denom - pi;
- best_numer = numer;
best_denom = denom;
- if numer/denom-pi > 0 denom = denom + 1;
else numer = numer + 1; end
End result:
numer = 3;
denom = 1;
best_diff = 3 - pi;
best_numer = 3;
best_denom = 1;
while numer < 10000
if abs(numer/denom - pi) < abs(best_diff)
best_diff = numer/denom - pi;
best_numer = numer;
best_denom = denom;
fprintf('%d/%d (%12.10f)\n',numer,denom,best_diff);
fflush(stdout);
end
if numer/denom-pi > 0
denom = denom + 1;
else
numer = numer + 1;
end
end
fprintf('\nAnswer: %d/%d with %12.10f error\n',best_numer,best_denom,best_diff);