function turns = snake()
% solves the wooden snake puzzle by brute force unraveling of states

global DEBUG_ON DIM CUBES DBG_SOL ALL_SOL;

% cube edge size and total cubicles
DIM =3;
CUBES = DIM*DIM*DIM;

% offsets of elbow cubicles (change of direction)
% should be increasing numbers
% elbows = [3 5 7 8]; this is a test flat side
% puzzle from this video www.youtube.com/watch?v=iTzVPgFjE9c
% the first cube is offset #1
elbows = [3 5 7 9 10 11 12 14 16 17 18 20 21 23 24 25]; % youtube
%elbows = [2 3 4 5 7 9 10 11 12 13 14 15 16 17 18 19 21 23 24 25 26]; % mine
NE = length(elbows);

if DEBUG_ON & ( ...
	elbows(end) >= CUBES | any(sort(elbows) ~= elbows) | ...
	length(unique(elbows)) < NE ) | elbows(1) <= 1
	error('elbows');
end
% elbows can be spaced at most 2 positions apart
if DEBUG_ON & ( ...
	elbows(1) > DIM | any(elbows(2:NE)-elbows(1:NE-1) >=DIM) )
	error('spacing');
end

% we have 3 dimensions X Y Z and a direction in each + or -
% each elbow can turn 4 ways, like the times of the clock 12, 3, 6 and 9
% turns are considered relative to the current orientation
XPLUS = 1; XMINU = 2;
YPLUS = 3; YMINU = 4;
ZPLUS = 5; ZMINU = 6;
T12 = 1; T3 = 2; T6 = 3; T9 = 4;

% when we turn an elbow, we switch orientation and direction
% row is the present orientation, column the turn applied
transit = [
ZPLUS YMINU ZMINU YPLUS
ZPLUS YPLUS ZMINU YMINU
ZPLUS XPLUS ZMINU XMINU
ZPLUS XMINU ZMINU XPLUS
YMINU XPLUS YPLUS XMINU
YMINU XMINU YPLUS XPLUS
];


% exhaustive enumeration of turns till we hit the solution
DBG_SOL = [];
% some solutions are hard to bend, so find all of them
ALL_SOL = [];
for a=1:4
    % first turn is always feasible, so try the next one immediately
	sol = recurse(a,elbows,transit);
	if length(sol), 
		break; 
	end;
end

if 0 % i don't think there's any chance to solve it backwards if no forward exists!
if 0==length(sol)
	% perhaps it's asymmetric and solvable only from the other side?
	e2 = sort((CUBES+1) - elbows);
	if any(e2 ~= elbows)
		disp('flipping snake and retrying...');
        DBG_SOL = [];
		for a=1:4
			sol = recurse(a,elbows,transit);
			if length(sol), 
				break; 
			end;
		end		
	end
    
    if 0==length(sol)
       	disp('no solution exists, best attempt was:');
        disp(DBG_SOL);
    end
end
end

if length(ALL_SOL)
    % pick a solution that's easy to solve
    % if this doesn't fold, pick the next one!
    sol = ALL_SOL(1,:);

    if DEBUG_ON
        xyz = 1;
        for a=1:length(sol)
            xyz = [xyz transit(xyz(end), sol(a))];
        end
        disp(xyz); % see turns in XYZ space
    end
    
    disp(sprintf('%d alternative folds exist', size(ALL_SOL,1)));
    
	% show clock turn equivalents
	disp('clock turns: ');
    turns = [];
    line = [];
	for i=1:length(sol)
		if sol(i)==1
			c=12;
		elseif sol(i)==2
			c=3;
		elseif sol(i)==3
			c=6;
		else
			c=9;
		end
		line = [line sprintf('%do ',c)];
        turns = [turns c];
	end
    disp(line);
    
    % on second thought, XYZ orientations are easier to fathom
    xyz = 1;
    for a=1:length(sol)
        xyz = [xyz transit(xyz(end), sol(a))];
    end
    disp('xyz arm orientation:'); % see turns in XYZ space
    line = [];
    for a=1:length(xyz)
        if xyz(a) < 3
            c = 'X';
        elseif xyz(a) < 5
            c = 'Y';
        else
            c = 'Z';
        end
 
        if mod(xyz(a),2)==1 
		    sgn = '+';
	    else 
		    sgn = '-';
	    end
        line = [line sprintf('%d:%c%c ',a,c,sgn)];
    end
    disp(line);    

