%{
%--- S U D O K U   L Ö S E R ---
%v.1.1
%by Peter Jüstel
%
%Eingabe: 9x9 Sudoku Matrix(Zahlen von 1 bis 9). Nullen symbolisieren leere Felder.
%
%Dieses Programm ist der Versuch, den Lösungsweg, den ein Mensch bei einem Sudoku geht mit der Maschine nachzuvollziehen.
%}
function ergebnis = SudokuFast(A)

tic
%{
%Keep in mind that tic and toc measure overall elapsed time. 
%Make sure that no other applications are running in the background 
%on your system that could affect the timing of your MATLAB programs.
%}

global fZ
global fS
global fM
global C
global bf
global firstrun
global B
global smv
global flago
global err1
global fundcounter_bf
global washere

washere = 0;
err1 = 0;           % sollte das Sudoku irgendwelche unmöglichkeiten enthalten ändert sich dieser wert
flago = 0;
smv = [0 1];        % ist ne Abbruchbedingung für die Rekursion(siehe f_finishcheck)
B = A;

C = num2cell(A);    % Erfassen der Matrix
fZ = cell(9,1);     % Feststehende Zahlen in den Zeilen (enthält später Vektoren)
fS = cell(1,9);     % Feststehende Zahlen in den Spalten (enthält später Vektoren)
fM = cell(3,3);     % Feststehende Zahlen in den Untermatrizen (enthält später Vektoren)
firstrun = 1;

bf = cell(1,9);


%% -1- Befehlskette zum lösen
    f_loose
    function f_loose
        if firstrun == 1
            f_hilfscell(B)
            firstrun = 0;
        end

        C = f_possum(C);
        washere = 0;
        if err1 == 1
            disp('error')
            return
        end
        
        bf = f_bitfields(C);
        
        for sb = 1:9        % durchlaufen der Bitfields für die Zahlen von 1 bis 9
            for sc = 1:9    % durchlaufen der Zeilen/Spalten von 1 bis 9
                C = f_check(C,sb,sc);  % prüfen von Zeilen und Spalten auf Singletons
            end
            C = f_check_m(C,sb);   % prüfen der Minimatrizen auf Singletons
        end
        
        f_ausgabe
    end

%% -2- Speichern der vorhandenen Zahlen in Arrays

    function f_hilfscell(G)
        
        kM1=1; kM2=1; kM3=1; kM4=1; kM5=1; kM6=1; kM7=1; kM8=1; kM9=1; %Längenvariablen reset
        
        for ii = 1:9        %"Zeilen" durchlaufen(bei allen *Z), "Spalten" durchlaufen(bei allen *S)
            kZ=length(fZ{ii,1})+1;            %k zurück setzen beim wechseln der Zeile/Spalte, ohne eventuell vorhandene Werte zu verlieren
            kS=length(fS{1,ii})+1;
        
            for jj = 1:9    %"Spalten" durchlaufen(bei allen *Z), "Zeilen" durchlaufen(bei allen *S)
        
                %Sammeln der feststehenden Zahlen in der Zeile    
                if G(ii,jj) ~= 0
                    fZ{ii,1}(kZ)= B(ii,jj);         %Speichere alle gefunden Zahlen der Zeile "ii" in der zugehörigen Zelle "ii" von fZ(und an der Stelle kz des enthaltenenen Vektors).      
                    kZ = kZ + 1;                    %anwählen des nächsten Elements im Vektor der Zeile(fZ)

                    %Sammeln der feststehenden Zahlen in der Untermatrix
                    switch ii && jj
                        case ii<4 && jj<4               %links oben
                            fM{1,1}(kM1) = G(ii,jj);
                            kM1 = kM1 + 1;
                        case ii<4 && jj>6               %rechts oben
                            fM{1,3}(kM3) = G(ii,jj);
                            kM3 = kM3 + 1;
                        case ii>6 && jj<4               %links unten
                            fM{3,1}(kM7) = G(ii,jj);
                            kM7 = kM7 + 1;
                        case ii>6 && jj>6               %rechts unten
                            fM{3,3}(kM9) = G(ii,jj);
                            kM9 = kM9 + 1;
        
                        case ii<4 && jj>3 && jj<7       %oben mitte
                            fM{1,2}(kM2) = G(ii,jj);
                            kM2 = kM2 + 1;
                        case ii>3 && ii<7 && jj<4       %links mitte
                            fM{2,1}(kM4) = G(ii,jj);
                            kM4 = kM4 + 1;
                        case ii>3 && ii<7 && jj>6       %rechts mitte
                            fM{2,3}(kM6) = G(ii,jj);
                            kM6 = kM6 + 1;
                        case ii>6 && jj>3 && jj<7       %unten mitte
                            fM{3,2}(kM8) = G(ii,jj);
                            kM8 = kM8 + 1;
        
                        case ii>3 && ii<7 && jj>3 && jj<7   % mitte
                            fM{2,2}(kM5) = G(ii,jj);
                            kM5 = kM5 + 1;
                    end
                    
                end
        
                %Sammeln der feststehenden Zahlen in der Spalte
                if G(jj,ii) ~= 0        %Das selbe wie für die Zeilen, nur für die Spalten.
                    fS{1,ii}(kS) = G(jj,ii);
                    kS = kS+1;
                end      
            end
        end
 %{
%Bedeutend kürzere/variablere version von dem switch(s.o.), braucht aber auch nen bisschen länger
 %       for lakai = [1 4 7]        
 %           for lakaj = [1 4 7]
 %               for hai = ((lakai-1)+(1:3))  
 %                   for haj = ((lakaj-1)+(1:3))
 %                       if B(hai,haj) ~= 0
 %                           ix = ceil(lakai/3);
 %                           iy = ceil(lakaj/3);
 %                           hak = length(fM{ix,iy});
 %                           fM{ix,iy}(hak+1) = B(hai,haj);
 %                       end
 %                   end
 %               end
 %           end
 %       end
 %}           
    end

