% Source Coding -- Problem 3.1 Arithmetic coding

function Arithmetic

clear;
EPS = 0.0000000001;             % define EPS for rounding
p=[0.08 0.4 0.3 0.2 0.02];      % vector of probabilities of each mark
low(1)=0.0;                     % set the low limit of mark 1
up(1)=p(1);                     % set the up limit of mark 1

for i=2:length(p)               
    low(i)=up(i-1);             % set the low limits from mark 2 to 5
    up(i)=low(i)+p(i);          % set the up limits from mark 2 to 5
end

l=0.0;                          % initialize the tag interval
u=1.0;                         
low_in=0.0;
up_in=1.0;

code=fileread('exam.txt');      % read in the codes
out=zeros(1,500);               % vector to store the 500 marks
count=1;                        % counter for marks
k=1;                            % pointer for the codes

while count<=500
    if k<=length(code)
        if code(k)=='0'                     % the situation of '0'
            up_in=(low_in+up_in)/2;         % modify the up limit
            up_in=round(up_in*1000000+EPS)/1000000;
        else                                % the situation of '1'
            low_in=(low_in+up_in)/2;        % modify the low limit
            low_in=round(low_in*1000000+EPS)/1000000;
        end
        k=k+1;
    else

    
        up_in=(low_in+up_in)/2;            % assume that several zeros exist at the end
        up_in=round(up_in*1000000+EPS)/1000000;
    end
   
    % to find the interval
    [mark,l_temp,u_temp]=search_interval(low,up,low_in,up_in);

    if mark==0                              % do not find the interval, repeat the loop
        continue;
        
    else
        out(count)=mark;                    % read out the mark
        count=count+1;
        l=l_temp;
        u=u_temp;
    end
    
    while (l >=0.5 | u <0.5)                % the position that both the limits are in one half of 1
        temp=round(u);
        l=2*l-temp;
        l=round(l*1000000+EPS)/1000000;
        u=2*u-temp;
        u=round(u*1000000+EPS)/1000000;
        low_in=2*low_in-temp;
        low_in=round(low_in*1000000+EPS)/1000000;
        up_in=2*up_in-temp;
        up_in=round(up_in*1000000+EPS)/1000000;
        
    end

    low(1)=l;                               % to recalculate the limits of each mark
    up(1)=l+(u-l)*p(1);
    
        for i=2:length(p)               
            low(i)=up(i-1);             
            low(i)=round(low(i)*1000000+EPS)/1000000;
            up(i)=low(i)+(u-l)*p(i);
            up(i)=round(up(i)*1000000+EPS)/1000000;
        end
        
end

out(9)
out(150)
out(376)
out(500)

sum(out==1)
sum(out==2)
sum(out==3)
sum(out==4)
sum(out==5)