/*
 **********************************************************************
 *
 *      CLP(R) Version 2.0	(Example Programs Release)
 *	(C) Copyright, March 1986, Monash University
 *
 **********************************************************************
 */

% use clpr -z 1e-20

/*
 * RKF45 implementation in CLP(R): See [Harland & Michaylov].
 *
 * Assumes that the function is evaluated by a predicate of the form
 *
 *	eval(T, Y, Yp), where T is the independent variable,
 *			Y the current vector of solutions,
 *			Yp the new vector derivative 
 *
 * e.g. eval(T, [Y1,Y2], [T*T + Y1*Y1, T*T - Y2*Y2]).
 *
 * Note: machine epsilon is available as #zero (pre-defined constant).
 *
 * Separate "mode" from the rest - may help efficiency
 */

/*
 * solve(Y, T, Tout, Incr, Relerr, Abserr)
 * Predicate called by the user - if he wants to specify a mode, he calls
 * solve/7 direct (only needed for single_step mode).
 *
 * Y - 		initial position vector
 * T -		starting point for independent variable
 * Tout - 	last point at which output is desired
 * Incr - 	default output increment 
 * Relerr - 	relative error tolerance
 * Abserr - 	absolute error tolerance
 */

solve(Y, Tstart, Tend, Incr, Relerr, Abserr) :-
	solve(Y, Tstart, Tend, Incr, Relerr, Abserr, normal).

/*
 *  Error checks
 */
solve(_, _, _, 0, _, _, normal) :-
	printf("Increment must be non-zero for normal mode\n",[]).
solve(_, _, _, _, Relerr, _, _) :-
	Relerr < 0,
	printf("Relative error must be non-negative\n",[]).
solve(_, _, _, _, _, Abserr, _) :-
	Abserr < 0,
	printf("Relative error must be non-negative\n",[]).
solve(_, T, T, _, _, _, _) :-
	printf("No interval\n",[]).

/*
 * Everything seems ok, so proceed
 */
solve(Y, T, Tend, Incr, Relerr, Abserr, Mode) :-
	set_yp(T, Y, Yp),		% Nfe = Nfe + 1
	output(T, Y),
	solve1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
	       work(Yp, 1, 0)).		% Work out the rest
	       %	Nfe = 1
	       %	Kop = 0

/*
 * Commence main iteration
 */
solve1(user(_, T, Tend, Incr, _, _, _), work(_, _, Nfe, _, _)) :-
	abs(T - Tend) < Incr / 4,
	printf("\nIteration finished\n------------------\n",[]),
	printf(" %d  derivative evaluations\n",[Nfe]).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, reset(Mode, relerr)), _) :-
	remin(Remin),
	Relerr < 2* #zero + Remin,
	printf("Warned you to reset it larger than %f\n", [Relerr]).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, reset(Mode, relerr)),
	work(Yp, H, Nfe, Kop, Tout)) :-
	remin(Remin),
	Relerr >= 2* #zero + Remin,
	solve1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
		work(Yp, H, Nfe, Kop, Tout)).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, reset(Mode, abserr, Oerr)), _) :-
	Abserr <= Oerr,
	printf("Abserr must be increased - not set to %f\n", [Abserr]).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, reset(Mode, abserr, Oerr)),
	work(Yp, H, Nfe, Kop, Tout)) :-
	Abserr > Oerr,
	solve1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
		work(Yp, H, Nfe, Kop, Tout)).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), work(Yp, Nfe, Kop)) :-
	% First time - need to set H
	set_tout(T, Tend, Incr, Tout),	% Next output point
	set_hinit(T, Tout-T, Y, Yp, Relerr, Abserr, H),
	set_kop(Kop, H, Tout-T, Newkop),
	stop_or_iter(user(Y, T, Tend, Incr, Relerr, Abserr, normal),
		     work(Yp, H, Nfe, Newkop, Tout)).