%% -3- Zusammensuchen der Wahrscheinlichen Zahlen für ein Feld

    function ergo = f_possum(G)
        fundcounter = 1;
        if washere == 0
        while fundcounter > 0 % Anzahl durchläufe zum suchen der Möglichkeiten
            fundcounter = 0;        % die schleife wird nicht mehr durchlaufen, solange nicht etwas neues gefunden wird.
            for mi = 1:9
                for mj = 1:9
                    if (G{mi,mj} == 0) | (length(G{mi,mj}) > 1)     %#ok<OR2> suchen nach Arbeit(leeres oder nicht bestimmtes Feld)
                        faz = fZ{mi,1};                             %Extrahieren der benötigten Vektoren aus den Hilfs Arrays
                        fas = fS{1,mj};
                        fam = fM{ceil(mi/3),ceil(mj/3)};
                        fubar = f_union(fam, f_union(fas,faz));     % bilden der Vereinigungsmenge aller bekannten Zahlen, die Einfluss auf das Feld haben
                        pvekt = setxor([1 2 3 4 5 6 7 8 9],fubar);   % bilden des Negativs der Vereinigungsmenge/bilden des Möglichkeitsvektors für das Feld
                        if isempty(pvekt) == true
                            err1 = 1;
                            ergo = G;
                            washere = 1;
                            return
                        else   
                        G{mi,mj} = pvekt;                            % Rückspeichern in der Möglichkeits-cell C
                        if length(pvekt) == 1                        % Speichern von Funden
                            fZ{mi,1}((length(faz)+1))= pvekt;
                            fS{1,mj}((length(fas)+1))= pvekt;
                            fM{ceil(mi/3),ceil(mj/3)}((length(fam)+1)) = pvekt;
                            B(mi,mj) = pvekt(1);
                            fundcounter = fundcounter + 1;          % aktivieren der Schleife, weil neue informationen existieren
                            fundcounter_bf = fundcounter_bf + 1;
                        end
                        end
                    end
                end
            end
        end
        end
        ergo = G;
    end
%% -4- Bitfields

    function ergo = f_bitfields(G)% erstellen der bitfields
        ef = cell(1,9);
        for efi = 1:9
            ef{efi} = zeros(9); % Möglichkeitsmatrizen/bitfields werden angelegt. Für jede Zahl eine.
        end
        for ki = 1:9
            for kj = 1:9
                l = length(G{ki,kj});
                if l > 1
                for kk = 1:l
                    gv = G{ki,kj}(kk);
                    switch gv
                        case 1
                            ef{1}(ki,kj)=1;
                        case 2
                            ef{2}(ki,kj)=1;
                        case 3
                            ef{3}(ki,kj)=1;
                        case 4
                            ef{4}(ki,kj)=1;
                        case 5
                            ef{5}(ki,kj)=1;
                        case 6
                            ef{6}(ki,kj)=1;
                        case 7
                            ef{7}(ki,kj)=1;
                        case 8
                            ef{8}(ki,kj)=1;
                        case 9
                            ef{9}(ki,kj)=1;
                    end
                end
                end
            end
        end
        ergo = ef;
    end
