%       program hd1dPni.m
%       Current version: July 3, 2017
%       Previous version: December 4, 2015
%       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 eval1dPn.m and thetabx.m
%
%    We approximate by continuous piecewise polynomials of degree <= deg, 
%     deg > 1, so the current theory does not apply.
%
%  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.
%  aa is the approximate matrix
%  muaa is the largest eigenvalue of aa
%  
%    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
      c = 0.0;
      d = 1.0;
      display(c)
      display(d)
      prompt = 'Input N ';  
      N = input(prompt)
%
% deg = the degree of the piecewise polynomial
      prompt = 'Input degree of piecewise polynomial ';  
      deg = input(prompt)
%
% Calculate the mesh size
% On each subinterval of length h, we will have meshpoints at the
% endpoints and at deg -1 equally spaced points in each subinterval.  Then the
% total number of mesh points xlen will be  deg*N + 1

      h = (d-c)/double(N)
      xlen = deg*N+1

%%%%%   The value of R must be input for each run
%%%%%   To do E[1,3], input R = 2.
      prompt = 'Input R ';  
      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 values of vector b in increasing order one at a time after each prompt.'; 
            b(blen+1-j) = input(prompt);
        end
        display(b)

% Index the mesh points in the domain by an index j
% of length xlen starting at 1.
%
        x = zeros(xlen,1);
        for j = 1: deg*N +1
	    x(j) = c + double(j-1)*h/deg;
        end
%
% 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.
%  There are N subintervals with each one having deg-1 interior mesh points.
%  Note that we must decide how to choose the subinterval 
%  when the points lie on the boundary of two subintervals
%
ix = zeros(xlen,blen);
Lx = ones(xlen,blen,deg+1);
for k = 1:xlen
   for l = 1:blen
      thetabxkl = thetabx(x(k),b(l));
      ix(k,l) = deg*floor((thetabxkl-c)/h) + 1;
%  If the index is outside the allowable limits, print the results
      if (ix(k,l) < 1)
        display(['(ix(k,l) < 1', k, l, ix(k,l), thetabxkl])
      end
      if (ix(k,l) > xlen)
        display(['ix(k,l) > xlen', k, l, ix(k,l), thetabxkl])
      end
% If the index is equal to xlen, then reset to deg less
      if (ix(k,l) == xlen)
         ix(k,l) = ix(k,l) - deg;
      end
%   We next compute the Lagrange basis functions Lx corresponding 
%   to each of these degrees of freedom at the point thetabx
%
     ixkl = ix(k,l);
     for i = 1: deg+1
       for j = [1:i-1 i+1:deg+1]
	     Lx(k,l,i) = (thetabxkl - x(ixkl + j-1))/((i-j)*h/deg)*Lx(k,l,i);
       end
     end
  end
end
%

%
% We compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%  

fun = @(s) log(eval1dPni(s, ix, x, b, aphib, xlen, blen, Lx,deg));
disp ('Value of Hausdorff dimension s.')
s = fzero(fun, [0.4,0.7])
mua = eval1dPni(s, ix, x, b, aphib, xlen, blen, Lx, deg);
% disp ('The eigenvalue of the matrix minus one.')
% disp ('The value should be approximately 0.')
mua1 = mua - 1;
%       end program hd1dPni.m
