the solution to Project Euler 100

computer_discs.m

function [a,b]=computer_discs()

b1=21;
a1=15;
b2=85+35;
a2=85;

[t1,s1]=computer_st_square(a1,b1);
[t2,s2]=computer_st_square(a2,b2);

%computer the coefficent of the recursive formula
A=[s1,2*t1;t1,s1];
y=[s2;t2];
x=inv(A)*y;
x=round(x);

%the discs number int total
discs=10^12;
A=[s2,2*t2;t2,s2];
while 1,
    %computer the next a,b
    st=A*x;
    s=st(1);
    t=st(2);
    [a,b]=computer_ab(s,t);
    fprintf('a is %d, b is %d\n',a,b);
    if b>discs,
        fprintf('The answer is :a is %d, b is %d, pr is %f\n',a,b,a*(a-1)/(b*(b-1)));
        return;
    else
        A=[s,2*t;t,s];
    end
end
end

function [t,s]=computer_st(a,b)
t=sqrt(1-2*b+2*b^2);
s=sqrt(2*t^2-1);
end

function [a,b]=computer_ab(s,t)
b=(1+sqrt(2*t^2-1))/2;
a=(1+sqrt(1-2*b+2*b^2))/2;
end