%% -5- Singletons

    function ergo = f_check(G,sb,sc)
        smz = 0;    % immer wenn die funktion aufgerufen wird, summiere ich ja eine neue Gegend. Deshalb müssen die variablen zurückggesetzt werden.
        sms = 0;
        
        for si = 1:9
            smz = smz + (bf{sb}(sc,si)); % summieren einer Zeile(sc) im Bitfield(sb)
            sms = sms + (bf{sb}(si,sc)); % summieren einer Spalte(sc) im Bitfield(sb)
        end
        
        % prüfen der Zeile auf einsame
        
        if smz == 1                     % wenn in der Zeile ein einsamer Einser gefunden wurde
            for sj = 1:9                % durchlaufe die Zeile
                if bf{sb}(sc,sj) == 1   % und finde den Ort des Einsamen
                    G{sc,sj} = sb;      % trage dann an dem zugehörigen Ort in G,B,f* die passende Zahl sb ein.
                    B(sc,sj) = sb;
                    
                    fZ{sc,1}((length(fZ{sc,1})+1))= sb;
                    fS{1,sj}((length(fS{1,sj})+1))= sb;
                    x1 = ceil(sc/3);
                    x2 = ceil(sj/3);
                    fM{x1,x2}(length(fM{x1,x2})+1) = sb;
                    
                    bf{sb}(sc,:) = 0;   % löschen der anderen möglichketen, die sich mit dem Fund schneiden
                    bf{sb}(:,sj) = 0;
                    for fi = 1:9
                        bf{fi}(sc,sj) = 0;
                    end
                    bf{sb} = f_clearmm(bf,sb,sc,sj);    %clear der Minimatrix(in bf{sb}) um punkt sc,sj
                    fundcounter_bf = fundcounter_bf + 1;
                    G = f_check(G,sb,sc);      % nach refresh rerun
                end
            end
        end
        % prüfen der Spalte auf einsame(fast das selbe wie for Zeilen,s.o.)
        if sms == 1
            for sj = 1:9
                if bf{sb}(sj,sc) == 1
                    G{sj,sc} = sb;
                    B(sj,sc) = sb;
                    
                    fZ{sj,1}((length(fZ{sj,1})+1))= sb;
                    fS{1,sc}((length(fS{1,sc})+1))= sb;
                    x1 = ceil(sj/3);
                    x2 = ceil(sc/3);
                    fM{x1,x2}(length(fM{x1,x2})+1) = sb;
                    
                    bf{sb}(sj,:) = 0;
                    bf{sb}(:,sc) = 0;
                    for fi = 1:9
                        bf{fi}(sj,sc) = 0;
                    end
                    bf{sb} = f_clearmm(bf,sb,sj,sc);
                    fundcounter_bf = fundcounter_bf + 1;
                    G = f_check(G,sb,sc);
                end
            end 
        end
        ergo = G;
    end
%% -6- Check für Minimatrizen
    function ergo = f_check_m(G,sb)
        
        for glai = 0:3:6        % diese beiden Laufvariablen realisieren die Grenzen zwischen den Minimatrizen
            for glaj = 0:3:6
                smm = 0;
                for mai = (glai+(1:3))  % diese beiden Laufvariablen durchlaufen die angewählte Minimatrix
                    for maj = (glaj+(1:3))
                        smm = smm +(bf{sb}(mai,maj)); % sucht nach einsamen in den Minimatrizen
                    end
                end
                
                if smm == 1
                    for mbi = (glai +(1:3))
                        for mbj = (glaj +(1:3))
                            if bf{sb}(mbi,mbj) == 1 % wenn der Ort mit der Einsamen gefunden ist
                            G{mbi,mbj} = sb;        % setze die aktuelle Bitfieldzahl an den gefundenen Ort.
                            B(mbi,mbj) = sb;
                            
                            fZ{mbi,1}((length(fZ{mbi,1})+1))= sb;   % Eintragen in die f's
                            fS{1,mbj}((length(fS{1,mbj})+1))= sb;
                            x1 = ceil(mbi/3);
                            x2 = ceil(mbj/3);
                            fM{x1,x2}(length(fM{x1,x2})+1) = sb;
                            
                            bf{sb}(mbi,:) = 0;      % clear der Spalten und Reihen
                            bf{sb}(:,mbj) = 0;
                            for fi = 1:9
                                bf{fi}(mbi,mbj) = 0;
                            end
                            bf{sb} = f_clearmm(bf,sb,mbi,mbj);  %clear der Minimatrix um Punkt mbi,mbj
                            fundcounter_bf = fundcounter_bf + 1;
                            G = f_check_m(G,sb);
                            end
                        end
                    end
                end
            end
        end
        ergo = G;
    end

