yap-6.3/Logtalk/examples/threads/functions/functions.lgt

:- protocol(find_rootp).

	:- info([
		version is 1.1,
		author is 'Paulo Moura and Paulo Nunes',
		date is 2006/11/26,
		comment is 'Default protocol for root find algorithms.']).

	:- public(find_root/5).
	:- mode(find_root(+object_identifier, +float, +float, +float, -float), one).
	:- info(find_root/5, [
		comment is 'Find the root of a function in the interval [A, B] given a maximum aproximation error.',
		argnames is ['Function', 'A', 'B', 'Error', 'Zero']]).

:- end_protocol.



:- protocol(functionp).

	:- info([
		version is 1.1,
		author is 'Paulo Moura and Paulo Nunes',
		date is 2006/11/26,
		comment is 'Default protocol for real functions of a single real variable.']).

	:- public(eval/2).
	:- mode(eval(+float, -float), one).
	:- info(eval/2, [
		comment is 'Calculates the function value.',
		argnames is ['X', 'Fx']]).

	:- public(evald/2).
	:- mode(evald(+float, -float), one).
	:- info(evald/2, [
		comment is 'Calculates the value of the function derivative.',
		argnames is ['X', 'DFx']]).

:- end_protocol.



:- object(f1,
	implements(functionp)).

	% x^2 - 4
	% 2.0 

	eval(X, Y) :-
		Y is X * X - 4.

	evald(X, Y) :-
		Y is 2 * X.

:- end_object.



:- object(f2,
	implements(functionp)).

	% x^7 + 9x^5 - 13x - 17
	% 1.29999999999945448
 
	eval(X, Y) :-
		Y is X**7 + 9*X**5 - 13*X - 17. 

	evald(X, Y) :-
		Y is 7*X**6 + 45*X**4 - 13. 

:- end_object.



:- object(f3,
	implements(functionp)).

	% (x - sqrt(2))^7
	% 1.41421356237309537
 
	eval(X, Y) :-
		Y is (X - sqrt(2.0))**8.

	evald(X, Y) :-
		Y is 8*(X - sqrt(2.0))**7. 

:- end_object.



:- object(f4,
	implements(functionp)).

	% x + x^2*sin(2.0/x)
	% 0.0
 
	eval(X, Y) :-
		Y is X + (X**2)*sin(2.0/X).

	evald(X, Y) :-
		Y is 1 + 2*X*sin(2.0/X) - 2*cos(2.0/X). 

:- end_object.



:- object(bisection,
	implements(find_rootp)).

	:- info([
		version is 1.2,
		author is 'Paulo Moura and Paulo Nunes',
		date is 2007/7/7,
		comment is 'Bisection algorithm.']).

	find_root(Function, A, B, Error, Zero) :-
		Function::eval(A, Fa),
		Function::eval(B, Fb),
		(	Fa > 0.0, Fb < 0.0 ->
			true
		;	Fa < 0.0, Fb > 0.0
		),
		X0 is (A + B) / 2.0,
		Function::eval(X0, F0),
		bisection(Function, A, B, X0, F0, Error, Zero).

	bisection(_, _, _, Xn, Fn, Error, Xn) :-
		abs(Fn) < Error,
		!.
	bisection(Function, An, Bn, _, _, Error, Zero) :-
		Xn1 is (An + Bn) / 2.0,
		Function::eval(Xn1, Fn1),
		Function::eval(An, FAn),
		(	Fn1*FAn < 0.0 ->
			An1 is An,
			Bn1 is Xn1
		;	An1 is Xn1,
			Bn1 is Bn
		),
		bisection(Function, An1, Bn1, Xn1, Fn1, Error, Zero).

:- end_object.



