Create expression trees from given sets of numbers and operations and find those that evaluate to a target number in Mathematica 8 or above

Question

Given a set of numbers and a set of binary operations, what is the fastest way to create random expression trees or exhaustively check every possible combination in Mathematica?

What I am trying to solve is given:

numbers={25,50,75,100,3,6}              (* each can ONLY be used ONCE  *)
operators={Plus,Subtract,Times,Divide}  (* each can be used repeatedly *)
target=99

find expression trees that would evaluate to target.

I have two solutions whose performances I give for the case where expression trees contain exactly 4 of the numbers and 3 operators:

1. random sample & choice: ~25K trees / second
2. exhaustive scan: 806400 trees in ~2.15 seconds

(timed on a laptop with: Intel(R) Core(TM)2 Duo CPU T9300 @ 2.50GHz, 3GB ram, no parallelization used yet but would be most welcome in answers)

My notebooks are a bit messy at the moment. So I would first love to pose the question and hope for original ideas and answers while I clean up my code for sharing.

Largest possible case is where every expression tree uses up all the (6) numbers and 'Length[numbers]-1' (5) operators.

Performance of methods in the largest case is:

1. random sample & choice: ~21K trees / second
2. exhaustive scan: 23224320 trees in ~100 seconds

Also I am using Mathematica 8.0.1 so I am more than all ears if there are any ways to do it in OpenCL or using compiled functions wiht CompilationTarget->"C", etc.

Answer 1

OK, this is not elegant or fast, and it's buggy, but it works (sometimes). It uses a monte carlo method, implementing the metropolis algorithm for a weight function that I (arbitrarily) selected just to see if this would work. This was some time ago for a similar problem; I suppose my mathematica skills have improved as it looks ugly now, but I have no time to fix it at the moment.

Execute this (it looks more reasonable when you paste it into a notebook):

ClearAll[swap];
swap[lst_, {p1_, p2_}] := 
 ReplacePart[
  lst, {p1 \[Rule] lst\[LeftDoubleBracket]p2\[RightDoubleBracket], 
   p2 \[Rule] lst\[LeftDoubleBracket]p1\[RightDoubleBracket]}]

ClearAll[evalops];
(*first element of opslst is Identity*)

evalops[opslst_, ord_, nums_] := 
 Module[{curval}, curval = First@nums;
  Do[curval = 
    opslst\[LeftDoubleBracket]p\[RightDoubleBracket][curval, 
     nums\[LeftDoubleBracket]ord\[LeftDoubleBracket]p\
\[RightDoubleBracket]\[RightDoubleBracket]], {p, 2, Length@nums}];
  curval]

ClearAll[randomizeOrder];
randomizeOrder[ordlst_] := 
 swap[ordlst, RandomInteger[{1, Length@ordlst}, 2]]

ClearAll[randomizeOps];
(*never touch the first element*)

randomizeOps[oplst_, allowedOps_] := 
 ReplacePart[
  oplst, {RandomInteger[{2, Length@oplst}] \[Rule] RandomChoice[ops]}]

ClearAll[takeMCstep];
takeMCstep[goal_, opslst_, ord_, nums_, allowedops_] := 
 Module[{curres, newres, newops, neword, p}, 
  curres = evalops[opslst, ord, nums];
  newops = randomizeOps[opslst, allowedops];
  neword = randomizeOrder[ord];
  newres = evalops[newops, neword, nums];
  Switch[Abs[newres - goal], 
   0, {newops, 
    neword}, _, (p = Abs[curres - goal]/Abs[newres - goal];
    If[RandomReal[] < p, {newops, neword}, {opslst, ord}])]]

then to solve your actual problem, do

ops = {Times, Plus, Subtract, Divide}
nums = {25, 50, 75, 100, 3, 6}
ord = Range[Length@nums]
(*the first element is identity to simplify the logic later*)
oplist = {Identity}~Join~RandomChoice[ops, Length@nums - 1]
out = NestList[
  takeMCstep[
    99, #\[LeftDoubleBracket]1\[RightDoubleBracket], #\
\[LeftDoubleBracket]2\[RightDoubleBracket], nums, ops] &, {oplist, 
   ord}, 10000]