solve1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), 
	% H already set - continue on
	work(Yp, H, Nfe, Kop, _)) :-
	set_tout(T, Tend, Incr, Tout),
	sign(Tout-T, S),
	set_kop(Kop, S*H, Tout-T, Newkop),
	stop_or_iter(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
		     work(Yp, S*H, Nfe, Newkop, Tout)).

/*
 * Perform certain error checks
 */
stop_or_iter(User, work(Yp, H, Nfe, Kop, Tout)) :-
	Kop >= 100,
	check_out(User, work(Yp, H, Nfe, 0, Tout), output_excess).
stop_or_iter(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), Work) :-
	remin(Remin),
	Relerr < 2* #zero + Remin,
	check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), Work, relerr).
stop_or_iter(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
	     work(Yp, H, Nfe, Kop, Tout)) :-
	Abs = abs(Tout-T), 
	Tabs = abs(T),
	% Use Euler approximation
	Abs <= 26* #zero*Tabs,
	set_euler(Tout-T, Y, Yp, NewY),
	% Nfe = Nfe + 1
	set_yp(Tout, NewY, NewYp),		
	print_or_it(user(NewY, Tout, Tend, Incr, Relerr, Abserr, 
		mode(Mode, output)), work(NewYp, H, Nfe+1, Kop, Tout)).
stop_or_iter(User, Work) :-
	% Go on with main loop .......
	iter(User, Work).

iter(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
     work(Yp, H, Nfe, Kop, Tout)) :-
	set_modeh(Tout-T, H, NewH, Out),
	iter1(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, Out, nofail)),
	     work(Yp, NewH, Nfe, Kop, Tout)).

iter1(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, Out, Hmode)),
      work(Yp, H, Nfe, Kop, Tout)) :-
	max_func_eval(Maxfunc),
	Nfe >= Maxfunc,
	check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
	     work(Yp, H, Nfe, Kop, Tout), max_func).
iter1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
      work(Yp, H, Nfe, Kop, Tout)) :-
      fehl(Y, T, Yp, H, Sum, Abserr, Relerr, Errest, Ind),
      stop_or_iter1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
      		    work(Yp, H, Nfe+5, Kop, Tout), Sum, Errest, Ind).

stop_or_iter1(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, Out, Hmode)),
      	      Work, _, _, vanished) :-
	check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), Work, vanished).
stop_or_iter1(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
      	      work(Yp, H, Nfe, Kop, Tout), Sum, Esttol, okay) :-
	Esttol > 1,
	set_smin(Esttol, S),
	check_h(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
		work(Yp, S*H, Nfe, Kop, Tout)).
stop_or_iter1(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, Out, Hmode)),
      	      work(Yp, H, Nfe, Kop, Tout), Sum, Esttol, okay) :-
	Esttol <= 1,		
	% Nfe = Nfe+1
	set_yp(T+H, Sum, NewYp),
	set_hmax(Esttol, Hmode, T, H, NewH),
	print_or_it(user(Sum, T+H, Tend, Incr, Relerr, Abserr, 
		mode(Mode, Out)), work(NewYp, NewH, Nfe+1, Kop, Tout)).

print_or_it(user(Y, T, Tend, Incr, Relerr, Abserr, 
		mode(normal,proceed)), Work) :-
	iter(user(Y, T, Tend, Incr, Relerr, Abserr, normal), Work).
print_or_it(user(Y, T, Tend, Incr, Relerr, Abserr, 
		mode(normal, output)), Work) :-
	output(T, Y),
	solve1(user(Y, T, Tend, Incr, Relerr, Abserr, normal), Work).
print_or_it(user(Y, T, Tend, Incr, Relerr, Abserr, 
		mode(single_step, _)), Work) :-
	check_out(user(Y, T, Tend, Incr, Relerr, Abserr, single_step), Work).

check_h(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, _, _)),
	work(Yp, H, Nfe, Kop, Tout)) :-
	Habs = abs(H), 
	Tabs = abs(T),
	Habs <= 26* #zero*Tabs,
	printf("Tabs = %e , Habs = %e , #zero  = %e \n",[Tabs,Habs,#zero]),
	check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode),
		  work(Yp, H, Nfe, Kop, Tout), step_small).
