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:

  1. Looks at the next numerator or denominator in sequence. Say it calls them numer and denom.
  2. Compares the absolute value of the difference of numer/denom and π with the value of a variable called best_diff.
  3. 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.
  4. 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:

  1. Initialize all variables: numer and best_numer get 3, denom and best_denom get 1, best_diff gets 3-π.
  2. Repeat the following while numer < 0
  3.    If abs(numer/denom-pi) < best_diff then
  4.       Set best_diff to abs(numer/denom-pi)
  5.       Set best_numer and best_denom
  6.    If the difference is positive, increment denom
       Otherwise increment numer

 

Translation to Matlab code: by line number

  1. numer = 3;
    denom = 1;
    best_numer = 3;
    best_denom = 1;
    best_diff = 3-pi;
  2. while numer < 10000
  3. if abs(numer/denom - pi) < abs(best_diff)
  4. best_diff = numer/denom - pi;
  5. best_numer = numer;
    best_denom = denom;
  6. 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);