function cantor  

% computes the function that matches numbers in [0,1] to the Cantor set
% via [x]_2  --> [y]_3  by changing  '1'   to '2' in [x]_2

x1=0;
x2=1;
M=10000;                % number of points in [0,1] to plot
del=(x2-x1)/M;
x=x1:del:x2;           % x(1), ..., x(M)  are the points in [0,1] to plot

b=2;    % base 2
d=3;    % base 3

for k=1:M                   %  [x(k)]_2 = .b1 b2 b3 .... b15
b1=floor(b*x(k));
b2=floor(b^2*x(k)-b*b1);
b3=floor(b^3*x(k)-b^2*b1-b*b2);
b4=floor(b^4*x(k)-b^3*b1-b^2*b2-b*b3);
b5=floor(b^5*x(k)-b^4*b1-b^3*b2-b^2*b3-b*b4);
b6=floor(b^6*x(k)-b^5*b1-b^4*b2-b^3*b3-b^2*b4-b*b5);
b7=floor(b^7*x(k)-b^6*b1-b^5*b2-b^4*b3-b^3*b4-b^2*b5-b*b6);
b8=floor(b^8*x(k)-b^7*b1-b^6*b2-b^5*b3-b^4*b4-b^3*b5-b^2*b6-b*b7);
b9=floor(b^9*x(k)-b^8*b1-b^7*b2-b^6*b3-b^5*b4-b^4*b5-b^3*b6-b^2*b7-b*b8);
b10=floor(b^10*x(k)-b^9*b1-b^8*b2-b^7*b3-b^6*b4-b^5*b5-b^4*b6-b^3*b7-b^2*b8-b*b9);
b11=floor(b^11*x(k)-b^10*b1-b^9*b2-b^8*b3-b^7*b4-b^6*b5-b^5*b6-b^4*b7-b^3*b8-b^2*b9-b*b10);
b12=floor(b^12*x(k)-b^11*b1-b^10*b2-b^9*b3-b^8*b4-b^7*b5-b^6*b6-b^5*b7-b^4*b8-b^3*b9-b^2*b10-b*b11);
b13=floor(b^13*x(k)-b^12*b1-b^11*b2-b^10*b3-b^9*b4-b^8*b5-b^7*b6-b^6*b7-b^5*b8-b^4*b9-b^3*b10-b^2*b11-b*b12);
b14=floor(b^14*x(k)-b^13*b1-b^12*b2-b^11*b3-b^10*b4-b^9*b5-b^8*b6-b^7*b7-b^6*b8-b^5*b9-b^4*b10-b^3*b11-b^2*b12-b*b13);
b15=floor(b^15*x(k)-b^14*b1-b^13*b2-b^12*b3-b^11*b4-b^10*b5-b^9*b6-b^8*b7-b^7*b8-b^6*b9-b^5*b10-b^4*b11-b^3*b12-b^2*b13-b*b14);

xb=[b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15];

%  Now change all '1'  to '2'

if b1==1 
    b1=2;
end
if b2==1
    b2=2;
end
if b3==1
    b3=2;
end
if b4==1
    b4=2;
end
if b5==1
    b5=2;
end
if b6==1
    b6=2;
end
if b7==1
    b7=2;
end
if b8==1
    b8=2;
end
if b9==1
    b9=2;
end
if b10==1
    b10=2;
end
if b11==1
    b11=2;
end
if b12==1
    b12=2;
end
if b13==1
    b13=2;
end
if b14==1
    b14=2;
end
if b15==1
    b15=2;
end

%  Now compute y

x(k)=b1/d + b2/(d^2) + b3/(d^3) + b4/(d^4) + b5/(d^5) + b6/(d^6) + ...
    b7/(d^7) + b8/(d^8) + b9/(d^9) + b10/(d^10) + b11/(d^11) + ...
    b12/(d^12) + b13/(d^13) + b14/(d^14) + b15/(d^15);

end


assignin('base','x',x);   %  this is the values of y