check_h(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode, Out, _)), Work) :-
	iter1(user(Y, T, Tend, Incr, Relerr, Abserr, mode(Mode,proceed,hfail)),
	      Work).

/*
 * check_out(user(Y, Tstart, Tend, Incr, Relerr, Abserr, Mode),
 *           work(Yp, H, Nfe, Kop, Tout), Flag)
 * Checks the status variable Flag, and takes appropriate action.
 * See if there is a user-defined procedure
 */
check_out(User, Work, Flag) :- 
	user_error(User, Work, Flag), 
	!.
% No user procedure, so use the default ones.
check_out(user(Y, T, Tend, Incr, Relerr, Abserr, single_step),
       Work, single_step) :-
	printf("(single_step) T: %10.6f      Y: %10.6f\n",[T,Y]),
	solve(Y, T, Tend, Incr, Relerr, Abserr, single_step),!.
check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), Work, relerr) :-
	% Iflag = 3
	remin(Remin),
	Rer = 3* #zero + Remin,
	printf("Relative error too small - reset to %f\n", [Rer]),
	solve1(user(Y, T, Tend, Incr, Rer, Abserr, reset(Mode, relerr)), Work).
check_out(User, work(Yp, H, Nfe, Kop, Tout), max_func) :-
	% Iflag = 4
	max_func_eval(Max),
	Nfe >= Max,
	printf("Lots of function evaluations - %f needed so far\n",[Max]),
	solve1(User, work(Yp, H, 0, Kop, Tout)).
check_out(user(Y, T, Tend, Incr, Relerr, Abserr, Mode), Work, vanished) :-
	% Iflag = 5
	abmin(Absmin),
	printf("Abserr reset to %f\n", [Absmin]),
	solve1(user(Y, T, Tend, Incr, Relerr, 10*Abserr + Absmin, 
		reset(Mode, abserr, Abserr)), Work).
check_out(user(Y, T, Tend, Incr, Relerr, Abserr, normal), Work, step_small) :-
	% Iflag = 6
	printf("Relerr too small; reset to %f\n", [10*Relerr]),
	solve1(user(Y, T, Tend, Incr, 10*Relerr, Abserr, normal), Work).
check_out(User, work(Yp, H, Nfe, Kop, Tout), output_excess) :-
	% Iflag = 7
	printf("Excessive output; this will be inefficient",[]),
	solve1(User, work(Yp, H, Nfe, 0, Tout)).
check_out(_, _, stuffed) :- 
	% Iflag = 8
	printf("Improper call - wrong number of arguments?\n").
check_out(_, _, missing) :-
	printf("User-supplied predicate eval is missing or wrong\n"),
	printf("Format is eval(T, Y, Yp),\n",[]),
	printf("where T is the independent variable, ",[]),
	printf("Y the vector, and Yp its derivative.",[]),
	printf("Y and Yp are CLP lists - [Y1, Y2, Y3], [Yp1, Yp2, Yp3] etc.",[]).
check_out(_, _, Other) :-
	printf("Damned if I know what''s wrong!\n",[]).

/*
 * output(T, Y).
 * output(user(Y, T, Tend, Incr, Relerr, Abserr), work(Yp, H, Nfe, Kop, Tout))
 * Output predicate. If there is a user_output predicate (matching the above
 * format), then it is used. Otherwise, T and Y are written to standard output.
 */
output(T, Y) :-
	user_output(T, Y), 
	!.
output(T, Y) :-
	printf("T: %10.6f      Y: ",[T]),
	write_list(Y).

write_list([L|List]) :-
	printf(" %10.6f ",[L]), 
	write_list(List).
write_list([]) :-
	nl.

/*
 * Set routines - setting values of variables given certain conditions etc.
 */

set_yp(T, Y, Yp) :-	
	eval(T, Y, Yp).
set_yp(T, [Y], [Yp]) :-	
	% For one var case with no list functor
	eval(T, Y, Yp).	

