Substitutions

Generation of Substitutions that Simplify Complex Expressions

Written by Kenneth J. Hill (kenneth.hill@autodesk.com)
Creation date: Wednesday April 24, 2002
Revision date: Wednesday April 25, 2002 (Added line-torus intersection example, cleaned up a few items)

Discussion

It is often useful, especially when using Mathematica for software development, to apply substitutions to complex expressions to reduce their form.  For large expressions, this task can become tedious.  Substitutions[] was designed help with the process of finding a useful set of substitutions to simplify an expression.

The output of  Substitutions[inExp] (call it outExp) is a list of expressions representing substitutions of increasing complexity.  The (i^th) element of outExp is referred to as S[i] in subsequent expressions.   The final expression, outExp[[-1]], represents inExp in reduced form.  Thus the expression

is an identity.

Package Definition

\( (*\ Substitution\  - \ 
        Generation\ of\ Substitutions\ that\ Simplify\ Complex\ Expressions\ *) 
 (*\ 
      Written\ by\ Kenneth\ J . 
          Hill\ \((kenneth . hill@autodesk . 
            com)\)\ *) 
 (*\ 
      Creation\ \(date : Wednesday\ April\ 24\), 
      2002\ *) 
 (*\ 
      Revision\ \(date : Revision\ \(date : Wednesday\ April\ 25\)\), \ 
      2002\ *) \  (*\ \ \ \ \((Added\ line - torus\ intersection\ example, 
        cleaned\ up\ a\ few\ items)\)\ *) 

\(BeginPackage["Substitution`"];\)
\)
Substitutions::usage = "Substitutions[exp] yields a list of \
substitutions such that the ith element of the output list is referred to as \
S[i] in subsequent elements of the output list.  The last element of the \
output list represents exp in reduced form. If outExp = Substitution[exp], \
then exp can be recovered with outExp[[-1]] //. S[i_]:>outExp[[i]].  \
Substitutions[exp, ingoreList] specifies a list of patterns to ignore while \
making the substitutions.";
S::usage  =  "S[i] represents the ith substitution as generated \
in Substitutions[]";

Begin["`Private`"];

Return all pairings of list elements, without regard to order

allPairsHelper[l_] := 
If[l === {}, {}, 
        Module[{f  =  First[l], 
            tail = Drop[l, 1]},  {Thread[pair[f, tail]], 
            allPairsHelper[tail]}]]

allPairs[l_] := Flatten[allPairsHelper[l]]

Compute all sub expressions of a given expression.  For orderless functions with more than 2 parameters, the sub expressions are the function applied to each pair of parameters.

subexpressions[h_Symbol[x___]] := 
      If[MemberQ[Attributes[h], Orderless] && Length[{x}] > 2, 
        allPairs[List[x]] /. pair -> h, 
        List[x]]
subexpressions[h_?NumberQ] := h; 
    subexpressions[h_Symbol] := h;