else
    disp('no solution exists!');
end

return;


% ------- recursive kernel
function sol = recurse(partial,elbow,transit)

global DEBUG_ON ALL_SOL;
if DEBUG_ON & length(partial) >= length(elbow)
	error('recurse fault');
end

% try all turns and further recurse
for a=1:4
	snak = [partial a];
	yep = feasible(snak,elbow,transit);
	if yep
		if length(snak)==length(elbow)
			% sol = snak;
            ALL_SOL = [ALL_SOL; snak];
			%return; find more possibilities
            disp('.');
            if size(ALL_SOL,1) > 100
                return;
            end
        else
    		% go deeper
		    sol = recurse(snak,elbow,transit);
		    if length(sol)
                if DEBUG_ON 
                	error('no solution should return here');
                end
			    break;
		    end
        end
	end
end

return;


% ------- is this snake feasible?
function yep = feasible(snak,elbow,transit)
global DEBUG_ON DIM CUBES DBG_SOL;

if DEBUG_ON & length(snak) > length(elbow)
	error('long snake');
end

elb2 = [1 elbow CUBES];
yep = 1;

% there are 2 fouls to consider, any dimension growing more than 3
% or overlapping cubicles (turn upon itself)

% the first arm is placed in the X+ direction
dir = 1; sgn = 1; % XPLUS
xyz = 1;

mn = ones(1,3)*DIM;
mx = ones(1,3)*DIM;
pt = ones(1,3)*DIM; % current position
% whole cube is represented expanded to cater for non-central arrangement
cube = zeros(2*DIM-1, 2*DIM-1, 2*DIM-1);
cube(DIM, DIM, DIM)=1;

for a=1:length(snak)+1
	% examine current arm (minus the elbow cubicle)
	len = elb2(a+1)-elb2(a);  
	p2 = pt(dir) + sgn*len;
	if mn(dir) > p2
		mn(dir) = p2;
	end
	if mx(dir) < p2
		mx(dir) = p2;
	end

	if mx(dir) - mn(dir) >= DIM
		% already out of bounds, can't fold into cube
		yep = 0;
		break;
	end

	% current position is already 1
	if DEBUG_ON & cube(pt(1), pt(2), pt(3)) ~= 1
		error('cube');
	end
	
	% near impossible to index in 3 dimensions using a direction only
	x = 1==dir;
	y = 2==dir;
	z = 3==dir;
	for i=1:len
		if cube(pt(1)+i*sgn*x, pt(2)+i*sgn*y, pt(3)+i*sgn*z)
			yep = 0;
			break; % spot already taken, can't fold!
		else
			cube(pt(1)+i*sgn*x, pt(2)+i*sgn*y, pt(3)+i*sgn*z) = 1;
		end
	end

	if 0==yep
		break;
	end

	pt(dir) = p2;

	if length(snak)+1 == a
		break;
	end

    % orientation of the next arm
	xyz = transit(xyz, snak(a));
	dir = floor((xyz+1) / 2);
	if DEBUG_ON & (dir <= 0 | dir > 6)
		error('direction mapping');
	end
	if mod(xyz,2)==1 
		sgn = 1;
	else 
		sgn = -1;
	end
end

if yep & length(snak)==length(elbow)
	yep = 2; % solution!
end

if DEBUG_ON & yep & length(snak) > length(DBG_SOL)
    DBG_SOL = snak; % keep best trial so far
end

return;