set_hinit(T, Dt, Y, Yp, Relerr, Abserr, H) :-
	H1 = abs(Dt),
	set_maxtol(H1, Y, Yp, Relerr, Abserr, Maxtol),
	set_h1(H1, T, Dt, Maxtol, H).

set_maxtol(H, [Y|Yrest], [Yp|Yprest], Relerr, Abserr, Maxtol) :-
	Yabs = abs(Y),
	Tol1 = Relerr*Yabs + Abserr,
	Tol1 > 0,
	Ypabs = abs(Yp),
	set_maxh(Ypabs, H, Tol1, H1),
	set_maxtol(H1, Yrest, Yprest, Relerr, Abserr, Tol2),
	Maxtol = max(Tol1, Tol2).
set_maxtol(H, [Y|Yrest], [Yp|Yprest], Relerr, Abserr, Maxtol) :-
	Yabs = abs(Y),
	Relerr*Yabs + Abserr <= 0,
	set_maxtol(H, Yrest, Yprest, Relerr, Abserr, Maxtol).
set_maxtol(_, [], [], _, _, 0).

set_maxh(0, H, _, H).
set_maxh(Ypabs, H, Tol, NewH) :-
	Ypabs > 0,
	NewH = min(H, pow(Tol/Ypabs, 1/5)).	% From abs(Yp)*h**5 <= Tol

