[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
for ii=1:length(numbers)-1;
    [vectors,terms]=expand(vectors,terms,allowfractions);
end
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
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=size(vectorsin,1);
newsize=size(vectorsin,2)*(4*nnum*(nnum-1)/2+nnum);
vectorsout=zeros(nnum-1,newsize);
termsout=cell(size(vectorsout));
index=1;
for ii=1:size(vectorsin,2);
    numbers=vectorsin(:,ii);
    terms=termsin(:,ii);
    for i1=1:nnum
        n1=numbers(i1);                                             % first number of pair
        term1=terms{i1};
        for i2=i1+1:nnum
            n2=numbers(i2);                                         % second number of pair
            term2=terms{i2};
            othernums=numbers([1:i1-1,i1+1:i2-1,i2+1:nnum]);
            otherterms=terms([1:i1-1,i1+1:i2-1,i2+1:nnum]);
            for operator='+-*/'
                switch operator
                    case '+'
                        answer=n1+n2;
                        term=sprintf('(%s+%s)',term1,term2);
                        vectorsout(:,index)=[answer; othernums];
                        termsout(:,index)=[{term}; otherterms];
                        index=index+1;
                    case '-'
                        if (n1>n2)           % ignore case n2=n1
                            answer=n1-n2;
                            term=sprintf('(%s-%s)',term1,term2);
                            vectorsout(:,index)=[answer; othernums];
                            termsout(:,index)=[{term}; otherterms];
                            index=index+1;
                        end
                    case '*'
                        if (n2>1)           % ignore trivial case of multiply by 1
                            answer=n1*n2;
                            term=sprintf('%sx%s',term1,term2);
                            vectorsout(:,index)=[answer; othernums];
                            termsout(:,index)=[{term}; otherterms];
                            index=index+1;
                        end
                    case '/'
                        if (n2>1)           % ignore trivial case of divide by 1
                            if mod(n1,n2)==0
                                    answer=round(n1/n2);   
                                    term=sprintf('%s/%s',term1,term2);
                                    vectorsout(:,index)=[answer; othernums];
                                    termsout(:,index)=[{term}; otherterms];
                                    index=index+1;
                            elseif allowfractions  % ignore cases where divide result is non interger
                                    answer=round(n1/n2);   
                                    term=sprintf('(%s/%s)',term1,term2);
                                    vectorsout(:,index)=[answer; othernums];
                                    termsout(:,index)=[{term}; otherterms];
                                    index=index+1;
                            end
                        end
                end % end switch
            end % next operator
        end % next n2
        othernums=numbers([1:i1-1,i1+1:nnum]);
        otherterms=terms([1:i1-1,i1+1:nnum]);
        vectorsout(:,index)=othernums;
        termsout(:,index)=otherterms;
        index=index+1;
    end % next n1
end % next vector to expand
[vectorsout,isort]=sort(vectorsout(:,1:index-1),1,'descend');     % sort new numbers in each vector with biggest first
% now copy all the terms with the vector sorting indices
terms=termsout(:,1:index-1);
for ii=1:index-1;
    terms(:,ii)=termsout(isort(:,ii),ii);
end
termcat=cell(index-1,1);
for ii=1:index-1
       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);