Re: Wick like theorem and "symbolic" compilation

*To*: mathgroup at smc.vnet.net*Subject*: [mg61131] Re: Wick like theorem and "symbolic" compilation*From*: Peter Pein <petsie at dordos.net>*Date*: Tue, 11 Oct 2005 03:20:29 -0400 (EDT)*References*: <diaaps$ic7$1@smc.vnet.net> <did2ub$qet$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Peter Pein schrieb: > Alexander schrieb: > >>Dear MathGroup! >> >>I wrote simple programm to perform "Wick" expanding. >>Here is the code: >> >>(*St1 - basic function to extract couples of fields from a given list*) >> >>st1[x_] := If[Length[ >> x] === 2, Wick[Sequence @@ x], Plus @@ >> ReplaceList[List @@ (x), {e___, a_, c___, >> b_, d___} -> Wick[a, b] e c d]]; >> >>(* Some transformation rules for this function *) >> >>st1[x_Wick y_] := x st1[y]; >>st1[a_ + b_] := st1[a] + st1[b]; >> >>(* Formatting rules *) >> >>\!\(\(Format[Wick[x_, y_]] := \(x\ y\)\&^;\)\) >> >>(* Final function *) >> >>WickTheorem[x_] := Expand[Plus @@ FixedPointList[st1, x]] //. a_ + >> b_Integer c_ -> a + c; >> >>Simple example, showing how it works: >> >>In[6]:= WickTheorem[a b c d] >> >>Out[6]:= \!\(\*FormBox[ >> RowBox[{\(a\ b\ c\ d\), "+", >> RowBox[{"c", " ", >> FormBox[\(\(a\ b\)\&^\), >> "TraditionalForm"], " ", "d"}], "+", >> RowBox[{"b", " ", >> FormBox[\(\(a\ c\)\&^\), >> "TraditionalForm"], " ", "d"}], "+", >> RowBox[{"a", " ", >> FormBox[\(\(b\ c\)\&^\), >> "TraditionalForm"], " ", "d"}], "+", >> RowBox[{"b", " ", "c", " ", >> FormBox[\(\(a\ d\)\&^\), >> "TraditionalForm"]}], "+", >> RowBox[{ >> FormBox[\(\(a\ d\)\&^\), >> "TraditionalForm"], " ", >> FormBox[\(\(b\ c\)\&^\), >> "TraditionalForm"]}], "+", >> RowBox[{"a", " ", "c", " ", >> FormBox[\(\(b\ d\)\&^\), >> "TraditionalForm"]}], "+", >> RowBox[{ >> FormBox[\(\(a\ c\)\&^\), >> "TraditionalForm"], " ", >> FormBox[\(\(b\ d\)\&^\), >> "TraditionalForm"]}], "+", >> RowBox[{"a", " ", "b", " ", >> FormBox[\(\(c\ d\)\&^\), >> "TraditionalForm"]}], "+", >> RowBox[{ >> FormBox[\(\(a\ b\)\&^\), >> "TraditionalForm"], " ", >> FormBox[\(\(c\ d\)\&^\), >> "TraditionalForm"]}]}], TraditionalForm]\) >> >>Maybe there is more simple and elegant realization of this theorem than >>my code ??? >> >>Second question is about optimization. >> >>In current realization, on my Cleleron 1700 (with 128Mb of ram) using >>Mathematica 5.1 for Windows, I have following timings: >> >>In[7]:= WickTheorem[a b c d e f w t] // Timing // Print["Time = ", #[[ >> 1]], ", Length = ", Length[#[[2]]]] & >>Out[7]:= Time = 10.688 Second, Length = 764 >> >>So, in final expression there was 764 terms and it takes ~10 seconds to >>evaluate them. Input expression contains only 8 fields, with 9 fields >>computation takes a relatively long time,about 600 seconds. >> >>Why so long? >>And what are the ways to reduce this time ? >> >>I have tried to use CompileEvaluate from Experimental package, but this >>doesn't reduce evaluation time greatly. >> >>Resuming, second question sounds like: >>Should I rewrite the program in more elegant and efficient way to >>reduce >>evaluation time or there is "symbolic" compilation technique to do >>that? >> >>Thanks for your help! >> >>Alexander. >> > > Hi Alexander, > > there's only one rule necessary, because Wick[x,y] is not a symbol. > > Off[General::spell1]; > wickrule={c___,a_Symbol,d___, b_Symbol ,e___}->{ c, d ,e,Wick[a,b]}; > WickTheorem[p_Times]:= > (Plus@@Times@@@Union@@ > NestWhileList[Flatten[(ReplaceList[#,wickrule]&)/@#,1]&, > {List@@p},{}=!=#&] > )/._Integer->1 > Format[Wick[x_,y_]] := OverHat[x*y]; > > In[5]:= WickTheorem[a*b*c*d] > Out[5]= > a*b*c*d + c*d*a^b + b*d*a^c + b*c*a^d + a*d*b^c + > a^d*b^c + a*c*b^d + a^c*b^d + a*b*c^d + a^b*c^d > > On an Athlon64 3000+ at 2130 MHz under Win2k with Mathematica 5.1 and > 512MB RAM, I get the following timings: > > In[6]:= > (Timing[Length[ > WickTheorem[Times @@ ToExpression["a"<>ToString[#]]&/@ Range[#]]]]&) > /@ Range[10] > Out[6]= > {{0.*Second, 1}, {0.*Second, 2}, {0.*Second, 4}, {0.*Second, 10}, > {0.*Second, 26}, {0.*Second, 76}, {0.016*Second, 232}, > {0.125*Second, 764}, {0.937*Second, 2620}, {8.469*Second, 9496}} > > I tried a rule that transforms products (a_Symbol b_Symbol > c_./;OrderedQ[a,b]->Wick[a,b] c), but it has not been as fast as this one. > > If you are just interested in the number of the resulting terms, > Sum[(2k)!/(k! 2^k) * Binomial[n,2k],{k,0,n/2}] will give the answer > almost intantly (see > http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=A000085). > > > Peter > For more than 6 fields the straightforward approach seems to be more efficient: Clear[st]; st[p:_Symbol | _Wick | _Integer] := p; st[p_Plus] := st[p] = st /@ p; st[p_Times] := st[p] = Module[{i, j}, Expand[p + Sum[Wick @@ p[[{i, j}]]*st[Delete[p, {{i}, {j}}]], {i, 1, Length[p] - 1}, {j, i + 1, Length[p]}]] /. _Integer -> 1] st[a*b*c*d] a*b*c*d + c*d*Wick[a, b] + b*d*Wick[a, c] + b*c*Wick[a, d] + a*d*Wick[b, c] + Wick[a, d]*Wick[b, c] + a*c*Wick[b, d] + Wick[a, c]*Wick[b, d] + a*b*Wick[c, d] + Wick[a, b]*Wick[c, d] Timing[Length[st[a*b*c*d*e*f*g*h*i*j]]] {2.282*Second, 9496} Peter