% started with a review of fl.#s
% along the line of the end of last day's diary (sep14):

>> format hex
>> 1
ans = 3ff0000000000000

>> (3*16+15)*16+15
ans = 408ff80000000000

>> format
>> ans
ans = 1023

>> format hex
>> 2
ans = 4000000000000000

>> inf
ans = 7ff0000000000000

>> -ans
ans = fff0000000000000

>> 2^1024
ans =  7fe0000000000000

>> ans*2
ans = 7ff0000000000000

% conclusion: the largest exponent available in matlab's fl.#s is 1024
% the next larger exponent, 1025, is reserved for inf = infinity
% Check it, by looking at (.1)_2 2^1024 = 2^1023, then multiplying that by 2:
>> format
>> 2^1023
ans = 8.9885e+307
>>  ans*2
ans = Inf

% catastrophic cancellation illustrated various ways in book and notes.
% Remedy: replace numerical subtraction of near-equal quantities by symbolic
% subtraction. In particular, replace numerical differentiation by symbolic
% (exact) differentiation. If the function is only known by way of a program,
% differentiate the program.

% Here is an example, evaluation of a polynomial in power form from its
% coefficients, using nested multiplication (aka Horner's method):

% Here is the function m-file for that calculation I started with. 
% Input is the coefficient vector, and the x at which the value is wanted:

>> type p

function v = p(a,x)

v = a(1);
for k=2:length(a)
    v = v*x + a(k);
end

% Now I'll differentiate this program with respect to x, using the fact that
% the derivative of  x  wrto  x   is  1, and making sure that the
% differentiated statement precedes the statement differentiated:

>> type p

function [v,dv] = p(a,x)

dv = 0;
v = a(1);
for k=2:length(a)
    dv = dv*x + v;
    v = v*x + a(k);
end

% Finally, I would like this function to work correctly even for a vector or
% matrix  x  , and this I can do by making certain that any multiplication
% (division, exponentiation; except there are none here) is done pointwise:


>> type p

function [v,dv] = p(a,x)

dv = 0;
v = repmat(a(1),size(x));
for k=2:length(a)
    dv = dv.*x + v;
    v = v.*x + repmat(a(k),size(x));
end


% We are ready to try it out.
% I choose the coefficient vector
>> a = [1 0 0];

% Hence, p(x) = 1*x^2 + 0*x + 0 = x^2, therefore p'(x) = 2x. We try it out,
% plotting the results:

>> x = -1:.01:1;

>> [v,dv] = p(a,x);
>> plot(x,v,x,dv)

% the plot shows that it works.

diary off