%% -8- Minifunktionen


%
% Da die Funktion union sehr oft angewandt werden, jedoch
% ziemlich klobig sind hab ich maßgeschneiderte Versionen von ihr
% geschrieben, die zwar nicht sicher ist, aber schnell.
%
    function ergo = f_union(d,g)
        ergo = unique([d,g]);       % so einfach kann das sein...
    end


    function ergo = f_finishcheck(G) % es wird geprüft, ob in einer Zelle von C auch nurnoch eine sichere Zahl steht.
        summe = 0;
        for ai = 1:9
            for aj = 1:9
                if G{ai,aj} ~= 0        % keine leeren Felder kommen in die Kalkulation
                    summe = summe + length(G{ai,aj});   % addieren der längen der Wahrscheinlichkeitsvektoren
                end
            end
        end

        smv(1) = smv(2);
        smv(2) = summe; 
        if summe == 81    % wenn überall nurnoch feststehende Zahlen stehen, und die Zeilen nurnoch einsame enthalten
            summe2 = 0;
            E = f_cell2num(G);
            for stri = 1:9
                a = E(stri,:);
                b = (E(:,stri))';
                strzeile = num2str(sort(a));
                strspalte = num2str(sort(b)); %#ok<UDIM>
                if strcmp(strzeile, '1  2  3  4  5  6  7  8  9') == 1 && strcmp(strspalte, '1  2  3  4  5  6  7  8  9') == 1
                    summe2 = summe2 + 1;
                end
            end
            if summe2 == 9
                ergo = true;
            else
                ergo = false;
            end
        else
            ergo = false; % andernfalls ist es alles für die katz
        end

    end


    function ergo = f_cell2num(G) % rückkonvertierung des gelösten Cellarrays in eine Matrix
        ergo = zeros(9);
        for ci = 1:9
            for cj = 1:9
                ergo(ci,cj) = G{ci,cj}(1);
            end
        end
    end

    
    function ergo = f_clearmm(ef,sb,y,x)      % clearen einer Minimatrix um (x,y) innerhalb des Bitfields der Zahl sb
        buffery = 0;        %links oben
        bufferx = 0;
        switch y && x
            case y > 6 && x < 4     % links unten
                buffery = 6;
            case y < 4 && x > 6     % rechts oben
                bufferx = 6;
            case y > 6 && x > 6     % rechts unten
                buffery = 6;
                bufferx = 6;
            
            case y > 3 && y < 7 && x < 4    % links mitte
                buffery = 3;
            case y > 3 && y < 7 && x > 6    % rechts mitte
                buffery = 3;
                bufferx = 6;
            case y < 4 && x > 3 && x < 7    % oben mitte
                bufferx = 3;
            case y > 6 && x > 3 && x < 7    % unten mitte
                bufferx = 3;
                buffery = 6;
                
            case y > 3 && y < 7 && x > 3 && x < 7   % mitte
                buffery = 3;
                bufferx = 3; 
        end
        for mmi = 1:3
            for mmj = 1:3
                ef{sb}((buffery + mmi), (bufferx + mmj)) = 0;
            end
        end
        ergo = ef{sb};
    end