and then to see that it worked,

ev = Map[evalops[#\[LeftDoubleBracket]1\[RightDoubleBracket], #\
\[LeftDoubleBracket]2\[RightDoubleBracket], nums] &, out];
ev // Last // N
ev // ListPlot[#, PlotMarkers \[Rule] None] &

giving

thus, it obtained the correct order of operators and numbers after around 2000 tries.

As I said, it's ugly, inefficient, and badly programmed as it was a quick-and-dirty adaptation of a quick-and-dirty hack. If you're interested I can clean up and explain the code.

Answer 2

This was a fun question. Here's my full solution:

ExprEval[nums_, ops_] := Fold[
  #2[[1]][#1, #2[[2]]] &,
  First@nums,
  Transpose[{ops, Rest@nums}]]

SymbolicEval[nums_, ops_] := ExprEval[nums, ToString /@ ops]

GetExpression[nums_, ops_, target_] := Select[
  Tuples[ops, Length@nums - 1],
  (target == ExprEval[nums, #]) &]

Usage example:

nums = {-1, 1, 2, 3};
ops = {Plus, Subtract, Times, Divide};
solutions = GetExpression[nums, ops, 3]

ExprEval[nums, #] & /@ solutions
SymbolicEval[nums, #] & /@ solutions

Outputs:

{{Plus, Times, Plus}, {Plus, Divide, Plus}, {Subtract, Plus, 
  Plus}, {Times, Plus, Times}, {Divide, Plus, Times}}

{3, 3, 3, 3, 3}

{"Plus"["Times"["Plus"[-1, 1], 2], 3], 
 "Plus"["Divide"["Plus"[-1, 1], 2], 3], 
 "Plus"["Plus"["Subtract"[-1, 1], 2], 3], 
 "Times"["Plus"["Times"[-1, 1], 2], 3], 
 "Times"["Plus"["Divide"[-1, 1], 2], 3]}

---

How it works

The ExprEval function takes in the numbers and operations, and applies them using (I think) RPN:

ExprEval[{1, 2, 3}, {Plus, Times}] == (1 + 2) * 3

It does this by continually folding pairs of numbers using the appropriate operation.

Now that I have a way to evaluate an expression tree, I just needed to generate them. Using Tuples, I'm able to generate all the different operators that I would intersperse between the numbers.

Once you get all possible operations, I used Select to pick out the the ones that evaluate to the target number.

---

Drawbacks

The solution above is really slow. Generating all the possible tuples is exponential in time. If there are k operations and n numbers, it's on the order of O(k^n).

For n = 10, it took 6 seconds to complete on Win 7 x64, Core i7 860, 12 GB RAM. The timings of the runs match the theoretical time complexity almost exactly:

Red line is the theoretical, blue is experimental. The x-axis is size of the nums input and the y-axis is the time in seconds to enumerate all solutions.

The above solution also solves the problem using a functional programming style. It looks pretty, but the thing also sucks up a butt ton of memory since it's storing the full results at nearly every step.

It doesn't even make use of parallelization, and I'm not entirely certain how you would even parallelize the solution I produced.

---

Some limitations

Mr. Wizard brought to my attention that this code only solves for only particular set of solutions. Given some input such as {a, b, c, d, e, ... } it only permutes the operators in between the numbers. It doesn't permute the ordering of the numbers. If it were to permute the numbers as well, the time complexity would rise up to O(k^n * n!) where k is the number of operators and n is the length of the input number array.

The following will produce the set of solutions for any permutation of the input numbers and operators:

(* generates a lists of the form 
{
 {number permutation, {{op order 1}, {op order 2}, ... }
 }, ...
}*)

GetAllExpressions[nums_, ops_, target_] :=
 ParallelMap[{#, GetExpression[#, ops, target]} &, 
  Tuples[nums, Length@nums]]