echo on
%%%%% ch01_pro.m %%%%%
% Run script and study the effects of each line.

format compact
%************************************************
%**** Some Computer Problems from Chapter 01 ****
%************************************************

%**** COMPUTER PROBLEM 1.8 *****
%Calculate the sum for a series: 1+1/2+1/3+...+1/i, i=1,2,....

echo off
suma=0;
fac=1;
i=1;
while  (fac >= eps*suma & i<10000)
   fac=1/i;
   suma = suma + fac;
 %  disp(sprintf('i=%d    1/i=%1.15e    sum=%1.15e',i,1/i,suma));
   i=i+1;
end;
disp(sprintf('i=%d    1/i=%1.15e    sum=%1.15e',i,1/i,suma));
echo on

%The series sum(1/n) is divergent, but in flp arithmetic it has the finite 
%sum because the contribution of 1/i becomes negligible relative to the 
%partial sum. The stopping criterion should be 1/i < eps*(partial sum).
%We would need about 10^6 iterations in sigle precision and much more 10^(13)
%in double precision!


%**** COMPUTER PROBLEM 1.9 *****
%Compute e^x with the series: 1+x+x^2/(1*2)+x^3/(1*2*3)+....

x=10;
exp(x)		%Expected value
echo off
suma=1;
cont=1;
fact=1;
i=1;
while  (cont >= eps*suma)
   fact=i*fact;
   cont=x^i/fact;
   suma = suma + cont;
   disp(sprintf('i=%d    cont=%1.15e    sum=%1.15e',i,cont,suma));
   i=i+1;
end;
echo on

%The stopping criterion is when the actual factor < eps*(partial sum).
%For x<0 results are not correct.
%We can use the formula e^-x = 1/e^x to avoid negative exponents, or
%Rearrange the summation (take first term and add last term, proceed with 
%second term and term before last term and so on.


%**** COMPUTER PROBLEM 1.10 *****
%Solve quadratic equation ax^2+bx+c=0 with the standard and alternate 
%formula for different values of coefficients.

a=1;
b=-10^10;
c=1;
echo off
% re-scale to avoid unnecessary underflow or overflow
if (abs(a)>abs(b) & abs(a)>abs(c)) b=b/a;c=c/a;a=1; 	%re-scale by a
elseif (abs(b)>abs(a) & abs(b)>abs(c)) a=a/b;c=c/b;b=1;	%re-scale by b
else a=a/c;b=b/c;c=1; %re-scale by c
end	

%complex square root ?
sq=b^2-(4*a*c);
if sq<0 disp(sprintf('Complex roots')); sq=abs(sq); % complex roots, no cancellation possible
	sqs=sq^0.5;
   xc= -b/(2*a); yc= sqs/(2*a);
   disp(sprintf('Complex roots: %g +- i*%g', xc, yc));
% real squre roots   
else if sq<eps disp(sprintf('Possible cancellation under square root: sq=%g', sq)); end  
	sqs=sq^0.5;
	x1=-b+sqs;	%first root
	if (abs(x1)<eps | a==0) disp(sprintf('Trying to find 1st root by the alternative formula (x1=%g).',x1));
   	x1=-b-sqs;
   	disp(sprintf('First root (alternative formula):  x1=%g', (2*c)/x1));
   else disp(sprintf('First root (standard formula):  x1=%g', x1/(2*a))); end
   if (a == 0) disp(sprintf('There is only a single root.'));
      else
      	x2=-b-sqs;	% second root
			if (abs(x2)<eps)disp(sprintf('Trying 2nd root by the alternative formula (x2=%g).',x2));
   			x2=-b+sqs;
      		disp(sprintf('Second root (alternative formula):  x2=%g', (2*c)/x2));
      	else disp(sprintf('Second root (standard formula):  x2=%g', x2/(2*a)));end
   end
end

echo on

%We should re-scale first the coefficients to avoid unnecessary underflow or overflow.
%Then we test if the square roots are complex (b^2-4ac<0).
%If not, the program test eventual unavoidable cancellation under square root
%and try to find both roots with the standard and alternative formulas.
%The decision criterion for the formula selection is the magnitude of fraction < eps.


%**** COMPUTER PROBLEM 1.12 *****
%calculate and print MatLab's std1 and my_std1 (two-pass formula)
% and my_std2 (one-pass formula). 

nn=100;		% Number of data elements
echo off
% MathLab's std
va=10^6*ones(nn,1);	% generate a vector of slightly different elements
vb=rand(nn,1);
for i=1:nn 
   if vb(i)>0.5 vb(i)=vb(i)-1; end; 
end;
vb=vb*10^-3;
vector=va+vb;			

std1= std(vector);	% MatLab's standard deviation 

% two sumation-pass formula
mean1=mean(vector);	% calculate mean (1.pass)
my_std1=0;
for i=1:nn 				% second pass
   my_std1=my_std1+(vector(i)-mean1)^2;
end
my_std1=sqrt((1/(nn-1))*my_std1);


% one summation-pass formula (cancellation by large similar numbers) 
my_std2=0;
mean2=0;  
for i=1:nn 
   mean2=mean2+vector(i);
   my_std2=my_std2+vector(i)^2;
end		% take abs for positive squared value
mean2=mean2^2/nn;
my_std2=sqrt(abs((1/(nn-1))*(my_std2-mean2))); 

disp('                           std');
disp(sprintf('Mat_lab''s:            %-15g %10d', std1));
disp(sprintf('Two-pass formula:      %-15g %10d', my_std1));
disp(sprintf('One-pass formula:      %-15g %10d', my_std2));
echo on

%If the elements of the input vector are app 10^6 and if they differ only for
%10^-3, one-pass formula gives incorrect results.