%       program hd1dP1i.m
%       Current version: June 27, 2017
   %       Past version: November 30, 2015 (hd1dP1.m)
%       Computes Hausdorff dimension of the Cantor set 
%       E_{a_1, a_2,...,a_R}
%       consisting of all points x in the interval (0,1) with
%       the continued fraction expansion of
%       x = [a_1, a_2, a_3, ...] =                  1
%                                        ----------------------
%                                          a_1 +     1
%                                                --------
%                                                a_2 +    1
%                                                     ---------
%                                                     a_3 + ...
%
%      where each a_n is one of the integers a_1, ... a_R.
%
%      Uses the routines mat1d.m  and thetabx.m
%
%      We approximate by continuous piecewise linear functions,
%
%
%  xlen = length of eigenvector
%  N+1 = number of mesh points;  Choose N to be divisable by 2
%  R = maximum size of b used in the computation
%  i.e., the number of a_i in the set E.
%  al is the lower matrix, aa is the approx matrix, au is the upper matrix
%  mual is the largest eigenvalue of al,
%  muaa is the largest eigenvalue of aa
%  muau is the largest eigenvalue of au
%  
%    If we find values s_1 for which mual >= 1 and s_2 for which muau <= 1,
%    then the true largest eigenvalue (giving the Hausdorff dimension of K)
%    satisfies  s_1 <= s <= s_2.
%
% Input data:
       format long
%      Note  N = int16(32766) still had a memory problem.
      prompt = 'Input N, the number of subintervals ';  
      N = input(prompt)
%
%
%%%%%   The value of R must be input for each run
%%%%%   To do E[1,3], input R = 2.
   prompt = 'Input R, the number of coefficients in the continued fraction ';  
      R = input(prompt)
      blen = R

%%%%% Input the values of the vector b in increasing order 
%%%%% one at a time at the prompt. To do E[1,3], input 1 at the first
%%%%% prompt and 3 at the second prompt.
        b = zeros(blen,1);
        for j =1: blen
           prompt = 'Input the values of the vector b in increasing order one at a time after each prompt.'; 
           b(blen+1-j) = input(prompt);
        end
        display(b)
       c = 0.0;
       d = 1.0/b(blen);
       display(c)
       display(d)

% Calculate the mesh size
      h = (d-c)/double(N)
      xlen = N+1
%
% Index the mesh points in the domain by an index j
% of length xlen starting at 1.
%
%%        x = zeros(xlen,1);
%%       for j = 1: N +1
%%	    x(j) = c + double(j-1)*h;
% %       end
       x = linspace(c, d, N+1);
%
% For each of the mesh points x and values of m, compute the functions
% aphib = 1/|x + m|
%
        aphib = zeros(xlen,blen);
        for k = 1: xlen
            for l = 1:blen 
		aphib(k,l) = 1/abs(x(k) + b(l));
            end
        end
%
%  For each of the points thetab(x) = thetabx, 
%  find the index number ix of the left mesh point
%  of the subinterval to which it belongs.
%  Note that we must decide how to choose the subinterval 
%  when the points lie on the boundary of two subintervals
%
	ix = zeros(xlen,blen);    
        for k = 1:xlen
            for l = 1:blen
		  thetabxkl = thetabx(x(k),b(l));
	          ix(k,l) = floor((thetabxkl-c)/h) + 1;
            end
%
%  If the index is outside the allowable limits, print the results
         if (ix(k,l) < 1)
	       display(['(ix(k,l) < 1'])
               display([k, l, ix(k,l), thetabxkl])
         end
         if (ix(k,l) > xlen)
	       display(['ix(k,l) > xlen'])
               display([k, l, ix(k,l), thetabxkl])
         end
% If the index is equal to xlen, then reset to one less
         if (ix(k,l) == xlen)
	   ix(k,l) = ix(k,l) -1;
          end

      end
%      display(ix)
%   We next compute the Lagrange basis functions Lx corresponding 
%   to each of these degrees of freedom at the point thetabx
%
for k = 1:xlen
  for l = 1: blen
    ixkl = ix(k,l);
    thetabxkl = thetabx(x(k),b(l));     
    Lx(k,l,2) = (thetabxkl - x(ixkl))/h;
    Lx(k,l,1) = 1.0 - Lx(k,l,2);
  end
end
%
%
% Compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%  
% First do the computation for the approximate matrix
mtype = 0;
funa = @(s) log(eigs(mat1d(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx),1));
sa = fzero(funa, [0.01, 0.99]);
muaa = eigs(mat1d(mtype, sa, h, ix, x, b, aphib, xlen, blen, Lx),1);
muaa1 = muaa - 1;
disp ('Approximation of Hausdorff dimension s.')
sa
disp ('The eigenvalue of the  approximate matrix minus one.')
muaa1
%
% Next do the computation for the lower matrix
mtype = -1;
funl = @(s) log(eigs(mat1d(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx),1));
sl = fzero(funl, sa);
As = mat1d(mtype, sl, h, ix, x, b, aphib, xlen, blen, Lx);
[val, mual] = eigs(As,1);
mual1 = mual - 1;
% Add loop to ensure mual > 0  and min(As*val - val) >=0
% First make sure that val is a strictly positive eigenvector
if (max(val) <= 0) val = -val;
end
minval = min(As*val - val);
while  (mual1 <=0 || minval < 0)  
   sl = sl - eps;
   As = mat1d(mtype, sl, h, ix, x, b, aphib, xlen, blen, Lx);
   [val, mual] = eigs(As,1);
   mual1 = mual - 1;
   if (max(val) <= 0) val = -val;
   end
minval =   min(As*val - val);
end
disp ('Lower bound on the Hausdorff dimension s.')
sl
disp ('The eigenvalue of the  lower matrix minus one.')
% disp ('The value should be >=0.')
mual1
% The following value should be >=0
minval
% Next do the computation for the upper matrix
mtype = 1;
funu = @(s) log(eigs(mat1d(mtype, s, h, ix, x, b, aphib, xlen, blen,Lx),1));
su = fzero(funu, sa);
Bs = mat1d(mtype, su, h, ix, x, b, aphib, xlen, blen,Lx);
[vau, muau] = eigs(Bs,1);
muau1 = muau - 1;
% Add loop to ensure muau1 <0 and min(vau - Bs*vau) >=0
% First make sure that vau is a strictly positive eigenvector
if (max(vau) <= 0) vau = -vau;
end
minvau = min(vau - Bs*vau);
while (muau1 >= 0 || minvau < 0)  
   su = su +  eps;   % eps is machine epsilon
   Bs = mat1d(mtype, su, h, ix, x, b, aphib, xlen, blen,Lx);
   [vau, muau] = eigs(Bs,1);
   muau1 = muau -1;
   if (max(vau) <= 0) vau = -vau;
   end
   minvau =  min(vau - Bs*vau);
end
disp ('Upper bound on the Hausdorff dimension s.')
su
disp ('The eigenvalue of the  upper matrix minus one.')
% disp ('The value should be <=0.')
muau1
% The following value should be >=0
minvau
%       end program hd1dP1i.m


