% use iterative taylor approximation to get three correct decimals of ln(2)

symbols;

%note: x=2
goal = 1 / (1 + 2 ^ 2)

n = 0;
x=sym("x");

for n=1:100
	n++;
	taylorexp = taylor(1 / (1 + x ^ 2), x, 0, n);
	result = vpa( subs(taylorexp, x, 2) );
	printf("%d	%d	%d\n", n, real(result), real(result) - goal);
end

%
% Produces a taylor approximation of the function to the given power
%
%	func = function to approximate
%	wrt = with repest to this variable
%	x_0 = base (a)
%	degree = desired degree of resulting polynomial
%
function Pt = taylor( func, wrt, x_0, degree )

	Pt = subs( func, wrt, x_0 );

	for k=1:degree
		dx = differentiate( func, wrt, k );

		# We create e here to force that GNC understands the exponent
		# as an integer and expands the expression
		e = vpa( sprintf( "%d", k ) );
		Pt += subs( dx, wrt, x_0 ) * expand( (wrt - x_0)^e ) / prod(1:k);
	end
	
endfunction