set_h1(H, T, Dt, Tol, 26* #zero*Max) :-
	Tol <= 0,
	Tabs = abs(T), 
	Dtabs = abs(Dt),
	Max = max(Tabs, Dtabs).
set_h1(H, T, Dt, Tol, NewH) :-
	Tol > 0,
	Tabs = abs(T), 
	Dtabs = abs(Dt),
	Max = max(Tabs, Dtabs),
	NewH = max(H, 26* #zero*Max).

set_modeh(Dt, H, H, proceed) :-
	Dtabs = abs(Dt), 
	Habs = abs(H),
	Dtabs >= 2*Habs.
set_modeh(Dt, H, Dt/2, proceed) :-
	Habs = abs(H), 
	Dtabs = abs(Dt),
	Dtabs < 2*Habs,
	Dtabs > Habs.
set_modeh(Dt, H, Dt, output) :-
	Habs = abs(H), 
	Dtabs = abs(Dt),
	Dtabs <= Habs.

set_kop(Kop, H, Dt, Kop+1) :-
	Habs = abs(H), 
	Dabs = abs(Dt),
	Habs >= 2*Dabs.
set_kop(Kop, H, Dt, Kop) :-
	Habs = abs(H), 
	Dabs = abs(Dt),
	Habs < 2*Dabs.

set_smin(Esttol, 0.9/(pow(Esttol, 1/5)) ) :- 
	Esttol < 59049.
set_smin(Esttol, 0.1) :- 
	Esttol >= 59049.

set_hmax(Esttol, Hmode, T, H, Sign*H1) :-
	set_smax(Esttol, Hmode, S),
	Habs = abs(H), 
	Tabs = abs(T),
	H1 = max(S*Habs, 26* #zero*Tabs),
	sign(H, Sign).

set_smax(Esttol, nofail, 5) :- 
	Esttol <= 0.0001889568.
set_smax(Esttol, nofail, 0.9/(pow(Esttol, 1/5)) ) :- 
	Esttol > 0.0001889568.
set_smax(Esttol, hfail, 1) :- 
	Esttol <= 0.0001889568.
set_smax(Esttol, hfail, S) :- 
	Esttol > 0.0001889568,
	S = min(0.9/(pow(Esttol, 1/5)), 1).

set_euler(Dt, [Y|Yrest], [Yp|Yprest], [NewY|NewYrest]) :-
	NewY = Y + Dt*Yp,
	set_euler(Dt, Yrest, Yprest, NewYrest).
set_euler(_, [], [], []).

/*
 * set_tout(T, Tend, Incr, Tout)
 * Tout is the next output point, T the current value, Tend the final point
 * and Incr the set increment.
 *
 * If there is a user-supplied predicate user_tout, then it is used.
 * Otherwise, we add Incr to T and proceed.
 * Note that Incr may be negative.
 */
set_tout(T, Tend, Incr, Tout) :-
	user_tout(T, Tend, Incr, Tout), 
	!.
set_tout(T, Tend, Incr, Tend) :-
	Abs = abs(Tend - (T + Incr)),
	Abs < Incr.
set_tout(T, Tend, Incr, T + Incr) :-
	Abs = abs(Tend - (T + Incr)),
	Abs >= Incr.

/*
 * fehl(Y, T, Yp, H, Sum, Abserr, Relerr, Errest, Errind)
 * Predicate to perform the evaluation of Fehlberg formulae
 * 
 * Sum is the new estimate for Y
 * Abserr is the absolute error.
 * Relerr is the relative error.
 * Errest is the Fehlberg estimate of the error.
 * Errind indicates whether Abs is too small or not.
 * If so, Errind is set to "vanished". Otherwise, it is "okay".
 */
fehl(Y, T, Yp, H, S, Abs, Rel, Err, Ind) :-
	Ch1 = H/4, 
	Ch2 = H*3/32, 
	Ch3 = H/2197,
	Ch4 = H/4104, 
	Ch5 = H/20520, 
	Ch6 = H/7618050,
	fehl1(Y, T, Yp, H, ch(Ch1, Ch2, Ch3, Ch4, Ch5, Ch6), 
		S, Abs, Rel, Err, Ind).

fehl1(Y, T, Yp, H, ch(Ch1, Ch2, Ch3, Ch4, Ch5, Ch6), S, Abs, Rel, Err, Ind) :-
	% calculate each of the Fi .....
	set_f1(Y, Yp, Ch1, P1),			set_yp(T+Ch1, P1, F1),
	set_f2(Y, Yp, Ch2, F1, P2), 		set_yp(T+ 3*H/8, P2, F2),
	set_f3(Y, Yp, Ch3, F1, F2, P3), 	set_yp(T + 12*H/13, P3, F3),
	set_f4(Y, Yp, Ch4, F1, F2, F3, P4),  	set_yp(T + H, P4, F4),
	set_f5(Y, Yp, Ch5, F1, F2, F3, F4, P5),	set_yp(T + H/2, P5, F5),
	% .. then the Y estimate ....
	set_sum(Y, Yp, Ch6, F2, F3, F4, F5, S),
	% .. and finally the error
	set_err(Y, Yp, Abs, Rel, S, f(F2, F3, F4, F5), H, 0, Err, Ind).

set_f1([Y|Yr], [Yp|Ypr], Ch, [P|Pr]) :-
	P = Y + Ch*Yp,
	set_f1(Yr, Ypr, Ch, Pr).
set_f1([], [], _, []).

set_f2([Y|Yr], [Yp|Ypr], Ch, [F1|F1r], [P|Pr]) :-
	P = Y + Ch*(Yp + 3*F1),
	set_f2(Yr, Ypr, Ch, F1r, Pr).
set_f2([], [], _, [], []).

set_f3([Y|Yr], [Yp|Ypr], Ch, [F1|F1r], [F2|F2r], [P|Pr]) :-
	P = Y + Ch*(1932*Yp + (7296*F2 - 7200*F1)),
	set_f3(Yr, Ypr, Ch, F1r, F2r, Pr).
set_f3([], [], _, [], [], []).

set_f4([Y|Yr], [Yp|Ypr], Ch, [F1|F1r], [F2|F2r], [F3|F3r], [P|Pr]) :-
	P = Y + Ch*((8341*Yp - 845*F3) + (29440*F2 - 32832*F1)),
	set_f4(Yr, Ypr, Ch, F1r, F2r, F3r, Pr).
set_f4([], [], _, [], [], [], []).

set_f5([Y|Yr], [Yp|Ypr], Ch, [F1|F1r], [F2|F2r], [F3|F3r], [F4|F4r], [P|Pr]) :-
	P = Y + Ch*((0-6080*Yp + (9295*F3 - 5643*F4)) + (41040*F1-28352*F2)),
	set_f5(Yr, Ypr, Ch, F1r, F2r, F3r, F4r, Pr).
set_f5([], [], _, [], [], [], [], []).

set_sum([Y|Yr], [Yp|Ypr], Ch, [F2|F2r], [F3|F3r], [F4|F4r], [F5|F5r], [S|Sr]) :-
	S = Y + Ch*((902880*Yp + (3855735*F3 - 1371249*F4)) + (3953664*F2 + 
		    277020*F5)),
	set_sum(Yr, Ypr, Ch, F2r, F3r, F4r, F5r, Sr).
set_sum([], [], _, [], [], [], [], []).

set_err([Y|Yr], Yp, Abs, Rel, [S|Sr], f(F2, F3, F4, F5), H, Oerr, Err, Ind) :-
	Yabs = abs(Y), 
	Sabs = abs(S),
	Et = Yabs + Sabs + Abs,
	check_err(Et, [Y|Yr], Yp, [S|Sr], f(F2, F3, F4, F5),
		Abs, Rel, H, Oerr, Err, Ind).
set_err([], [], _, Rel, [], _, H, Err, Habs*Err*(2/Rel)/752400, okay) :-
	Habs = abs(H).

check_err(Et, Y, Yp, S, _, _, _, _, _, _, vanished) :-
	Et <= 0.
check_err(Et, [Y|Yr], [Yp|Ypr], [S|Sr], f([F2|R2], [F3|R3], [F4|R4], [F5|R5]),
		Abs, Rel, H, Oerr, Err, Ind) :-
	Et > 0,
	Ee = abs( (0-2090*Yp + (21970*F3 - 15048*F4)) + (22528*F2 - 27360*F5)),
	Nerr = max(Oerr, Ee/Et),
	set_err(Yr, Ypr, Abs, Rel, Sr, f(R2, R3, R4, R5), H, Nerr, Err, Ind).

% Low level stuff - primitives etc.

sign(V, 1) :-	
	V >= 0.
sign(V, 0-1) :-	
	V < 0.

% Magic numbers
remin(0.000000000001).	% 1.0e-012 - from FM&M.
abmin(0.000000001).	% 1.0e-009 - from FM&M.
max_func_eval(3000).	% Also from FF&M.

eval(T, [Y1, Y2, Y3, Y4], [Y3, Y4, 0-Y1/R, 0-Y2/R]) :-
	R1 = Y1*Y1 + Y2*Y2,
	R = R1*pow(R1, 0.5)/((#p/4)*(#p/4)).

user_output(T, [Y1, Y2, Y3, Y4]) :-
	printf("Point %10.5f : %10.5f %10.5f\n",[T,Y1,Y2]).

user_error(user(Y, T, Tend, Incr, Rer, Abs), 
	work(Yp, H, Nfe, Kop, Tout), max_func) :-
	printf("Too hard for me - try something more accurate\n",[]).

user_error(user(Y, T, Tend, Incr, Rer, Abs), 
	work(Yp, H, Nfe, Kop, Tout), output_excess) :-
	printf("Outputs too often - try another method\n",[]).

user_error(user(Y, T, Tend, Incr, Rer, Abs), 
	work(Yp, H, Nfe, Kop, Tout), step_small) :-
	printf("Can''t achieve required accuracy\n",[]).

go:- Ecc = 0.25, solve([1-Ecc, 0, 0, (#p/4)*pow((1+Ecc)/(1-Ecc), 0.5)],
			0, 3, 0.5, 0.000000001, 0).

% Output:
%
%  Point    0.00000 :    0.75000    0.00000
%  Point    0.50000 :    0.61977    0.47779
%  Point    1.00000 :    0.29442    0.81218
%  Point    1.50000 :   -0.10518    0.95804
%  Point    2.00000 :   -0.49030    0.93987
%  Point    2.50000 :   -0.81394    0.79959
%  Point    3.00000 :   -1.05403    0.57571
%  
%  Iteration finished
%  ------------------
%   439  derivative evaluations

?- printf("\n>>> Sample goal: go/0\n", []).