[equations,answers]=numbersgame(target,numbers,maxerror,allowfractions)

function [equations,answers]=numbersgame(target,numbers,maxerror,allowfractions)
%  function [equations,answers]=numbersgame(target,numbers,[maxerror],[allowfractions])
% solves the countdown numbers game
% e.g.a=numbersgame(715,[1 2 8 6 4 7]
%
%
if nargin<4
    allowfractions=0;
end
if nargin<3
    maxerror=10;
end
numbers=fliplr(sort(numbers(:)'));
vectors=numbers';
for ii=1:length(numbers);
    terms(ii,1)={num2str(numbers(ii))};
end

% now recursively expand to more vectors with fewer terms each by joining pairs of terms
for ii=1:length(numbers)-2;
    [vectors,terms]=expand(vectors,terms,allowfractions);
end

% break the final expansion up into 4 smaller pieces for memory reasons
st=floor(size(vectors,2)/4);
termsA=terms(:,1:st);
vectorsA=vectors(:,1:st);
termsB=terms(:,st+1:2*st);
vectorsB=vectors(:,st+1:2*st);
termsC=terms(:,2*st+1:3*st);
vectorsC=vectors(:,2*st+1:3*st);
terms=terms(:,3*st+1:end);
vectors=vectors(:,3*st+1:end);

% expand first quarter
[vectors,terms]=expand(vectors,terms,allowfractions);
equations={};
maxe=200000*eps;
while isempty(equations) & (maxe<=maxerror+100*eps)
    equations=terms(abs(target-vectors)<=maxe)';
    answers=vectors(abs(target-vectors)<=maxe);
    maxe=maxe+1;
end
clear terms vectors

% expand next quarter
[vectorsA,termsA]=expand(vectorsA,termsA,allowfractions);
equationsA={};
maxeA=200000*eps;
while isempty(equationsA) & (maxeA<=maxe-1+100*eps)
    equationsA=termsA(abs(target-vectorsA)<=maxeA)';
    answersA=vectorsA(abs(target-vectorsA)<=maxeA);
    maxeA=maxeA+1;
end
clear termsA vectorsA
if maxeA==maxe
    equations=unique([equations;equationsA]);
    answers=[answers answersA];
else
    equations=equationsA;
    answers=answersA;
end


% expand next quarter
[vectorsB,termsB]=expand(vectorsB,termsB,allowfractions);
equationsB={};
maxeB=200000*eps;
while isempty(equationsB) & (maxeB<=maxeA-1+100*eps)
    equationsB=termsB(abs(target-vectorsB)<=maxeB)';
    answersB=vectorsB(abs(target-vectorsB)<=maxeB);
    maxeB=maxeB+1;
end
clear termsB vectorsB
if maxeB==maxeA
    equations=unique([equations;equationsB]);
    answers=[answers answersB];
else
    equations=equationsB;
    answers=answersB;
end


% expand last quarter
[vectorsC,termsC]=expand(vectorsC,termsC,allowfractions);
equationsC={};
maxeC=200000*eps;
while isempty(equationsC) & (maxeC<=maxeB-1+100*eps)
    equationsC=termsC(abs(target-vectorsC)<=maxeC)';
    answersC=vectorsC(abs(target-vectorsC)<=maxeC);
    maxeC=maxeC+1;
end
if maxeC==maxeB
    equations=unique([equations;equationsC]);
    answers=[answers answersC];
else
    equations=equationsC;
    answers=answersC;
end



if ~isempty(answers) && all(abs(answers-answers(1))<200000*eps)
    answers=answers(1);                     % make only one answer, if all same
end 

return


%
% function which takes all possible pairs of terms in the 
% input matrix and combines in all possible way to get an
% more vectors each with one less term.
% 
function [vectorsout,termsout]=expand(vectorsin,termsin,allowfractions)
[nnum,nvec]=size(vectorsin);
newsize=nvec*(4*nnum*(nnum-1)/2+nnum);
vectorsout=zeros(nnum-1,newsize);
termsout=cell(size(vectorsout));
index=0;
pluses=repmat('+',1,nvec);
minuses=repmat('-',1,nvec);
lbrackets=repmat('(',1,nvec);
rbrackets=repmat(')',1,nvec);
multiplies=repmat('x',1,nvec);
divides=repmat('/',1,nvec);

    for i1=1:nnum
        n1=vectorsin(i1,:);                                             % first number of pair
        term1=termsin(i1,:);
        term1=char(term1{:}).';
        for i2=i1+1:nnum
            n2=vectorsin(i2,:);                                         % second number of pair
            term2=termsin(i2,:);
            term2=char(term2{:}).';
            othernums=vectorsin([1:i1-1,i1+1:i2-1,i2+1:nnum],:);
            otherterms=termsin([1:i1-1,i1+1:i2-1,i2+1:nnum],:);
            for operator='+-*/'
                switch operator
                    case '+'
                        answer=n1+n2;
                        term=cellstr([lbrackets;term1;pluses;term2;rbrackets]');
                        term=strrep(term,' ','')';
                        valid=true(1,size(term,2));
                    case '-'
                        answer=n1-n2;
                        term=cellstr([lbrackets;term1;minuses;term2;rbrackets]');
                        valid=(answer>0);
                        term=strrep(term(valid),' ','')';
                        answer=answer(valid);     % remove non integer results
                    case '*'
                        answer=n1.*n2;
                        term=cellstr([term1;multiplies;term2]');
                        valid=(n2>1);
                        term=strrep(term(valid),' ','')';
                        answer=answer(valid);     % remove non integer results
                    case '/'
                        answer=n1./n2;
                        term=cellstr([term1;divides;term2]');
                        valid=(allowfractions|((answer-floor(answer))==0))&(n2>1);
                        term=strrep(term(valid),' ','')';
                        answer=answer(valid);     % remove non integer results
                end     % end switch
                nvl=size(term,2);
                if nvl>0
                    vectorsout(:,index+1:index+nvl)=[answer; othernums(:,valid)];
                    termsout(:,index+1:index+nvl)=[term; otherterms(:,valid)];
                    index=index+nvl;
                end
            end         % next operator
        end             % next n2
        clear term2;
        othernums=vectorsin([1:i1-1,i1+1:nnum],:);
        otherterms=termsin([1:i1-1,i1+1:nnum],:);
        vectorsout(:,index+1:index+nvec)=othernums;
        termsout(:,index+1:index+nvec)=otherterms;
        index=index+nvec;
    end % next n1
clear term1 otherterms termsin;
% sort new numbers in each vector with biggest first
[vectorsout,isort]=sort(vectorsout(:,1:index),1,'descend');
% now copy all the terms with the vector sorting indices
terms=termsout(:,1:index);
for ii=1:index;
    terms(:,ii)=termsout(isort(:,ii),ii);
end
termcat=cell(size(terms,2),1);
for ii=1:size(terms,2);
    termcat(ii)={[terms{:,ii}]};                 % concatenate all terms in each vector to allow UNIQUE sets of terms to be found
end
[garbage,isort]=unique(termcat);
vectorsout=vectorsout(:,isort);
termsout=terms(:,isort);