%% this file demonstrates the implementation of integer arithmetic
%Exa = 1;

switch Exa
case 1
x = linspace(0,2,1001);
f = @(x)(1-0.8*x.^2);

y = f(x);
figure(1)
plot(x,y);
title('original function');
xlabel('x'); ylabel('y = 1-0.8*x*x');
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

case 2
%% use the above graph to verify that -3<y<1

%% we aim for an 8-bit implementation
%% thus the 8-bit integer result y8 should satisfy -127<y8<127
%% as consequence we compute y8 = 2^5*y = 32*y

%% the values of x vary from 0 to 2, thus we pass x8 = 64*x=2^6*x
%% as parameters to the result
x8 = int8(64*x);

%% compute r1=x*x 
r1 = int16(int16(x8).*int16(x8));
%% verify that the limit -2^15 and 2^15 are respected
lim1 = [min(r1),max(r1),2^15]
deviation1 = max(abs(double(r1)-x.^2*2^28))
figure(11)
plot(x,x.^2,x,double(r1)/2^12)
xlabel('x'); ylabel('x^2');
title('test1')
legend('x^2','approx')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

case 3
%% use only the high byte
%% then multiply by 0.8, since 0.8*256=204.8, multiply by 205
r2   = int16(205*bitshift(r1,-8));
lim2 = [min(r2),max(r2),2^15]
deviation1 = max(abs(double(r2)-x.^2*2^12))
figure(12)
plot(x,0.8*x.^2,x,double(r2)/2^12)
xlabel('x'); ylabel('x^2');
title('test4')
legend('0.8*x^2','approx')

%% rescale to obtain an 8bit result
r3   = int8(bitshift(r2,-7));
lim3 = [min(r3),max(r3),2^7]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

case 4
%% subtract the results r3 from 1, resp from 1*32
r4   = int8(32-r3);
lim4 = [min(r4),max(r4),2^7]
deviation4 = max(abs(double(r2)-(1-0.8*x.^2)*2^12))
figure(14)
plot(x,1-0.8*x.^2,x,double(r4)/2^5)
xlabel('x'); ylabel('x^2');
title('test4')
legend('1-0.8*x^2','approx')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

case 5
%% a plot of the error
r5 = single(r4)/2^5;
figure(2);
plot(x,r5-y)
title('arithmetic error')
xlabel('x'); ylabel('error'); grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

case 6
%% compute the relative error
relErr       = max(abs(r5-y))/max(abs(y))
bitError     = log2(relErr)
bitErrorMean = log2(mean(abs(r5-y))/max(abs(y)))
%% the result shows a worst case accuracy of 5 bit 
%% and a mean accuracy of 6 bit
end%switch
