(************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 28341, 746]*) (*NotebookOutlinePosition[ 29078, 771]*) (* CellTagsIndexPosition[ 29034, 767]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Substitutions", "Title"], Cell["\<\ Generation of Substitutions that Simplify Complex Expressions\ \>", "Subtitle"], Cell["\<\ 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)\ \>", "Subsubtitle"], Cell[CellGroupData[{ Cell["Discussion", "Section"], Cell[TextData[{ "It is often useful, especially when using ", StyleBox["Mathematica", FontSlant->"Italic"], " 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." }], "Text"], Cell[TextData[{ "The output of Substitutions[inExp] (call it outExp) is a list of \ expressions representing substitutions of increasing complexity. The ", Cell[BoxData[ \(TraditionalForm\`i\^th\)]], " element of outExp is referred to as ", Cell[BoxData[ \(TraditionalForm\`S[i]\)]], " in subsequent expressions. The final expression, outExp[[-1]], \ represents inExp in reduced form. Thus the expression" }], "Text"], Cell[BoxData[ \(\((outExp[\([\(-1\)]\)] //. S[i] \[RightArrow] outExp[\([i]\)])\) === inExp\)], "Output"], Cell["is an identity.", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Package Definition", "Section", InitializationCell->True], Cell[BoxData[{ \( (*\ Substitution\ - \ Generation\ of\ Substitutions\ that\ Simplify\ Complex\ Expressions\ *) \ \[IndentingNewLine] (*\ Written\ by\ Kenneth\ J . Hill\ \((kenneth . hill@autodesk . com)\)\ *) \[IndentingNewLine] (*\ Creation\ \(date : Wednesday\ April\ 24\), 2002\ *) \[IndentingNewLine] (*\ Revision\ \(date : Revision\ \(date : Wednesday\ April\ 25\)\), \ 2002\ *) \ (*\ \ \ \ \((Added\ line - torus\ intersection\ example, cleaned\ up\ a\ few\ items)\)\ *) \n\ \[IndentingNewLine]\(BeginPackage["\"];\)\[IndentingNewLine]\)\ , "\n", \(\(Substitutions::usage = "\outExp[[i]]. \ Substitutions[exp, ingoreList] specifies a list of patterns to ignore while \ making the substitutions.\>";\)\), "\n", \(\(S::usage\ = \ "\";\)\[IndentingNewLine]\), "\[IndentingNewLine]", \(\(Begin["\<`Private`\>"];\)\)}], "Input", InitializationCell->True], Cell["\<\ Return all pairings of list elements, without regard to order\ \>", "Text"], Cell[BoxData[ \(allPairsHelper[l_] := \[IndentingNewLine]If[l === {}, {}, Module[{f\ = \ First[l], tail = Drop[l, 1]}, \ {Thread[pair[f, tail]], allPairsHelper[tail]}]]\)], "Input", InitializationCell->True], Cell[BoxData[ \(allPairs[l_] := Flatten[allPairsHelper[l]]\)], "Input", InitializationCell->True], Cell["\<\ 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.\ \>", "Text"], Cell[BoxData[{ \(subexpressions[h_Symbol[x___]] := If[MemberQ[Attributes[h], Orderless] && Length[{x}] > 2, allPairs[List[x]] /. pair \[Rule] h, List[x]]\), "\[IndentingNewLine]", \(subexpressions[h_?NumberQ] := h; subexpressions[h_Symbol] := h;\)}], "Input", InitializationCell->True], Cell[BoxData[ \(inIgnoreList[exp_, ignore_] := Or\ @@ \ \((\(MatchQ[exp, #] &\)\ /@ \ ignore)\)\)], "Input", InitializationCell->True], Cell["\<\ Remove all invalid candidates from the given list of candidates\ \>", "Text"], Cell[BoxData[ \(filterCandidates[candidates_List, ignore_List] := Select[candidates, \[IndentingNewLine]\((\(! \((NumberQ[#]\ || \ Head[#] === \ Symbol\ || \ Head[#] === SS\ || \ inIgnoreList[#, ignore])\)\))\) &]\)], "Input", InitializationCell->True], Cell["\<\ 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.\ \>", "Text"], Cell[BoxData[ \(\(substHelper[exp_List, candidates_List, ignore_List, r_Integer] := Module[\[IndentingNewLine]{\[IndentingNewLine]i, n\ = \ Length[candidates], \[IndentingNewLine]index\ = \ Length[exp], \[IndentingNewLine]newExp\[IndentingNewLine]}, \ \[IndentingNewLine]\[IndentingNewLine]If[ candidates\ === {}\ || \ r \[Equal] 0, \ Return[exp]]; \[IndentingNewLine]\[IndentingNewLine] (*\ Try\ each\ subexpression\ as\ a\ pattern\ \ *) \[IndentingNewLine]For[i = 1, \ i\ <= \ n, \ \(++i\), \[IndentingNewLine]newExp\ = \ \((exp\ //. \ candidates\[LeftDoubleBracket] i\[RightDoubleBracket] \[Rule] SS[index])\); \[IndentingNewLine]\[IndentingNewLine] (*\ Test\ if\ newExp\ has\ this\ symbol\ more\ than\ once\ \ *) \[IndentingNewLine]If[ Length[Position[newExp, SS[index]]] \[GreaterEqual] 2, \[IndentingNewLine]\(Return[ Append[newExp, \ candidates\[LeftDoubleBracket] i\[RightDoubleBracket]]];\)\[IndentingNewLine]];\ \[IndentingNewLine]]; \[IndentingNewLine]\[IndentingNewLine] (*\ No\ substitutions\ were\ found\ at\ this\ level, \ so\ recurse\ *) \[IndentingNewLine]substHelper[exp, filterCandidates[ Flatten[\(subexpressions[#] &\) /@ candidates, 1], ignore], ignore, r - 1]\[IndentingNewLine]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(Substitutions[exp_, ignore_List: {}] := Module[\[IndentingNewLine]{lastExp\ = \ {}, curExp, outExp, n, i}, \[IndentingNewLine]\[IndentingNewLine]\ curExp\ = \ substHelper[{exp}, filterCandidates[Flatten[subexpressions[{exp}], 1], ignore], ignore, 255]; \[IndentingNewLine]While[ curExp\ \[NotEqual] \ lastExp, \[IndentingNewLine]lastExp\ = \ curExp; \[IndentingNewLine]curExp\ = \ \ substHelper[lastExp, filterCandidates[Flatten[subexpressions[lastExp], 1], ignore], ignore, 255];\[IndentingNewLine]]; \[IndentingNewLine]\ \[IndentingNewLine] (*\ Order\ curExp\ according\ to\ its\ dependency\ graph\ \ *) \[IndentingNewLine]outExp\ = \ {}; \ \[IndentingNewLine]\[IndentingNewLine]n\ = \ Length[curExp]; \[IndentingNewLine]While[ Length[outExp] < n, \[IndentingNewLine] (*\ Find\ an\ element\ of\ curExp\ that\ doesn' t\ depend\ on\ SS\ *) \[IndentingNewLine]\(For[i\ = \ 1, \ i\ <= \ n, \ \(++i\), \[IndentingNewLine]\(If[ Position[curExp\[LeftDoubleBracket]i\[RightDoubleBracket], SS] === {}, \[IndentingNewLine]AppendTo[outExp, curExp\[LeftDoubleBracket] i\[RightDoubleBracket]]; \[IndentingNewLine]curExp\ \[LeftDoubleBracket]i\[RightDoubleBracket]\ = \ SS; \[IndentingNewLine]curExp\ = \ \((curExp\ /. \ SS[i - 1] \[Rule] S[Length[ outExp]])\)\[IndentingNewLine]];\)\ \[IndentingNewLine]];\)\[IndentingNewLine]]; \[IndentingNewLine]\ \[IndentingNewLine]outExp\[IndentingNewLine]]\)], "Input", InitializationCell->True], Cell[BoxData[{ \(\(End[];\)\), "\[IndentingNewLine]", \(\(EndPackage[];\)\)}], "Input", InitializationCell->True] }, Closed]], Cell[CellGroupData[{ Cell["Example: Substitutions for Solving the Quatic Polynomial", "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " can symbolically solve the quartic (remove the \";\" to see the full \ expression--its large)" }], "Text"], Cell[BoxData[ \(\(quartic\ = \ Solve[a\ + \ b\ x\ + \ c\ x\^2 + d\ x\^3 + e\ x\^4 \[Equal] 0, x];\)\)], "Input"], Cell[TextData[{ "Substitutions creates a list ", StyleBox["sub", FontSlant->"Italic"], " such that when S[i] \[RightArrow] sub[[i]], Last[sub] becomes the \ original expression (or possibly its equivalent)." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(sub\ = \ Substitutions[quartic]\)], "Input"], Cell[BoxData[ \({1\/e, d\^2, 1\/e\^2, c\ S[1], S[1]\/3, b\ d, a\ e, \(-\(S[1]\/3\)\), \(-\(1\/4\)\)\ d\ S[1], S[2]\ S[3], c\^2 - 3\ S[6] + 12\ S[7], 2\ c\^3 - 9\ c\ S[6], 27\ b\^2\ e - 72\ c\ S[7], 27\ a\ S[2] + S[12], S[13] + S[14], \(-\(\(4\ S[4]\)\/3\)\) + S[10]\/2, 2\^\(1/3\)\ S[11], S[15] + \@\(\(-4\)\ S[11]\^3 + S[15]\^2\), S[18]\^\(1/3\)\/2\^\(1/3\), S[17]\/S[18]\^\(1/3\), S[8]\ S[19] + S[8]\ S[20], S[16] + S[21], \(-\(\(2\ S[4]\)\/3\)\) + S[10]\/4 + S[5]\ S[19] + S[5]\ S[20], \(\(-\(d\^3\/e\^3\)\) - 8\ b\ S[1] + 4\ c\ d\ \ S[3]\)\/\@S[23], \@\(S[22] + S[24]\/4\), \@S[23], S[9] + S[26]\/2, \@\(S[22] - S[24]\/4\), S[9] - S[26]\/2, {{x \[Rule] \(-\(S[28]\/2\)\) + S[29]}, {x \[Rule] S[28]\/2 + S[29]}, {x \[Rule] \(-\(S[25]\/2\)\) + S[27]}, {x \[Rule] S[25]\/2 + S[27]}}}\)], "Output"] }, Open ]], Cell["Here are the individual substitutions:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(GridBox[Table[{S[i], sub[\([i]\)]}, {i, 1, Length[sub]}], RowLines \[Rule] True, ColumnLines \[Rule] True, RowSpacings \[Rule] 2] // \ DisplayForm\)], "Input"], Cell[BoxData[ TagBox[GridBox[{ {\(S[1]\), \(1\/e\)}, {\(S[2]\), \(d\^2\)}, {\(S[3]\), \(1\/e\^2\)}, {\(S[4]\), \(c\ S[1]\)}, {\(S[5]\), \(S[1]\/3\)}, {\(S[6]\), \(b\ d\)}, {\(S[7]\), \(a\ e\)}, {\(S[8]\), \(-\(S[1]\/3\)\)}, {\(S[9]\), \(\(-\(1\/4\)\)\ d\ S[1]\)}, {\(S[10]\), \(S[2]\ S[3]\)}, {\(S[11]\), \(c\^2 - 3\ S[6] + 12\ S[7]\)}, {\(S[12]\), \(2\ c\^3 - 9\ c\ S[6]\)}, {\(S[13]\), \(27\ b\^2\ e - 72\ c\ S[7]\)}, {\(S[14]\), \(27\ a\ S[2] + S[12]\)}, {\(S[15]\), \(S[13] + S[14]\)}, {\(S[16]\), \(\(-\(\(4\ S[4]\)\/3\)\) + S[10]\/2\)}, {\(S[17]\), \(2\^\(1/3\)\ S[11]\)}, {\(S[18]\), \(S[15] + \@\(\(-4\)\ S[11]\^3 + S[15]\^2\)\)}, {\(S[19]\), \(S[18]\^\(1/3\)\/2\^\(1/3\)\)}, {\(S[20]\), \(S[17]\/S[18]\^\(1/3\)\)}, {\(S[21]\), \(S[8]\ S[19] + S[8]\ S[20]\)}, {\(S[22]\), \(S[16] + S[21]\)}, {\(S[23]\), \(\(-\(\(2\ S[4]\)\/3\)\) + S[10]\/4 + S[5]\ S[19] + S[5]\ S[20]\)}, {\(S[ 24]\), \(\(\(-\(d\^3\/e\^3\)\) - 8\ b\ S[1] + 4\ c\ d\ S[3]\)\/\@S[23]\)}, {\(S[25]\), \(\@\(S[22] + S[24]\/4\)\)}, {\(S[26]\), \(\@S[23]\)}, {\(S[27]\), \(S[9] + S[26]\/2\)}, {\(S[28]\), \(\@\(S[22] - S[24]\/4\)\)}, {\(S[29]\), \(S[9] - S[26]\/2\)}, {\(S[ 30]\), \({{x \[Rule] \(-\(S[28]\/2\)\) + S[29]}, {x \[Rule] S[28]\/2 + S[29]}, {x \[Rule] \(-\(S[25]\/2\)\) + S[27]}, {x \[Rule] S[25]\/2 + S[27]}}\)} }, RowSpacings->2, RowLines->True, ColumnLines->True], DisplayForm]], "Output"] }, Open ]], Cell["\<\ We should be able to recover the quartic using the substitutions:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\((Last[sub] //. \ S[i_] :> sub\[LeftDoubleBracket]i\[RightDoubleBracket])\)\ === \ quartic\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Open ]], Cell[BoxData[ \(\(ClearAll[quartic, sub];\)\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Example: Substitutions to Compute Line-Torus Intersections", "Section"], Cell["\<\ 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.\ \>", "Text"], Cell["\<\ The implicit form or a torus with center at the origin, and axis of rotation \ along z:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(torus\ = \ \((\(-r\^2\) + R\^2 + x\^2 + y\^2 + z\^2)\)\^2 == 4\ R\^2\ \((x\^2 + y\^2)\)\)], "Input"], Cell[BoxData[ \(\((\(-r\^2\) + R\^2 + x\^2 + y\^2 + z\^2)\)\^2 == 4\ R\^2\ \((x\^2 + y\^2)\)\)], "Output"] }, Open ]], Cell["\<\ Substituting the parametric equation of a line (point-vector form)\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(intersectionEq[1]\ = \ torus\ /. \ {x \[Rule] px + t\ vx, \ y \[Rule] py\ + \ t\ vy, \ z \[Rule] pz\ + t\ vz}\)], "Input"], Cell[BoxData[ \(\((\(-r\^2\) + R\^2 + \((px + t\ vx)\)\^2 + \((py + t\ vy)\)\^2 + \((pz \ + t\ vz)\)\^2)\)\^2 == 4\ R\^2\ \((\((px + t\ vx)\)\^2 + \((py + t\ vy)\)\^2)\)\)], "Output"] }, Open ]], Cell["Putting all variables on the left:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(intersectionEq[2]\ = \ Thread[intersectionEq[ 1] - \(intersectionEq[ 1]\)\[LeftDoubleBracket]2\[RightDoubleBracket], Equal]\)], "Input"], Cell[BoxData[ \(\(-4\)\ R\^2\ \((\((px + t\ vx)\)\^2 + \((py + t\ vy)\)\^2)\) + \((\(-r\ \^2\) + R\^2 + \((px + t\ vx)\)\^2 + \((py + t\ vy)\)\^2 + \((pz + t\ \ vz)\)\^2)\)\^2 == 0\)], "Output"] }, Open ]], Cell["Expanding", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(intersectionEq[3]\ = \ ExpandAll[intersectionEq[2]]\)], "Input"], Cell[BoxData[ \(px\^4 + 2\ px\^2\ py\^2 + py\^4 + 2\ px\^2\ pz\^2 + 2\ py\^2\ pz\^2 + pz\^4 - 2\ px\^2\ r\^2 - 2\ py\^2\ r\^2 - 2\ pz\^2\ r\^2 + r\^4 - 2\ px\^2\ R\^2 - 2\ py\^2\ R\^2 + 2\ pz\^2\ R\^2 - 2\ r\^2\ R\^2 + R\^4 + 4\ px\^3\ t\ vx + 4\ px\ py\^2\ t\ vx + 4\ px\ pz\^2\ t\ vx - 4\ px\ r\^2\ t\ vx - 4\ px\ R\^2\ t\ vx + 6\ px\^2\ t\^2\ vx\^2 + 2\ py\^2\ t\^2\ vx\^2 + 2\ pz\^2\ t\^2\ vx\^2 - 2\ r\^2\ t\^2\ vx\^2 - 2\ R\^2\ t\^2\ vx\^2 + 4\ px\ t\^3\ vx\^3 + t\^4\ vx\^4 + 4\ px\^2\ py\ t\ vy + 4\ py\^3\ t\ vy + 4\ py\ pz\^2\ t\ vy - 4\ py\ r\^2\ t\ vy - 4\ py\ R\^2\ t\ vy + 8\ px\ py\ t\^2\ vx\ vy + 4\ py\ t\^3\ vx\^2\ vy + 2\ px\^2\ t\^2\ vy\^2 + 6\ py\^2\ t\^2\ vy\^2 + 2\ pz\^2\ t\^2\ vy\^2 - 2\ r\^2\ t\^2\ vy\^2 - 2\ R\^2\ t\^2\ vy\^2 + 4\ px\ t\^3\ vx\ vy\^2 + 2\ t\^4\ vx\^2\ vy\^2 + 4\ py\ t\^3\ vy\^3 + t\^4\ vy\^4 + 4\ px\^2\ pz\ t\ vz + 4\ py\^2\ pz\ t\ vz + 4\ pz\^3\ t\ vz - 4\ pz\ r\^2\ t\ vz + 4\ pz\ R\^2\ t\ vz + 8\ px\ pz\ t\^2\ vx\ vz + 4\ pz\ t\^3\ vx\^2\ vz + 8\ py\ pz\ t\^2\ vy\ vz + 4\ pz\ t\^3\ vy\^2\ vz + 2\ px\^2\ t\^2\ vz\^2 + 2\ py\^2\ t\^2\ vz\^2 + 6\ pz\^2\ t\^2\ vz\^2 - 2\ r\^2\ t\^2\ vz\^2 + 2\ R\^2\ t\^2\ vz\^2 + 4\ px\ t\^3\ vx\ vz\^2 + 2\ t\^4\ vx\^2\ vz\^2 + 4\ py\ t\^3\ vy\ vz\^2 + 2\ t\^4\ vy\^2\ vz\^2 + 4\ pz\ t\^3\ vz\^3 + t\^4\ vz\^4 == 0\)], "Output"] }, Open ]], Cell["Extracting coefficients", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(coefs[1]\ = CoefficientList[\(intersectionEq[ 3]\)\[LeftDoubleBracket]1\[RightDoubleBracket], t]\)], "Input"], Cell[BoxData[ \({px\^4 + 2\ px\^2\ py\^2 + py\^4 + 2\ px\^2\ pz\^2 + 2\ py\^2\ pz\^2 + pz\^4 - 2\ px\^2\ r\^2 - 2\ py\^2\ r\^2 - 2\ pz\^2\ r\^2 + r\^4 - 2\ px\^2\ R\^2 - 2\ py\^2\ R\^2 + 2\ pz\^2\ R\^2 - 2\ r\^2\ R\^2 + R\^4, 4\ px\^3\ vx + 4\ px\ py\^2\ vx + 4\ px\ pz\^2\ vx - 4\ px\ r\^2\ vx - 4\ px\ R\^2\ vx + 4\ px\^2\ py\ vy + 4\ py\^3\ vy + 4\ py\ pz\^2\ vy - 4\ py\ r\^2\ vy - 4\ py\ R\^2\ vy + 4\ px\^2\ pz\ vz + 4\ py\^2\ pz\ vz + 4\ pz\^3\ vz - 4\ pz\ r\^2\ vz + 4\ pz\ R\^2\ vz, 6\ px\^2\ vx\^2 + 2\ py\^2\ vx\^2 + 2\ pz\^2\ vx\^2 - 2\ r\^2\ vx\^2 - 2\ R\^2\ vx\^2 + 8\ px\ py\ vx\ vy + 2\ px\^2\ vy\^2 + 6\ py\^2\ vy\^2 + 2\ pz\^2\ vy\^2 - 2\ r\^2\ vy\^2 - 2\ R\^2\ vy\^2 + 8\ px\ pz\ vx\ vz + 8\ py\ pz\ vy\ vz + 2\ px\^2\ vz\^2 + 2\ py\^2\ vz\^2 + 6\ pz\^2\ vz\^2 - 2\ r\^2\ vz\^2 + 2\ R\^2\ vz\^2, 4\ px\ vx\^3 + 4\ py\ vx\^2\ vy + 4\ px\ vx\ vy\^2 + 4\ py\ vy\^3 + 4\ pz\ vx\^2\ vz + 4\ pz\ vy\^2\ vz + 4\ px\ vx\ vz\^2 + 4\ py\ vy\ vz\^2 + 4\ pz\ vz\^3, vx\^4 + 2\ vx\^2\ vy\^2 + vy\^4 + 2\ vx\^2\ vz\^2 + 2\ vy\^2\ vz\^2 + vz\^4}\)], "Output"] }, Open ]], Cell["Horner simplifies computation and increases accuracy:", "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(\(<< Algebra`Horner`;\)\), "\[IndentingNewLine]", \(coefs[2]\ = \ Horner\ /@ \ coefs[1]\)}], "Input"], Cell[BoxData[ \({R\^4 + py\^2\ \((py\^2 + 2\ pz\^2 - 2\ r\^2 - 2\ R\^2)\) + px\^2\ \((px\^2 + 2\ py\^2 + 2\ pz\^2 - 2\ r\^2 - 2\ R\^2)\) + r\^2\ \((r\^2 - 2\ R\^2)\) + pz\^2\ \((pz\^2 - 2\ r\^2 + 2\ R\^2)\), pz\ \((4\ pz\^2\ vz - 4\ r\^2\ vz + 4\ R\^2\ vz)\) + py\ \((4\ pz\^2\ vy - 4\ r\^2\ vy - 4\ R\^2\ vy + py\ \((4\ py\ vy + 4\ pz\ vz)\))\) + px\ \((4\ py\^2\ vx + 4\ pz\^2\ vx - 4\ r\^2\ vx - 4\ R\^2\ vx + px\ \((4\ px\ vx + 4\ py\ vy + 4\ pz\ vz)\))\), r\^2\ \((\(-2\)\ vx\^2 - 2\ vy\^2 - 2\ vz\^2)\) + R\^2\ \((\(-2\)\ vx\^2 - 2\ vy\^2 + 2\ vz\^2)\) + pz\^2\ \((2\ vx\^2 + 2\ vy\^2 + 6\ vz\^2)\) + px\ \((8\ py\ vx\ vy + 8\ pz\ vx\ vz + px\ \((6\ vx\^2 + 2\ vy\^2 + 2\ vz\^2)\))\) + py\ \((8\ pz\ vy\ vz + py\ \((2\ vx\^2 + 6\ vy\^2 + 2\ vz\^2)\))\), px\ vx\ \((4\ vx\^2 + 4\ vy\^2 + 4\ vz\^2)\) + pz\ \((4\ vx\^2\ vz + 4\ vy\^2\ vz + 4\ vz\^3)\) + py\ \((4\ vx\^2\ vy + vy\ \((4\ vy\^2 + 4\ vz\^2)\))\), vz\^4 + vy\^2\ \((vy\^2 + 2\ vz\^2)\) + vx\^2\ \((vx\^2 + 2\ vy\^2 + 2\ vz\^2)\)}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(subs[1] = Substitutions[coefs[2]]\)], "Input"], Cell[BoxData[ \(General::"spell1" \(\(:\)\(\ \)\) "Possible spelling error: new symbol name \"\!\(subs\)\" is similar to \ existing symbol \"\!\(sub\)\"."\)], "Message"], Cell[BoxData[ \({py\^2, r\^2, pz\^2, px\^2, R\^2, px\ vx, vy\^2, vx\^2, \(-2\)\ S[5], \(-2\)\ S[2], \(-2\)\ S[7] - 2\ S[8], vz\^2, 4\ S[7], 4\ S[8], 2\ S[7], 8\ pz, 2\ S[8], vz\ S[16], 4\ vy, 4\ vz, 4\ vx, \(-4\)\ vy, \(-4\)\ S[2], 2\ S[3] + S[10], S[9] + S[24], 2\ S[12], 4\ S[12] + S[13], S[15] + S[26], py\ S[19] + pz\ S[20], {R\^4 + S[2]\ \((S[2] + S[9])\) + S[3]\ \((S[3] + 2\ S[5] + S[10])\) + S[1]\ \((S[1] + S[25])\) + S[4]\ \((2\ S[1] + S[4] + S[25])\), pz\ \((S[3]\ S[20] + S[5]\ S[20] + vz\ S[23])\) + py\ \((S[3]\ S[19] + S[2]\ S[22] + S[5]\ S[22] + py\ S[29])\) + px\ \((\(-4\)\ vx\ S[5] + S[1]\ S[21] + S[3]\ S[21] + vx\ S[23] + px\ \((4\ S[6] + S[29])\))\), S[2]\ \((S[11] - 2\ S[12])\) + S[3]\ \((6\ S[12] + S[15] + S[17])\) + S[5]\ \((S[11] + S[26])\) + py\ \((vy\ S[18] + py\ \((6\ S[7] + S[17] + S[26])\))\) + px\ \((8\ py\ vx\ vy + vx\ S[18] + px\ \((6\ S[8] + S[28])\))\), pz\ \((4\ vz\^3 + vz\ S[13] + vz\ S[14])\) + S[6]\ \((S[14] + S[27])\) + py\ \((vy\ S[14] + vy\ S[27])\), vz\^4 + S[7]\ \((S[7] + S[26])\) + S[8]\ \((S[8] + S[28])\)}}\)], "Output"] }, Open ]], Cell[TextData[{ "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, ", Cell[BoxData[ \(TraditionalForm\`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." }], "Text"], Cell[BoxData[ \(\(countMultAdd[expr_] := Module[\[IndentingNewLine]{\[IndentingNewLine]opCount\ \[IndentingNewLine]}, \[IndentingNewLine]\[IndentingNewLine]opCount\ = \ Which[\[IndentingNewLine]Head[expr]\ === \ Times, {Length[expr] - 1, 0}, \[IndentingNewLine]Head[expr]\ === \ Plus, \ {0, Length[expr] - 1}, \[IndentingNewLine]Head[expr]\ === \ Power, \ \[IndentingNewLine]Which[\[IndentingNewLine]IntegerQ[ expr\[LeftDoubleBracket]2\[RightDoubleBracket]]\ && \ expr\[LeftDoubleBracket]2\[RightDoubleBracket] > 0, {expr\[LeftDoubleBracket]2\[RightDoubleBracket] - 1, 0}, \[IndentingNewLine]IntegerQ[ expr\[LeftDoubleBracket]2\[RightDoubleBracket]]\ && \ expr\[LeftDoubleBracket]2\[RightDoubleBracket] < 0, \ {expr\[LeftDoubleBracket]2\[RightDoubleBracket], 0}, \[IndentingNewLine]True, {0, 0}\[IndentingNewLine]], \[IndentingNewLine]True, \ {0, 0}\[IndentingNewLine]]; \ \[IndentingNewLine]\[IndentingNewLine]If[ Length[expr] > 0, \[IndentingNewLine]opCount += \ Plus\ @@ \ \((\(countMultAdd[#] &\)\ /@ \ \((List\ @@ \ expr)\))\)]; \ \[IndentingNewLine]\[IndentingNewLine]opCount\[IndentingNewLine]];\)\)], \ "Input"], Cell["Before applying Horner:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(Print["\", #1, "\< Additions: \>", #2] &\)\ @@ \ countMultAdd[coefs[1]]\)], "Input"], Cell[BoxData[ InterpretationBox[\("Multiplications: "\[InvisibleSpace]244\ \[InvisibleSpace]" Additions: "\[InvisibleSpace]58\), SequenceForm[ "Multiplications: ", 244, " Additions: ", 58], Editable->False]], "Print"] }, Open ]], Cell["After applying Horner:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(Print["\", #1, "\< Additions: \>", #2] &\)\ @@ \ countMultAdd[coefs[2]]\)], "Input"], Cell[BoxData[ InterpretationBox[\("Multiplications: "\[InvisibleSpace]171\ \[InvisibleSpace]" Additions: "\[InvisibleSpace]58\), SequenceForm[ "Multiplications: ", 171, " Additions: ", 58], Editable->False]], "Print"] }, Open ]], Cell["After applying Substitutions:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(Print["\", #1, "\< Additions: \>", #2] &\)\ @@ \ countMultAdd[subs[1]]\)], "Input"], Cell[BoxData[ InterpretationBox[\("Multiplications: "\[InvisibleSpace]86\ \[InvisibleSpace]" Additions: "\[InvisibleSpace]52\), SequenceForm[ "Multiplications: ", 86, " Additions: ", 52], Editable->False]], "Print"] }, Open ]], Cell["\<\ We should be able to recover coefs[2]. The next expression should be an \ identity:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\((Last[subs[1]] //. \ S[i_] :> \(subs[1]\)\[LeftDoubleBracket]i\[RightDoubleBracket])\) === coefs[2]\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Open ]], Cell[TextData[{ "In ", StyleBox["Mathematica", FontSlant->"Italic"], " there is no real point to using subs[1] to evaluate the coefficients. \ However, we do so for expository purposes:" }], "Text"], Cell["\<\ 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}:\ \>", "Text"], Cell[BoxData[ \(R\ = \ 2; \ r\ = \ 1; \ p\ = \ {\(-3\), \(-3\), 1}; \ v\ = \ {1, 1, \(-1\)/10};\)], "Input"], Cell[BoxData[ \(\(subs[2]\ = \ subs[1]\ /. \ {px \[Rule] p[\([1]\)], py \[Rule] p[\([2]\)], pz \[Rule] p[\([3]\)], vx \[Rule] v[\([1]\)], vy \[Rule] v[\([2]\)], vz \[Rule] v[\([3]\)]};\)\)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(coefs[ 3]\ = \ \(subs[ 2]\)\[LeftDoubleBracket]\(-1\)\[RightDoubleBracket]\ //. \ S[i_] :> \(subs[2]\)[\([i]\)]\)], "Input"], Cell[BoxData[ \({196, \(-\(1724\/5\)\), 5132\/25, \(-\(12261\/250\)\), 40401\/10000}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(intersectionEq[4] = coefs[3] . Table[t\^i, {i, 0, 4}] \[Equal] 0\)], "Input"], Cell[BoxData[ \(196 - \(1724\ t\)\/5 + \(5132\ t\^2\)\/25 - \(12261\ t\^3\)\/250 + \ \(40401\ t\^4\)\/10000 == 0\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(solnT\ = \ NSolve[intersectionEq[4], t]\)], "Input"], Cell[BoxData[ \({{t \[Rule] 1.2442018976458646`}, {t \[Rule] 2.0110944956415104`}, {t \[Rule] 3.8562938769624013`}, {t \[Rule] 5.02771321233729`}}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(points\ = \ Point\ /@ \ \((p\ + \ t\ v\ /. \ solnT)\)\)], "Input"], Cell[BoxData[ \({Point[{\(-1.7557981023541354`\), \(-1.7557981023541354`\), 0.8755798102354135`}], Point[{\(-0.9889055043584896`\), \(-0.9889055043584896`\), 0.798890550435849`}], Point[{0.8562938769624013`, 0.8562938769624013`, 0.6143706123037599`}], Point[{2.0277132123372903`, 2.0277132123372903`, 0.49722867876627097`}]}\)], "Output"] }, Open ]], Cell["\<\ Showing the solutions graphically (execute the next statement to see the \ image):\ \>", "Text"], Cell[BoxData[{ \(\(<< Graphics`Shapes`;\)\), "\[IndentingNewLine]", \(Show[ Graphics3D[{Torus[R, r], {RGBColor[0, 1, 0], Line[{p, p + 6\ v}]}, {RGBColor[1, 0, 0], PointSize[ .02], points}}], Boxed \[Rule] False]\)}], "Input"], Cell[BoxData[ \(\(ClearAll[torus, intersectionEq, coefs, subs, p, v, solnT, points];\)\)], "Input"] }, Closed]] }, Open ]] }, FrontEndVersion->"4.1 for Microsoft Windows", ScreenRectangle->{{0, 1280}, {0, 951}}, AutoGeneratedPackage->Automatic, WindowSize->{930, 740}, WindowMargins->{{48, Automatic}, {Automatic, 36}}, ShowSelection->True, StyleDefinitions -> "PastelColor.nb" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1727, 52, 30, 0, 98, "Title"], Cell[1760, 54, 89, 2, 41, "Subtitle"], Cell[1852, 58, 231, 5, 58, "Subsubtitle"], Cell[CellGroupData[{ Cell[2108, 67, 29, 0, 51, "Section"], Cell[2140, 69, 406, 8, 50, "Text"], Cell[2549, 79, 444, 10, 50, "Text"], Cell[2996, 91, 116, 2, 37, "Output"], Cell[3115, 95, 31, 0, 29, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[3183, 100, 65, 1, 51, "Section", InitializationCell->True], Cell[3251, 103, 1372, 24, 358, "Input", InitializationCell->True], Cell[4626, 129, 85, 2, 29, "Text"], Cell[4714, 133, 249, 5, 58, "Input", InitializationCell->True], Cell[4966, 140, 103, 2, 38, "Input", InitializationCell->True], Cell[5072, 144, 205, 4, 50, "Text"], Cell[5280, 150, 326, 7, 98, "Input", InitializationCell->True], Cell[5609, 159, 149, 3, 38, "Input", InitializationCell->True], Cell[5761, 164, 87, 2, 29, "Text"], Cell[5851, 168, 306, 5, 58, "Input", InitializationCell->True], Cell[6160, 175, 180, 4, 29, "Text"], Cell[6343, 181, 1595, 28, 458, "Input", InitializationCell->True], Cell[7941, 211, 1836, 32, 538, "Input", InitializationCell->True], Cell[9780, 245, 124, 3, 58, "Input", InitializationCell->True] }, Closed]], Cell[CellGroupData[{ Cell[9941, 253, 75, 0, 35, "Section"], Cell[10019, 255, 178, 5, 29, "Text"], Cell[10200, 262, 139, 3, 39, "Input"], Cell[10342, 267, 229, 6, 29, "Text"], Cell[CellGroupData[{ Cell[10596, 277, 65, 1, 38, "Input"], Cell[10664, 280, 902, 14, 212, "Output"] }, Open ]], Cell[11581, 297, 54, 0, 29, "Text"], Cell[CellGroupData[{ Cell[11660, 301, 197, 3, 58, "Input"], Cell[11860, 306, 1820, 42, 943, "Output"] }, Open ]], Cell[13695, 351, 89, 2, 29, "Text"], Cell[CellGroupData[{ Cell[13809, 357, 143, 3, 38, "Input"], Cell[13955, 362, 38, 1, 37, "Output"] }, Open ]], Cell[14008, 366, 60, 1, 38, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[14105, 372, 77, 0, 35, "Section"], Cell[14185, 374, 331, 5, 50, "Text"], Cell[14519, 381, 111, 3, 29, "Text"], Cell[CellGroupData[{ Cell[14655, 388, 130, 2, 40, "Input"], Cell[14788, 392, 117, 2, 38, "Output"] }, Open ]], Cell[14920, 397, 90, 2, 29, "Text"], Cell[CellGroupData[{ Cell[15035, 403, 164, 3, 38, "Input"], Cell[15202, 408, 194, 3, 38, "Output"] }, Open ]], Cell[15411, 414, 50, 0, 29, "Text"], Cell[CellGroupData[{ Cell[15486, 418, 199, 5, 38, "Input"], Cell[15688, 425, 198, 3, 38, "Output"] }, Open ]], Cell[15901, 431, 25, 0, 29, "Text"], Cell[CellGroupData[{ Cell[15951, 435, 85, 1, 38, "Input"], Cell[16039, 438, 1506, 22, 132, "Output"] }, Open ]], Cell[17560, 463, 39, 0, 29, "Text"], Cell[CellGroupData[{ Cell[17624, 467, 149, 3, 38, "Input"], Cell[17776, 472, 1214, 18, 132, "Output"] }, Open ]], Cell[19005, 493, 69, 0, 29, "Text"], Cell[CellGroupData[{ Cell[19099, 497, 129, 2, 58, "Input"], Cell[19231, 501, 1185, 19, 132, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[20453, 525, 66, 1, 38, "Input"], Cell[20522, 528, 180, 3, 22, "Message"], Cell[20705, 533, 1275, 21, 151, "Output"] }, Open ]], Cell[21995, 557, 495, 9, 71, "Text"], Cell[22493, 568, 1468, 24, 458, "Input"], Cell[23964, 594, 39, 0, 29, "Text"], Cell[CellGroupData[{ Cell[24028, 598, 134, 2, 38, "Input"], Cell[24165, 602, 233, 4, 39, "Print"] }, Open ]], Cell[24413, 609, 38, 0, 29, "Text"], Cell[CellGroupData[{ Cell[24476, 613, 134, 2, 38, "Input"], Cell[24613, 617, 233, 4, 39, "Print"] }, Open ]], Cell[24861, 624, 45, 0, 29, "Text"], Cell[CellGroupData[{ Cell[24931, 628, 133, 2, 38, "Input"], Cell[25067, 632, 231, 4, 39, "Print"] }, Open ]], Cell[25313, 639, 108, 3, 29, "Text"], Cell[CellGroupData[{ Cell[25446, 646, 152, 3, 38, "Input"], Cell[25601, 651, 38, 1, 37, "Output"] }, Open ]], Cell[25654, 655, 209, 6, 29, "Text"], Cell[25866, 663, 164, 3, 29, "Text"], Cell[26033, 668, 121, 2, 38, "Input"], Cell[26157, 672, 240, 4, 38, "Input"], Cell[CellGroupData[{ Cell[26422, 680, 175, 4, 38, "Input"], Cell[26600, 686, 109, 2, 50, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[26746, 693, 104, 2, 39, "Input"], Cell[26853, 697, 130, 2, 53, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[27020, 704, 73, 1, 38, "Input"], Cell[27096, 707, 186, 3, 37, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[27319, 715, 90, 1, 38, "Input"], Cell[27412, 718, 396, 7, 56, "Output"] }, Open ]], Cell[27823, 728, 106, 3, 29, "Text"], Cell[27932, 733, 266, 5, 78, "Input"], Cell[28201, 740, 112, 2, 38, "Input"] }, Closed]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)