Error in an algorithm for continued fractions

Hello everyone,

I use Matlab R2020a to create an algorithm that gives all the digits of the continued fraction expansion of a real number (I know there is a function that does it more or less, but I haven't been able to use it).

function A=frac_continu(r,nmax)
    % r is a real number
    % nmax is the "limit rank"
    A=zeros(1,nmax);
    j=1;
    alpha=r;
    while j<= nmax && alpha-floor(alpha) ~= 0
        A(j)=floor(alpha);
        alpha=1/(alpha-A(j));
        j=j+1;
    end
end

What bothers me is that when I apply it to different real numbers, I don't get the right results from a certain rank.

For example, if I take the golden ratio (i.e. r = (1+sqrt(5))/2 and nmax=50), I should only get 1's, but the outgoing vector eventually gives me two 2's, etc.

So I wanted to know if the problem is algorithmic or if it's just due to a rounding problem or something else.

Thanks in advance!

[Steven Lord]

FYI you might find this post on Cleve's blog interesting.

[Kermatoni]

It's really interesting, thank you !

Answer 1

[Rik]

I suspect you are indeed running into rounding issues. Using an awful approach to unwind this continued fraction, I found this for nmax=30:

r = (1+sqrt(5))/2;
nmax=30;
fprintf('%.20f (true value)\n%.20f (approximate value)\n',r,unwind(frac_continu(r,30)))

1.61803398874989490253 (true value)
1.61803398875054060824 (approximate value)

function A=frac_continu(r,nmax)
    % r is a real number
    % nmax is the "limit rank"
    A=zeros(1,nmax);j=1;alpha=r;
    while j<= nmax && alpha-floor(alpha) ~= 0
        A(j)=floor(alpha);        alpha=1/(alpha-A(j));        j=j+1;
    end
end
function v=unwind(A)
%this is a TERRIBLE idea, don't ever use this, even a nested anonymous function is a better idea
if numel(A)==1,v=A;else,str=[sprintf('%d',A(1)),sprintf('+1/(%d',A(2:end)),repmat(')',1,numel(A)-1)];v=eval(str);end
end

After that many inversions, the rounding errors just keep on building up in your alpha, as you don't recalculate it from the A so far.

[Kermatoni]

Okay, thank you for your answer !