%-----------------------------------------------------------------------------
%               ------< Problem 2.2 (Text 5.2.7) >------
%-----------------------------------------------------------------------------
% script hw2_2
%             driver for function MatPower
% asks for a matrix dimension, m, and a power, n. Defines a random mXm 
% matrix A, prints the binary expansion of n (invokes to script bin.m), 
% and estimates the number of flops that would be required to compute
% A^n. Then it calls the function MatPower, which computes the product
% efficiently.
% NOTE: the estimates given here were NOT REQUIRED! Just for fun...
m = input('enter matrix dimension:');
A = randn(m,m);
n = input('enter power:');
disp('--------------------------------')
disp(sprintf(' binary expansion for  %6.0f ',n))
y = bin(n)
disp('--------------------------------')
disp(sprintf('estimate of work required to raise %4.0f X %4.0f matrix to %5.0f-th power',m,m,n))
k = 2*m^3*n;
disp(sprintf('  (1)   By direct computation, flops ~ %10.0f ',k))
k=length(y);
s=0;
for i = 1:k
   s = s+y(i);
end
k = s+k;
k = 2*m^3*(k-1);
disp(sprintf('  (2)   By efficient computation, flops ~ %10.0f ', k ))

% Direct computation of A^n;
   B = eye(m,m);
   flops(0);
for i = 1:n
   B = B*A;
end
   f1 = flops;
% Computation via MatPower
   flops(0);
   C = MatPower(A,n);
   f2 = flops;
disp(sprintf(' Actual cost: direct %10.0f vs. MatPower %10.0f ',f1,f2))
E=C-B;
error = norm(E,inf)/norm(C,inf);
disp(sprintf('relative error in computation %22.15f ',error))
%------------------------------------------------------------------------------
% sample output of script:
%>> hw2_2
%enter matrix dimension:67
%enter power:55
%--------------------------------
% binary expansion for      55
%
%y =
%
%     1     1     0     1     1     1
%
%--------------------------------
%estimate of work required to raise   67 X   67 matrix to    55-th power
%  (1)   By direct computation, flops ~   33083930
%  (2)   By efficient computation, flops ~    6015260
% Actual cost: direct   33083930 vs. MatPower    6015307
%relative error in computation      0.000000000000003
%------------------------------------------------------------------------------
%   Functions called by script; included for completeness
%------------------------------------------------------------------------------
%              ------< function y = bin(n) >------
%------------------------------------------------------------------------------
% function y = bin(n)
% m=1;
% remainder=0;
% s = 0;
% while (remainder ~= n)
%    s = s+1;
%    x = rem((n-remainder)/m,2);
%    remainder = remainder + x*m;
%    m = m*2;
%    y(s) = x;
% end
% z=y;
% for i=1:s
%  y(s+1-i)=z(i);
% end
%------------------------------------------------------------------------------
%              ------< function C = MatPower(A,n) >------
%------------------------------------------------------------------------------
% function C = MatPower(A,n)
% p=size(A);
% l=p(1,1);
% m=1;
% remainder=0;
% B=A;
% C=diag(ones(1,l));
% while (remainder ~= n)
%    x = rem((n-remainder)/m,2);
%    remainder = remainder + x*m;
%    if (x ~= 0) C = B*C;end
%    if (remainder ~= n)
%       B = B*B;m=2*m;
%    end
% end
%------------------------------------------------------------------------------
