echo off;

% Global constants
T = 10^(-5);
Y = 1.5;


% Quadratic Interpolation
h  = 1;
i  = 0;
a  = 0.5;
b  = 0.75;
fa = gamma(a) - Y;

disp(sprintf('\n\ni   x         f(x)'));
disp(sprintf('---------------------'));

while( abs(h) > T )
	fb = gamma(b);
	disp(sprintf('%d   %1.5f   %1.5f', i, b, fb));
    fb = fb - Y;
    h  = b - a;
    a  = b;
    b  = b - fb * h / (fb - fa);
	fa = fb;
    i  = i + 1;   
end


% Inverse Quadratic Interpolation
h  = 1;
i  = 0;
a  = 0.5;
b  = 0.75;
c  = 1.0;
fa = gamma(a) - Y;
fc = gamma(c) - Y;

disp(sprintf('\n\ni   x         f(x)'));
disp(sprintf('---------------------'));

while( abs(h) > T )
	fb = gamma(b);
    disp(sprintf('%d   %1.5f   %1.5f', i, b, fb));
    fb = fb - Y;
	u  = fb / fc;
	v  = fb / fa;
	w  = fa / fc;
    fc = fa;
	fa = fb;
	p  = v * (w * (u - w) * (c - b) - (1 - u) * (b - a));
	q  = (w - 1) * (u - 1) * (v - 1);
	h  = p / q;
    c  = a;
	a  = b;
	b  = b + h;
    i  = i + 1;
end


% Linear Fractional Interpolation
h  = 1;
i  = 0;
a  = 0.5;
b  = 0.75;
c  = 1.0;

disp(sprintf('\n\ni   x         f(x)'));
disp(sprintf('---------------------'));

while( abs(h) > T )
	fa = gamma(a) - Y;
    fb = gamma(b) - Y;
    fc = gamma(c);
    disp(sprintf('%d   %1.5f   %1.5f', i, c, fc));
    fc = fc - Y;
	p  = (a - c) * (b - c) * (fa - fb) * fc;
	q  = (a - c) * (fc - fb) * fa - (b - c) * (fc - fa) * fb;
	h  = p / q;
	c  = c + h;
    i  = i + 1;
end