%% -9- Try and error 
    function ergo = f_try(G)

        if flago == 0       % nur beim ersten mal ausführen von try
            if B == zeros(9)
                ergo = num2cell([2 6 9 3 8 5 1 4 7; 1 8 4 9 2 7 5 6 3; 3 7 5 1 6 4 2 9 8; 9 4 3 2 7 6 8 1 5; 8 2 1 5 4 3 9 7 6; 7 5 6 8 9 1 4 3 2; 6 1 2 7 5 9 3 8 4; 4 9 8 6 3 2 7 5 1; 5 3 7 4 1 8 6 2 9]);
                return
            else
                summae = 0;
                for poi = 1:9
                    summae = [summae, fZ{poi,1}]; %#ok<AGROW>
                end
                lsummae = length(summae);        % so bekommt man den wert, wieviele Zahlen tatsächlich im Sudoku stehen
                if lsummae < 22         % an dieser stelle könnten doubletons implementiert werden
                    disp('Es ist ein Sudoku LÖSER, kein Sudoku GENERATOR. Es kann mit den implementierten Methoden nicht weiter gelöst werden.')
                    ergo = num2cell(B);
                    return
                else
                    flago = 1;
                    G = f_try(G);   % zu brute force übergehen
                end
            end
        end


        bkpG = G;       %Backups erstellen
        bkpZ = fZ;
        bkpS = fS;
        bkpM = fM;
        bkpB = B;
        bkpbf = bf;
        
        flag = 2;
        counter = 1;
        while counter < 501 && flag < 5       % Er fängt bei den zweier vektoren an, und zählt hoch bis zu den neunern. Außerdem durchläuft er maximal 500 versuche.
            try
                for ti = 1:9
                    for tj = 1:9
                        if length(G{ti,tj}) == flag
                            zufall = ceil(flag*rand(1));     %Er wählt zufällig einen der möglichen aus dem vektor
    
                            G{ti,tj} = bkpG{ti,tj}(zufall);     % speichert ihn überall ab
                            B(ti,tj) = G{ti,tj};
                            fZ{ti,1}((length(fZ{ti,1})+1))= G{ti,tj};
                            fS{1,tj}((length(fS{1,tj})+1))= G{ti,tj};
                            x1 = ceil(ti/3);
                            x2 = ceil(tj/3);
                            fM{x1,x2}(length(fM{x1,x2})+1) = G{ti,tj};
                            
                            fundcounter_bf = 1;
                            while fundcounter_bf > 0    % solange er was neues findet mit der Zufallszahl, macht er weiter
                                fundcounter_bf = 0;
                                G = f_possum(G);
                                washere = 0;
                                if err1 == 1        %wird gestoppt, wenn error auftritt
                                    G = bkpG;   % Geht zurück zur Wurzel des Möglichkeitenbaumes
                                    fZ = bkpZ;
                                    fS = bkpS;
                                    fM = bkpM;
                                    B = bkpB;
                                    bf = bkpbf;
                                    err1 = 0;
                                    counter = counter + 1;
                                    error(warnix)       % springt aus den schleifen zu catch
                                else                %mach weiter, wenn alles ok scheint
                                    durchgang = 0;
                                    bf = f_bitfields(G);
                                    while durchgang < 10
                                        for sb = 1:9        % durchlaufen der Bitfields für die Zahlen von 1 bis 9
                                            for sc = 1:9    % durchlaufen der Zeilen/Spalten von 1 bis 9
                                                G = f_check(G,sb,sc);  % prüfen von Zeilen und Spalten auf Singletons
                                            end
                                            G = f_check_m(G,sb);   % prüfen der Minimatrizen auf Singletons
                                        end
                                        durchgang = durchgang +1;
                                    end
                                end
                            end
                        end
                    end
                end
    
                flag = flag +1;
            catch warnix
                clear warnix
            end
        end
        
        if f_finishcheck(G) == true
            if washere == 0     % sonst zeigt er es 2 mal an
                disp(['-> brute force beendet nach ' num2str(counter) ' Versuchen.' ])
                washere = 1;
            end
            ergo = G;
        else
            if washere == 0
                disp('Sudoku weder mit Algorithmus, noch mit "brute force" gelöst.')
                washere = 1;
            end
            ergo = num2cell(zeros(9));
        end
    end
%% -10- Ausgabe
    function f_ausgabe
        a = f_finishcheck(C);
        if a == true                    % wenn das Sudoku gelöst ist,
            ergebnis = f_cell2num(C);   % wird es zurückkonvertiert in eine Matrix, und als ausgabewert festgelegt.
            t = toc
        else
            if smv(1) == smv(2)    % Wenn neue Zahlen dazu gekommen sind, ändert sich die summe. er probiert solange weiter, bis keine neuen zahlen mehr gefunden werden. 
                disp(' starte brute force . . .')
                ergebnis = f_cell2num(f_try(C)); 
                t = toc
            else
                f_loose
            end
        end
    end
end