# Rational Approximation for e with Matlab

I've written a program that calculates the best rational approximation for e with a 1, 2, 3, 4, 5, and 6 digit denominator respectively (on Matlab). For a 5-digit denominator it takes about a minute. For a 6-digit denominator it takes 2 hours. I've tried to improve the code, but I am a beginner. Suggestions are kindly appreciated.

tic
clear
p=1;
q=1;
e=exp(1);
digits=1;
trial=1;
while numel(num2str(q))<=6
while numel(num2str(q))==digits
p=p+1;
if p/q>e
q=q+1;
end
A(trial)=p/q;
B(trial,1:2)=[p,q];
trial=trial+1;
end
for i=1:trial-1
error(i)=abs(A(i)-e);
end
[b,r]=min(error);
p_q(digits,1:2)=B(r,1:2);
digits=digits+1;
end
p_q
toc


This is the output:

p_q =
19           7
106          39
1264         465
25946        9545
271801       99990
1084483      398959

Elapsed time is 7232.153581 seconds.


Edit: Thanks to all who answered. Dennis and Jonas were right. It was the strings and growing matrices that slowed everything down. Besides I got two wrong values for p and q. The new code is this:

tic
p_best=1;
q_best=1;
e=exp(1);
for i=1:6
for q=10^(i-1):10^i-1
p=round(q*e);
if abs(p/q-e)<abs(p_best/q_best-e)
p_best=p;
q_best=q;
end
end
disp([p_best,q_best])
end
toc


The output is:

      19           7
193          71
1457         536
25946        9545
271801       99990
1084483      398959


Elapsed time is 0.054986 seconds.

It's faster and actually correct. I assumed I knew e like Dennis said. To Joni and hardmath: I can't use continued fractions because the teacher wanted us to use brute force, but thanks.

• I think this question may be better off at codereview or math. Two pointers: working with strings is usually quite slow, and it can be a lot faster if you think more about the properties of e. – Dennis Jaheruddin Feb 5 '13 at 16:39
• What really slows you down is that you grow all your arrays (A, B, etc) inside the loop. Preassign them to a sufficiently large array, and you'll get decent speed-ups – Jonas Feb 5 '13 at 17:13
• Indeed there's a nice characterization of the best rational approximations for real numbers in terms of the convergents of their continued fraction expansions. Obviously trial and error works, but if you are aiming for much faster answers as the digits allowed in the denominator increase, some theory can really speed things up. I'm pretty sure I remember Knuth covering this in vol. 2 of AofCP (Seminumerical Algorithms). – hardmath Feb 5 '13 at 18:19
• It's even more tempting to go the continued fraction route if you know that e (unlike pi) has a regular continued fraction expansion. – hardmath Feb 5 '13 at 18:24

Continued fractions are an easy way to calculate good rational approximations to real numbers. You can find a disussion of an algorithm to convert a floating point number into a fraction in my blog, with code in JavaScript which should be easy to convert to matlab:

function float2rat(x) {
var tolerance = 1.0E-6;
var h1=1; var h2=0;
var k1=0; var k2=1;
var b = x;
do {
var a = Math.floor(b);
var aux = h1; h1 = a*h1+h2; h2 = aux;
aux = k1; k1 = a*k1+k2; k2 = aux;
b = 1/(b-a);
} while (Math.abs(x-h1/k1) > x*tolerance);
return h1+"/"+k1;
}


If you're specifically interested in e you can take a shortcut and use the regularities in its continued fraction representation to produce rational approximations that are far more accurate than the machine epsilon. This method is described here.

In the evening I thought a bit about the question, and though it is a bit more radical than giving pointers on the code, I have created an example on how to approach this problem the matlab way.

Judging from your code the value of e can be considered known. With this I have created a small fragment to make over and under estimations and see which is the best one.

x=1e6:1e7-1; %To guarantee that the denominator has 6 digits
e = exp(1);
ylow = floor(x*e); %Candidates for over estimation
yhigh = ceil(x*e); %Candidates for under estimation
lowerr = e - ylow./x ;
higherr = yhigh./x - e;
lowfind = find(lowerr == min(lowerr));
highfind = find(higherr == min(higherr));
goodfractions = [ylow(lowfind),yhigh(highfind)
x(lowfind),x(highfind)]
errors = [lowerr(lowfind) higherr(highfind)]


Note that the runtime is quite fast but that the errors are or order 10^-14 which indicates that unless some measures are taken rounding errors may become relevant when you try to upscale this solution one or two digits.