%% Sunpack.m % Soma_osf.m unpacks itself. fin = fopen("Soma_osf.m",'r'); mkdir("Soma"); L = fgetl(fin); while ~isequal(L,'%% EOF') F = ['Soma/' L(4:end)]; disp(F) fout = fopen(F,'w'); L = fgetl(fin); while length(L) < 2 || ~isequal(L(1:2),'%%') fprintf(fout,'%s\n',L); L = fgetl(fin); end fclose(fout); end fclose(fin); %% Soma_Solve.m function Soma_Solve(idx) if nargin == 0 somafigure % McKeeman, MATLAB Central File Exchange #26346, 2010. obj = somasol(); sols = obj.solutions(); set(gcf,'userdata',sols) else sols = get(gcf,'userdata'); idxstr = get(gca,'userdata'); idxstr.String = int2str(idx); somashow(colorvals(sols{idx}),idx); end end %% Spack.m function Spack fout = fopen('../Soma_osf.m','w'); d = dir; d = d(2:end); d = string({d.name})'; d(d=="Sunpack.m") = []; d(1) = "Sunpack.m"; for F = d' disp(F) fprintf(fout,'%%%% %s\n',F); fin = fopen(F,'r'); L = fgetl(fin); while ~isequal(L,-1) fprintf(fout,'%s\n',L); L = fgetl(fin); end fclose(fin); end fprintf(fout,'%%%% %s\n','EOF'); fclose(fout); end %% colorvals.m function cvals = colorvals(sol) % somashow-readable solution cvals = zeros(3,3,3); perm = [ 3 2 4 5 6 7 1]; for iv = 1:7 bits = sol(iv); for iu = 0:31 x = mod(iu,3); y = mod(floor(iu/3), 3); z = mod(floor(iu/9), 3); if (bitand(bits, bitshift(uint32(1), iu))) cvals(x+1,y+1,z+1) = perm(iv); end end end end %% somafigure.m function somafigure % Initialize figure window for Soma cube clf set(gcf,'name','Soma', ... 'numbertitle','off', ... 'color',0.75*[1 1 1]) uicontrol('Style','pushbutton', ... 'String','+', ... 'Units','normalized', ... 'Position',[.56 .93 .08 .04], ... 'Fontsize',12, ... 'Fontweight','bold', ... 'Callback',@somacb) uicontrol('Style','pushbutton', ... 'String','-', ... 'Units','normalized', ... 'Position',[.36 .93 .08 .04], ... 'Fontweight','bold', ... 'Fontsize',12, ... 'Callback',@somacb) idxstr = uicontrol( ... 'units','normalized', ... 'position',[.46 .93 .08 .04], ... 'style','pushbutton', ... 'string',0, ... 'Fontweight','bold', ... 'Fontsize',9, ... 'callback',@somacb); set(gca,'position',[0 0 1 1], ... 'userdata',idxstr, ... 'clipping','off') axis(14*[-1 1 -1 1 -1 1]) axis('square','off','vis3d') rotate3d on view(3) shg function somacb(arg,~) idxstr = get(gca,'userdata'); idx = str2double(idxstr.String); switch arg.String case '+', idx = mod(idx,240)+1; case '-', idx = mod(idx-2,240)+1; otherwise, idx = randi(240); end idxstr.String = num2str(idx); Soma_Solve(idx) end end %% somashow.m function somashow(cvals,idx) % somashow(cvals) Display solution colorvals. cla Q = Q0; F = [1 5 7 3 3 7 8 4 1 3 4 2 2 4 8 6 1 2 6 5 5 6 8 7]; C = get(0,'defaultaxescolororder'); idxstr = get(gca,'userdata'); if nargin < 2 idx = str2double(idxstr.String); idx = idx + 1; idxstr.String = int2str(idx); else idxstr.String = int2str(idx); end % Cube for j = 1:27 patch( ... Vertices = Q{j}, ... Faces = F, ... FaceColor = C(cvals(j),:)) end % Pieces for k = 1:7 switch k case 1, jk = [ 1 2 4]; % V case 2, jk = [ 1 2 3 4]; % L case 3, jk = [ 1 2 3 5]; % T case 4, jk = [ 1 2 5 6]; % Z case 5, jk = [ 1 2 5 14]; % S case 6, jk = [ 1 2 11 14]; % R case 7, jk = [ 1 2 4 10]; % Y end t = 2*pi*k/7; shift = 15*[cos(t) sin(t) .15]; for j = 1:length(jk) patch( ... Vertices = .75*Q{jk(j)} + shift, ... Faces = F, ... FaceColor = C(k,:)) end end drawnow function Q = Q0 % From Qube, the Digital Rubik's Cube. v = 0.90*qzero; Q = cell(3,3,3); for z = [-2 0 2] for y = [-2 0 2] for x = [-2 0 2] Q{x/2+2,y/2+2,z/2+2} = v + [x y z]; end end end end function q0 = qzero % Unit cubelet. q0 = [-1 -1 -1 -1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 -1 1 1 1 -1 1 1 1]; end end %% somasol.m function obj = somasol() % Bill McKeeman, MATLAB Central File Exchange #26346, 2010. % CBM, insert line 364, June 12, 2022. % % SOMASOL Compute all the unique solutions to the soma cube. % Examples: % >> obj = somasol(); % >> sols = obj.solutions(); % >> solct = numel(sols); % >> as_vec = obj.vecsol(sols{1}(1)); % it will be the T % >> as_bits = obj.bitsol(sols{1}(1)); % T as bit vector % >> howto = obj.layers(sols{1}); % how to build solution 1 % % bot mid top % TLL RVS RRS The first solution % TTL RVV YSS % TZL YZZ YYZ % % The soma cube puzzle consists of 7 pieces (named V L T Z R S Y) which can % be assembled into a 3x3x3 cube (consisting of 27 labeled cubelets). See % the soma demo in MATLAB(TM) for more details. % % Historical note: This computation was originally done in V4 MATLAB. The % generated tables were used in the MATLAB soma demo but the table % generating code itself was lost. I have recreated the table computation % here. In V4 on a mac it took several hours. In R2009b on a PC it took 3 % seconds. % % The mapping of cubelet labels to bit positions is rectangular, the lower % left corner being 0, increasing by 1 in the x direction, by 3 in the y % direction, and 9 in the z direction. Converting a label to base-3 % notation gives the xyz coordinates of the cubelet. %{ bottom middle top 6 7 8 151617 242526 y 3 4 5 121314 212223 0 1 2 91011 181920 x %} % somasol computes 240 solutions, each of which is a vector of 7 piece % positions. % % The computational strategy is to consider the cube as a set, represented % by the leftmost 27 1-bits in a uint32 value. A piece position is another % uint32 value (containing only 3 or 4 1-bits). Two piece positions collide % when their intersection (bitand) is non-zero. Any collection consisting % of one each of the seven pieces, the union (bitor) of which fills the 27 % bits is a solution. % % The first task is computing all 688 piece positions. The second task is % classifying them into mutually compatible subsets (theory needed here). % Ref: McKeeman, W.M., A Formal Model for Space-filling Puzzles, % Machine Intelligence vol.8, Elcock&Michie (ed.), Wiley (1977) pp 86-93 % Conway, Berlekamp and Guy, Winning Ways, vol. 2. See also % http://www.cs.dartmouth.edu/~mckeeman % % The third task is to remove equivalent solutions (rotation and % reflection) before doing the computation (special knowledge of the soma % cube needed here). Finally backtracking is applied to the remaining % search space. % 3x3x3 rotation permutation vectors. r = uint8([0 3 6]); r = [r+2 r+1 r]; r = [r r+9 r+18]; xy = r; r = uint8([6 7 8]); r = [r r+9 r+18]; r = [r r-3 r-6]; yz = r; r = uint8([18 9 0]); r = [r r+3 r+6]; r = [r r+1 r+2]; zx = r; % vefc categorized pieces positions and vefc categories [V,vs] = Vsets(); if size(vs,1) ~= 6; error('bad V'); end [L,ls] = Lsets(); if size(ls,1) ~= 4; error('bad L'); end %[T,ts] = Tsets(); % use T piece to remove rotations T = {v2b(uint8([0 1 2 4]))}; ts = [2 1 1 0]; [Z,zs] = Zsets(); if size(zs,1) ~= 2; error('bad Z'); end [R,rs] = Rsets(); if size(rs,1) ~= 4; error('bad R'); end [S,ss] = Ssets(); if size(ss,1) ~= 4; error('bad S'); end [Y,ys] = Ysets(); if size(ys,1) ~= 4; error('bad Y'); end % remove reflections using Y piece positions ed = uint8([0 3 6]); leftface = v2b([ed; ed+9; ed+18]); for ih = 1:numel(Y); Y{ih} = remove(Y{ih}, leftface); end % compute mutually consistent sets of piece positions compat = makeTrials(vs, ls, ts, zs, rs, ss, ys); % solution strategy: % place T in its fixed position first. % place V last (it is easiest). sols = {}; for iz=1:size(compat,1) idx = compat(iz,:); % idx order V L T Z S R Y trial = {T{1}, L{idx(2)}, Z{idx(4)}, ... R{idx(5)}, S{idx(6)}, Y{idx(7)}, V{idx(1)}}; sols = [sols, solve(trial)]; end obj = public(); return; % that's all folks % --- nested function definitions --- % the public interface function o = public() o = struct; o.solutions = sols; o.vecsol = @b2v; o.bitsol = @b2str; o.layers = @b2layers; end % all piece positions in terms of labels for V L R Z R S Y function res = Vvecs() p = uint8([0 1 3]); % init V p = [p; p+1]; p = [p; p+3]; p = [p; p+9; p+18]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Lvecs() p = uint8([0 1 2 3]); % init L p = [p; p+3]; p = [p; p+9; p+18]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Tvecs() p = uint8([0 1 2 4]); % init T p = [p; p+3]; p = [p; p+9; p+18]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Zvecs() p = uint8([0 1 4 5]); % init Z p = [p; p+3]; p = [p; p+9; p+18]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Rvecs() p = uint8([0 1 3 12]); % init R p = [p; p+1]; p = [p; p+3]; p = [p; p+9]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Svecs() p = uint8([0 1 3 10]); % init S p = [p; p+1]; p = [p; p+3]; p = [p; p+9]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end function res = Yvecs() p = uint8([0 1 3 9]); % init Y p = [p; p+1]; p = [p; p+3]; p = [p; p+9]; p = rXYZ(p); res = unique(sort(p,2), 'rows'); end % get all piece positions as bit vectors function bits = Vbits(); bits = bitRep(Vvecs()); end function bits = Lbits(); bits = bitRep(Lvecs()); end function bits = Tbits(); bits = bitRep(Tvecs()); end function bits = Zbits(); bits = bitRep(Zvecs()); end function bits = Rbits(); bits = bitRep(Rvecs()); end function bits = Sbits(); bits = bitRep(Svecs()); end function bits = Ybits(); bits = bitRep(Yvecs()); end % get signatures for each piece position function sigs = Vsigs(); sigs = makeSigs(Vbits()); end function sigs = Lsigs(); sigs = makeSigs(Lbits()); end function sigs = Tsigs(); sigs = makeSigs(Tbits()); end function sigs = Zsigs(); sigs = makeSigs(Zbits()); end function sigs = Rsigs(); sigs = makeSigs(Rbits()); end function sigs = Ssigs(); sigs = makeSigs(Sbits()); end function sigs = Ysigs(); sigs = makeSigs(Ybits()); end % separate piece positions by signature function [sets, s] = Vsets(); [sets, s] = makeSets(Vbits(), Vsigs()); end function [sets, s] = Lsets(); [sets, s] = makeSets(Lbits(), Lsigs()); end function [sets, s] = Tsets(); [sets, s] = makeSets(Tbits(), Tsigs()); end function [sets, s] = Zsets(); [sets, s] = makeSets(Zbits(), Zsigs()); end function [sets, s] = Rsets(); [sets, s] = makeSets(Rbits(), Rsigs()); end function [sets, s] = Ssets(); [sets, s] = makeSets(Sbits(), Ssigs()); end function [sets, s] = Ysets(); [sets, s] = makeSets(Ybits(), Ysigs()); end % make vector of signatures parallel to vector of piece positions function sigs = makeSigs(pp) vc = v2b(vCat()); ec = v2b(eCat()); fc = v2b(fCat()); cc = v2b(cCat()); sigs = zeros(size(pp,1), 4); for ig=1:size(pp,1) t(1) = bitsum(bitand(pp(ig,:), vc)); t(2) = bitsum(bitand(pp(ig,:), ec)); t(3) = bitsum(bitand(pp(ig,:), fc)); t(4) = bitsum(bitand(pp(ig,:), cc)); sigs(ig,:) = t; % signature for this piece position end end % given piece position & signature, separate by signature function [sets, s] = makeSets(pp, sigs) s = unique(sigs, 'rows'); sets = cell(size(s,1), 1); for jm=1:size(s,1); sets{jm} = uint32([]); end for im=1:size(pp,1) for jm=1:size(s,1) if all(s(jm,:)==sigs(im,:)) sets{jm} = [sets{jm}; pp(im)]; end end end end % solve integer linear equation function compat = makeTrials(v1,v2,v3,v4,v5,v6,v7) compat = zeros(0,7,'uint8'); for i1=1:size(v1,1) for i2=1:size(v2,1) for i3=1:size(v3,1) for i4=1:size(v4,1) for i5=1:size(v5,1) for i6=1:size(v6,1) for i7=1:size(v7,1) s = sum([v1(i1,:);v2(i2,:);v3(i3,:);v4(i4,:); v5(i5,:);v6(i6,:);v7(i7,:)]); % v e f c if all(s==[8 12 6 1]) % a solution compat = [compat; [i1,i2,i3,i4,i5,i6,i7]]; end end end end end end end end end % rotationally symmetric catgories function res = vCat() % vertex res = rXYZ(uint8(0)); % 0 2 6 8 18 20 24 26 end function res = eCat() % edge res = rXYZ(uint8(1)); % 1 3 5 7 9 11 15 17 19 21 23 25 end function res = fCat() % face res = rXYZ(uint8(4)); % 4 10 12 14 16 22 end function res = cCat() % center res = uint8(13); % 13 end % rotation actions function res = rXY(p) res = xy(p+1); % 90 degrees X into Y if size(p,2)==1; res = res'; end end function res = rYZ(p) res = yz(p+1); % 90 degrees Y into Z if size(p,2)==1; res = res'; end end function res = rZX(p) res = zx(p+1); % 90 degrees Z into X if size(p,2)==1; res = res'; end end function res = rXYZ(p) % all 24 rotations p = unique([p; rXY(p); rXY(rXY(p)); rXY(rXY(rXY(p)))], 'rows'); p = unique([p; rYZ(p); rYZ(rYZ(p)); rYZ(rYZ(rYZ(p)))], 'rows'); p = unique([p; rZX(p); rZX(rZX(p)); rZX(rZX(rZX(p)))], 'rows'); res = p; end % vector to bits function bits = v2b(vec) bits = uint32(0); for ib=1:numel(vec) bits = bitor(bits, bitshift(uint32(1),vec(ib))); end end % bits to vector function vec = b2v(bits) vec = zeros(1,0, 'uint8'); for iv=0:31 if (bitand(bits, bitshift(uint32(1), iv))) vec(end+1) = iv; end end end % for printing function str = b2str(bits) str = char(zeros(1,32)); for iv=0:31 if (bitand(bits, bitshift(uint32(1), iv))) str(iv+1) = '1'; else str(iv+1) = '0'; end end end % human-readable solution function str = b2layers(sol) str = char(zeros(3,3)); names = 'TLZSRYV'; % the order in the solve input for iv=1:numel(sol) bits = sol(iv); for iu=0:31 x = mod(iu,3); y = mod(floor(iu/3), 3); z = mod(floor(iu/9), 3); if (bitand(bits, bitshift(uint32(1), iu))) str(x+1,y+1,z+1) = names(iv); end end end end % add up bits function n = bitsum(bits) n = 0; for i=0:31 if bitand(bits, uint32(1)); n = n + 1; end bits = bitshift(bits, -1); end end % change vector sets into bit sets function bits = bitRep(matrix) bits = zeros(size(matrix, 1), 1, 'uint32'); for i=1:size(matrix,1) bits(i) = v2b(matrix(i,:)); end end % Q is a vector of piece positions % F is any part of the 3x3x3 cube. % Discard all Q that intersect F. function Q = remove(Q, F) for ij = size(Q,1):-1:1 % each piece if bitand(Q(ij), F) ~= 0 Q(ij) = []; % discard pieces that hit forbid end end end % given signature-compatible inputs, get solutions function sols = solve(ppsets) sols = {}; cube = v2b(0:26); recur(uint32(0), uint32([]), 1); % nothing so far... function recur(sofar, partial, n) if n > numel(ppsets) % cannot try more if sofar == cube % found solution sols = [sols, partial]; somashow(colorvals(partial)) % cbm, 2022 end return; end pps = ppsets{n}; for ir = 1:size(pps,1) pp = pps(ir); if bitand(sofar,pp) == 0 % no collisions %fprintf('%d %d %s\n', n, ir, b2str(sofar)) recur(bitor(sofar,pp), [partial, pp], n+1); end end end end end %% EOF