:- object(newton,
	implements(find_rootp)).

	:- info([
		version is 1.2,
		author is 'Paul Crocker... No More Coffee',
		date is 2007/07/06,
		comment is 'Newton algorithm.']).

    find_root(Function, Xa, Xb, Deviation, Zero) :-
		Ac is (Xb - Xa) / 2,		
		newton(Function, Xa, Ac, Deviation, Zero).

	newton(_, Zero, Ac, Deviation, Zero) :-
		abs(Ac) < Deviation,
		!.
    newton(Function, X0, Ac, Deviation, Zero):-
        Xn1 is X0 + Ac,
        Function::eval(Xn1, Fx),
        Function::evald(Xn1, DFx),
        Ac2 is -(Fx/DFx),
        newton(Function, Xn1, Ac2, Deviation, Zero).	

:- end_object.


:- object(muller,
	implements(find_rootp)).

	:- info([
		version is 1.2,
		author is 'Paulo Moura and Paulo Nunes',
		date is 2006/11/26,
		comment is 'Muller algorithm.']).

	find_root(Function, Xa, Xb, Deviation, Zero) :-
		Xc is (Xa + Xb) / 2.0, 
		muller(Function, Xa, Xc, Xb, Deviation, Zero).

	muller(Function, Xa, Xb, Xc, Deviation, Zero) :-
		Function::eval(Xa, Ya),
		Function::eval(Xb, Yb),
		Function::eval(Xc, Yc),
		H1 is Xb - Xa,
		DDba is (Yb - Ya) / H1,
		Ac is Deviation + 1.0,
		muller(Function, Xa, Xb, Xc, Deviation, Ya, Yb, Yc, Ac, H1, DDba, Zero).

	muller(_, _, _, Xc, Deviation, _, _, _, Ac, _, _, Xc) :-
		abs(Ac) < Deviation,  
		!.
	muller(Function, Xa, Xb, Xc, Deviation, _, Yb, Yc, _, _, DDba, Zero) :-
		H2n is Xc - Xb,
		DDcbn is (Yc - Yb) / H2n,
		Cn is (DDcbn - DDba) / (Xc - Xa),
		Bn is DDcbn + H2n * Cn,
		Rn is Bn * Bn - 4.0 * Yc * Cn,
		(	Rn < 0.0 ->
			fail
		;	V is sqrt(Rn)
		),
		(	Bn > 0.0 ->
			Dn is Bn + V
		;	Dn is Bn - V
		),
		Acn is -(2 * Yc / Dn),
		Xan is Xb,
		Xbn is Xc,
		Xcn is Xc + Acn,
		Yan is Yb,
		Ybn is Yc,
		Function::eval(Xcn, Ycn),
		H1n is H2n,
		DDban is DDcbn,
		muller(Function, Xan, Xbn, Xcn, Deviation, Yan, Ybn, Ycn, Acn, H1n, DDban, Zero).

:- end_object.



:- object(function_root,
	implements(find_rootp)).

	:- info([
		version is 2.0,
		author is 'Paulo Moura and Paulo Nunes',
		date is 2007/07/05,
		comment is 'Multi-threading interface to root finding algorithms.']).

	:- threaded.

	:- public(find_root/6).
	:- mode(find_root(+object_identifier, +float, +float, +float, -float, -object_identifier), one).
	:- info(find_root/6, [
		comment is 'Finds the root of a function in the interval [A, B] given a maximum aproximation error. Returns the method used.',
		argnames is ['Function', 'A', 'B', 'Error', 'Zero', 'Method']]).

	find_root(Function, A, B, Error, Zero, Algorithm) :-
		threaded((
				(bisection::find_root(Function, A, B, Error, Zero), Algorithm = bisection)
			;	(newton::find_root(Function, A, B, Error, Zero), Algorithm = newton)
			;	(muller::find_root(Function, A, B, Error, Zero), Algorithm = muller)
			)).

%	an alternative, possibly better definition would be to make the methods simply fail in case of error:
%
%	find_root(Function, A, B, Error, Zero, Algorithm) :-
%		threaded((
%				(catch(bisection::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = bisection)
%			;	(catch(newton::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = newton)
%			;	(catch(muller::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = muller)
%			)).

	find_root(Function, A, B, Error, Zero) :-
		find_root(Function, A, B, Error, Zero, _).

:- end_object.