inIgnoreList[exp_, ignore_] := 
      Or  @@  (MatchQ[exp, #] &  /@  ignore)

Remove all invalid candidates from the given list of candidates

filterCandidates[candidates_List, ignore_List] := 
      Select[candidates, 
(! (NumberQ[#]  ||  
                  Head[#] ===  Symbol  ||  Head[#] === SS  ||  
                  inIgnoreList[#, ignore])) &]

The main part of the algorithm: Attempt to use each sub expression at the current level as a substitution.  Only use the ones that occur more than once.

substHelper[exp_List, candidates_List, ignore_List, r_Integer] := 
        Module[
{
i, 
            n  =  Length[candidates], 
index  =  
              Length[exp], 
newExp
}, 

If[
            candidates  === {}  ||  r == 0,  
            Return[exp]]; 

 (* 
            Try each subexpression as a pattern *) 
For[i = 1,  
            i  <=  
              n,  ++i, 
newExp  =  (exp  //.  
                  candidates[[
                      i]] -> 
                    SS[index]); 

 (* 
              Test if newExp has this symbol more than once *) 
If[
              Length[Position[newExp, SS[index]]] ≥ 
                2, 
Return[
                  Append[newExp,  
                    candidates[[
                      i]]]];
];
]; 

 (* 
            No substitutions were found at this level,  
            so recurse *) 
substHelper[exp, 
            filterCandidates[
              Flatten[subexpressions[#] & /@ candidates, 1], ignore], 
            ignore, r - 1]
];

Substitutions[exp_, ignore_List: {}] := 
      Module[
{lastExp  =  {}, curExp, outExp, n, 
          i}, 

 
        curExp  =  
          substHelper[{exp}, 
            filterCandidates[Flatten[subexpressions[{exp}], 1], ignore], 
            ignore, 255]; 
While[
          curExp  ≠  lastExp, 
lastExp  =  
            curExp; 
curExp  =   substHelper[lastExp, 
              filterCandidates[Flatten[subexpressions[lastExp], 1], ignore], 
              ignore, 255];
]; 

 (* 
          Order curExp according to its dependency graph *) 
outExp  =  {}; 

n  =  
          Length[curExp]; 
While[
          Length[outExp] < n, 
 (* 
            Find an element of curExp that doesn' 
              t depend on SS *) 
For[i  =  1,  
              i  <=  n,  ++i, 
If[
                  Position[curExp[[i]], 
                      SS] === {}, 
AppendTo[outExp, 
                    curExp[[
                      i]]]; 
curExp[[i]]  =  
                    SS; 
curExp  =  (curExp  /.  
                        SS[i - 1] -> 
                          S[Length[
                              outExp]])
];
];
]; 

outExp
]

End[];
EndPackage[];

Example: Substitutions for Solving the Quatic Polynomial

Mathematica can symbolically solve the quartic (remove the ";" to see the full expression--its large)

quartic  =  
        Solve[a  +  b x  +  c (x^2) + d (x^3) + e (x^4) == 0, 
          x];

Substitutions creates a list sub such that when S[i] → sub[[i]], Last[sub] becomes the original expression (or possibly its equivalent).

sub  =  Substitutions[quartic]

Here are the individual substitutions:

GridBox[Table[{S[i], sub[[i]]}, {i, 1, Length[sub]}], 
        RowLines -> True, ColumnLines -> True, 
        RowSpacings -> 2] //  DisplayForm

We should be able to recover the quartic using the substitutions:

(Last[sub] //.  
          S[i_] :> sub[[i]])  ===  
      quartic

ClearAll[quartic, sub];

Example: Substitutions to Compute Line-Torus Intersections

Suppose one wishes to write a procedure to compute intersections between a line and a torus.  The answer involves computing the roots of a polynomial.  We'll assume the software package has a root finder, so the task is to find the the coefficients of the polynomial in a form that is easy to compute.

The implicit form or a torus with center at the origin, and axis of rotation along z:

torus  =  (((-(r^2) + (R^2) + (x^2) + (y^2) + (z^2)))^2) == 
        4 (R^2) ((x^2) + (y^2))

Substituting the parametric equation of a line (point-vector form)

intersectionEq[1]  =  
      torus  /.  {x -> px + t vx,  y -> py  +  t vy,  
          z -> pz  + t vz}

Putting all variables on the left:

intersectionEq[2]  =  
      Thread[intersectionEq[
            1] - intersectionEq[
              1][[2]], 
        Equal]

Expanding

intersectionEq[3]  =  ExpandAll[intersectionEq[2]]

Extracting coefficients

coefs[1]  = 
      CoefficientList[intersectionEq[
            3][[1]], t]

Horner simplifies computation and increases accuracy:

<< Algebra`Horner`;
coefs[2]  =  Horner  /@  coefs[1]

subs[1] = Substitutions[coefs[2]]

Here's a little procedure to count the number of multiplications and additions a C program would need to compute an expression.  It doesn't correspond perfectly to the operations a complier would actually preform (e.g. -a is counted as a product, (a^(m - n))can't be evaluated at all since the number of multiplications is unknown), but it gives a reasonable idea of the complexity of the expressions.

countMultAdd[expr_] := 
        Module[
{
opCount
}, 

opCount  =  
            Which[
Head[expr]  ===  
                Times, {Length[expr] - 1, 
                0}, 
Head[expr]  ===  Plus,  {0, 
                Length[expr] - 1}, 
Head[expr]  ===  
                Power,  
Which[
IntegerQ[
                    expr[[2]]]  &&  
                  expr[[2]] > 
                    0, {expr[[2]] - 1, 
                  0}, 
IntegerQ[
                    expr[[2]]]  &&  
                  expr[[2]] < 
                    0,  {expr[[2]], 
                  0}, 
True, {0, 
                  0}
], 
True,  {0, 
                0}
]; 

If[
            Length[expr] > 0, 
opCount +=  
              Plus  @@  (countMultAdd[#] &  /@  (List  @@  
                        expr))]; 

opCount
];

Before applying Horner:

Print["Multiplications: ", #1, " Additions: ", #2] &  @@  
      countMultAdd[coefs[1]]

After applying Horner:

Print["Multiplications: ", #1, " Additions: ", #2] &  @@  
      countMultAdd[coefs[2]]

After applying Substitutions:

Print["Multiplications: ", #1, " Additions: ", #2] &  @@  
      countMultAdd[subs[1]]

We should be able to recover coefs[2].  The next expression should be an identity:

(Last[subs[1]] //.  
          S[i_] :> subs[1][[i]]) === 
      coefs[2]

In Mathematica there is no real point to using subs[1] to evaluate the coefficients.  However, we do so for expository purposes:

Letting the major radius R = 2, and the minor radius r = 1, and letting the line be defined by the point {-3,-3,1} and vector {1,1,-1/10}:

R  =  2;  r  =  1;  p  =  {-3, -3, 1};  
    v  =  {1, 1, -1/10};

subs[2]  =  
        subs[1]  /.  {px -> p[[1]], py -> p[[2]], 
            pz -> p[[3]], vx -> v[[1]], 
            vy -> v[[2]], vz -> v[[3]]};

coefs[
        3]  =  subs[
            2][[-1]]  //.  
        S[i_] :> subs[2][[i]]

intersectionEq[4] = 
      coefs[3] . Table[(t^i), {i, 0, 4}] == 0

solnT  =  NSolve[intersectionEq[4], t]

points  =  Point  /@  (p  +  t v  /.  solnT)

Showing the solutions graphically (execute the next statement to see the image):

<< Graphics`Shapes`;
Show[
      Graphics3D[{Torus[R, r], {RGBColor[0, 1, 0], 
            Line[{p, p + 6 v}]}, {RGBColor[1, 0, 0], PointSize[ .02], 
            points}}], Boxed -> False]

ClearAll[torus, intersectionEq, coefs, subs, p, v, solnT, 
        points];