home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / COMPLEXE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  40.8 KB  |  1,241 lines
  1.  
  2. (* ****************************************************************
  3. *
  4. *     ComplexExpand, Roman Maeder, V2.0
  5. *     Debugged/enhanced for V2.1 by Daniel Lichtblau
  6. ***************************************************************** *)
  7.  
  8. BeginPackage["System`ComplexExpand`"]
  9.  
  10. (*
  11. This is not too dangerous since the System`ComplexExpand` context is removed at
  12. the end
  13. *)
  14.  
  15. Off[ General::shdw]
  16.  
  17. Unprotect[ {ComplexExpand, TargetFunctions, System`ComplexExpand`AbsExpr,
  18.             System`ComplexExpand`ArgExpr, System`ComplexExpand`ConjugateExpr,
  19.             System`ComplexExpand`ReImExpr, System`ComplexExpand`SignExpr}]
  20.  
  21.  
  22. (********************************************************************** *)
  23.  
  24. (* ****************************************************************
  25. *
  26. *     Rigorous Definitions Of The Elementary Inverse Functions
  27. *
  28. *Definitions:
  29. *
  30. *(1) Log[z_]:=If[And[SameQ[Im[z],0],Re[z]<0],Log[-z]+Pi*I,Log[z]]
  31. *(2) Power[u_,v_]:=E^(v*Log[u])
  32. *(3) Sqrt[z_]:=Power[z,1/2]
  33. *(4) ArcSin[z_]:=-I*Log[I*z+Sqrt[1-z^2]]
  34. *(5) ArcCos[z_]:=Pi/2-ArcSin[z]
  35. *(6) ArcTan[z_]:=-I*Log[(1+I*z)/Sqrt[1+z^2]]
  36. *    ArcTan[u_,v_]:=-I*Log[(u+I*v)/Sqrt[u^2+v^2]]
  37. *(7) ArcCsc[z_]:=ArcSin[1/z]
  38. *(8) ArcSec[z_]:=ArcCos[1/z]
  39. *(9) ArcCot[z_]:=If[SameQ[z,0],Pi/2,ArcTan[1/z]]
  40. *(10) ArcSinh[z_]:=Log[z+Sqrt[1+z^2]]
  41. *(11) ArcCosh[z_]:=Log[z+(z+1)*Sqrt[(z-1)/(z+1)]]
  42. *(12) ArcTanh[z_]:=Log[(1+z)/Sqrt[1-z^2]]
  43. *(13) ArcCsch[z_]:=ArcSinh[1/z]
  44. *(14) ArcSech[z_]:=ArcCosh[1/z]
  45. *(15) ArcCoth[z_]:=If[SameQ[z,0],Pi*I/2,ArcTanh[1/z]]
  46. *
  47. *Branch Cuts:
  48. *     
  49. *(1) Log
  50. *     (-Infinity,0) attaches to quadrant II.
  51. *     Im[Log[z]] == Pi     for z in (-Infinity,0)
  52. *(2) Power  (Considering Power for constant v.)
  53. *     (-Infinity,0) attaches to quadrant II, for u not an integer.
  54. *     Power[u,v] == (-u)^v*(Cos[v]+I*Sin[v])     for v in (-Infinity,0)
  55. *(3) Sqrt
  56. *     (-Infinity,0) attaches to quadrant II.
  57. *     Sign[Im[Sqrt[z]]] > 0     for z in (-Infinity,0)
  58. *(4) ArcSin
  59. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  60. *     Re[ArcSin[z]] == -Pi/2     for z in (-Infinity,-1)
  61. *     Re[ArcSin[z]] == Pi/2      for z in (1,Infinity)
  62. *(5) ArcCos
  63. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  64. *     Re[ArcCos[z]] == Pi     for z in (-Infinity,-1)
  65. *     Re[ArcCos[z]] == 0      for z in (1,Infinity)
  66. *(6) ArcTan
  67. *     (-I*Infinity,-I) attaches to quadrant IV and (I,I*Infinity) to quadrant II.
  68. *     Re[ArcTan[z]] == -Pi/2     for z in (-I*Infinity,-I)
  69. *     Re[ArcTan[z]] == Pi/2      for z in (I,I*Infinity)
  70. *(7) ArcCsc
  71. *     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
  72. *     Re[ArcCsc[z]] == -Pi/2     for z in (-1,0).
  73. *     Re[ArcCsc[z]] == Pi/2      for z in (0,1).
  74. *(8) ArcSec
  75. *     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
  76. *     Re[ArcSec[z]] == Pi     for z in (-1,0).
  77. *     Re[ArcSec[z]] == 0      for z in (0,1).
  78. *(9) ArcCot
  79. *     (-I,0) attaches to quadrant III and (0,I) to quadrant I.
  80. *     Re[ArcCot[z]] == Pi/2     for z in (-I,0).
  81. *     Re[ArcCot[z]] == -Pi/2    for z in (0,I).
  82. *(10) ArcSinh
  83. *     (-I*Infinity,-I) attaches to quadrant I and (I,I*Infinity) to quadrant III.
  84. *     Im[ArcSinh[z]] == -Pi/2   for z in (-I*Infinity,-I)
  85. *     Im[ArcSinh[z]] == Pi/2    for z in (I,I*Infinity)
  86. *(11) ArcCosh
  87. *     (-Infinity,1) attaches to quadrants I and II.
  88. *     Sign[Im[ArcCosh[z]]] > 0     for z in (-Infinity,1).
  89. *(12) ArcTanh
  90. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  91. *     Im[ArcTanh[z]]] == Pi/2      for z in (-Infinity,-1).
  92. *     Im[ArcTanh[z]]] == -Pi/2     for z in (1,Infinity).
  93. *(13) ArcCsch
  94. *     (-I,0) attaches to quadrant IV and (0,I) to quadrant I.
  95. *     Sign[Re[ArcCsch[z]]] > 0     for z in (-I,0)
  96. *     Sign[Re[ArcCsch[z]]] < 0     for z in (0,I)
  97. *(14) ArcSech
  98. *     (-Infinity,0) attaches to quadrant III and (1,Infinity) to quadrant IV.
  99. *     Sign[Im[ArcSech[z]]] > 0     for z in (-Infinity,0)
  100. *     Sign[Im[ArcSech[z]]] > 0     for z in (1,Infinity)
  101. *(15) ArcCoth
  102. *     (-1,0] attaches to quadrant III and (0,1) to quadrant I.
  103. *     Im[ArcCoth[z]] == Pi/2        for z in (-1,0]
  104. *     Im[ArcCoth[z]] == -Pi/2      for z in (0,1)
  105. ***************************************************************** *)
  106.  
  107. (* ****************************************************************
  108. *
  109. *     Useful Trigonometric Identities
  110. *
  111. *    Sinh[X] == -I*Sin[I*X] == -Sinh[-X] == I*Sin[-I*X]
  112. *    Cosh[X] == Cos[I*X] == Cosh[-X] == Cos[-I*X]
  113. *    Tanh[X] == -I*Tan[I*X] == -Tanh[-X] == I*Tan[-I*X]
  114. *    ArcSinh[X] == I*ArcSin[-I*X] == Log[X+Sqrt[1+X^2]]
  115. *    Sin[X] == -I*Sinh[I*X] == -Sin[-X] == I*Sinh[-I*X]
  116. *    Cos[X] == Cosh[I*X] == Cos[-X] == Cosh[-I*X]    
  117. *    Tan[X] == -I*Tanh[I*X] == -Tan[-X] == I*Tanh[-I*X]
  118. *    ArcSin[X] == -I*ArcSinh[I*X] == -I*Log[I*X+Sqrt[1-X^2]]
  119. *
  120. *Note: Be very careful with ArcSinh and ArcSin.  The above equations
  121. *are correct.  Cannot assume X is real and must pay attention to branch
  122. *cuts.
  123. ***************************************************************** *)
  124.  
  125. (* ****************************************************************
  126. *
  127. *     Abs, Arg, Conjugate, Im, Re, Sign Identities
  128. *
  129. *For X!=0:
  130. *
  131. *     Abs[X] == Sqrt[X*Conjugate[X]]
  132. *     Arg[X] == -I*Log[X/Sqrt[X*Conjugate[X]]]
  133. *     Im[X] == (X-Conjugate[X])/2
  134. *     Re[X] == (X+Conjugate[X])/2
  135. *     Sign[X] == X/Sqrt[X*Conjugate[X]]
  136. *
  137. *All other relations follow from these relations.  Examples:
  138. *
  139. *     Conjugate[X]=Abs[X]^2/X
  140. *     Re[X]=X/2+X/(2*Sign[X]^2)
  141. *     Im[X]=X*(1-E^(-2*I*Arg[X]))/2
  142. *
  143. *Etc.
  144. ***************************************************************** *)
  145.  
  146. (* ****************************************************************
  147. *
  148. *     Leaf Functions
  149. *
  150. ***************************************************************** *)
  151.  
  152. LeafAbs[f_?NumberQ, True] := Abs[f] /; IsOK[Abs,f]
  153. LeafAbs[f_?ConstantQ, toplevel_] := ConstantAbs[f, toplevel]
  154.  
  155. LeafAbs[f_, True] := ReImFail /; IsOK[Abs,f]
  156. (* LeafAbs[f_, True] :=Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im},f] *)
  157.  
  158. LeafAbs[f_?RealValQ, False] := Sqrt[f^2]
  159.  
  160. LeafAbs[f_, False] := Abs[f] /; IsOK[Abs,f]
  161.  
  162. (* cartesian *)
  163. LeafAbs[f_, False] := Module[{x,y,norm},
  164.     {x,y}=ReImExpr[f,False];
  165.     norm=GetNorm[x,y];
  166.     Which[norm===1,1,
  167.         True, Sqrt[x^2 + y^2]]
  168.     ] /; IsOK[{Re, Im},f]
  169.     
  170. LeafAbs[f_, True]:= Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im},f]
  171.  
  172. (* polar: n/a *)
  173.  
  174. (* other cases *)
  175. LeafAbs[f_, False]:=
  176.     Which[  IsOK[Arg, f], f/E^(I*ArgExpr[f, False]),
  177.         IsOK[Conjugate, f], Sqrt[f*ConjugateExpr[f, False]],
  178.         IsOK[Im, f], Sqrt[f^2-2*I*f*ReImExpr[f,False][[2]]],
  179.         IsOK[Re, f], Sqrt[2*f*ReImExpr[f, False][[1]]-f^2],
  180.         IsOK[Sign, f], f/SignExpr[f, False],
  181.         True, ReImFail]
  182.  
  183. LeafAbs[f_, True]:=
  184.     Which[    IsOK[Arg, f], f/E^(I*Arg[f]),
  185.         IsOK[Conjugate, f], Sqrt[f*Conjugate[f]],
  186.         IsOK[Im, f], Sqrt[2*f*Re[f]-f^2],
  187.         IsOK[Re, f], Sqrt[f^2-2*I*f*Im[f]],
  188.         IsOK[Sign, f], f/Sign[f],
  189.         True, ReImFail]
  190.  
  191.  
  192. LeafArg[f_?NumberQ, _] := Arg[f]
  193.  
  194. LeafArg[f_?ConstantQ, toplevel_] := ConstantArg[f, toplevel]
  195.  
  196. LeafArg[f_, True]:= ReImFail /; IsOK[Arg,f]
  197.  
  198. LeafArg[f_, False]:= Arg[f] /; IsOK[Arg,f] 
  199.  
  200. (* cartesian *)
  201. LeafArg[f_, True]:= ArcTan[Re[f],Im[f]] /; IsOK[{Re, Im},f]
  202. LeafArg[f_, False]:= 
  203.   Module[{x,y,q},
  204.     {x,y}=ReImExpr[f, False];
  205.     (*q=y/x;
  206.     If [Head[q]===Tan && N[-Pi < N[q[[1]]] <= Pi], q[[1]],]*)     
  207.     MyArcTan[x, y, False]
  208.        ]  /; IsOK[{Re, Im},f] 
  209.  
  210. (* polar: n/a *)
  211.  
  212. (* other cases *)
  213. LeafArg[f_, False] := 
  214.     Which[  IsOK[Im, f], With[{im=ReImExpr[f, False][[2]]},
  215.             MyArcTan[f-I*im, im, False]],
  216.         IsOK[Re, f], With[{re=ReImExpr[f, False][[1]]},
  217.             MyArcTan[re, -I*(f-re), False]],
  218.         IsOK[Sign, f], -I*Log[SignExpr[f, False]],
  219.         IsOK[Abs, f], -I*Log[f/AbsExpr[f, False]],
  220.         IsOK[Conjugate, f], -I*Log[f/Sqrt[f*ConjugateExpr[f, False]]],
  221.         True, ReImFail]
  222.  
  223. LeafArg[f_, True] :=
  224.     Which[  IsOK[Im, f], ArcTan[f-I*Im[f],Im[f]],
  225.         IsOK[Re, f], ArcTan[Re[f],-I*(f-Re[f])],
  226.         IsOK[Sign, f], -I*Log[Sign[f]],
  227.         IsOK[Abs, f], -I*Log[f/Abs[f]],
  228.         IsOK[Conjugate, f], -I*Log[f/Sqrt[f*Conjugate[f]]],
  229.         True, ReImFail]
  230.  
  231.  
  232. LeafConjugate[f_?RealValQ, False] := f
  233.  
  234. LeafConjugate[f_?NumberQ, _] := Conjugate[f]
  235. LeafConjugate[g_?ConstantQ, True] := Re[g] - I Im[g]
  236.  
  237. LeafConjugate[f_, True]:= ReImFail /; IsOK[Conjugate,f]
  238. LeafConjugate[f_, False]:= Conjugate[f] /; IsOK[Conjugate,f]
  239.  
  240. (* cartesian *)
  241. LeafConjugate[f_, True]:= Re[f]-I*Im[f] /; IsOK[{Re, Im}]
  242. LeafConjugate[f_, False]:= Module[{x,y},
  243.     {x,y} = ReImExpr[f, False];
  244.     x-I*y ]  /; IsOK[{Re, Im},f] 
  245.  
  246. (* polar *)
  247. LeafConjugate[f_, True]:= Abs[f] Exp[-I Arg[f]] /; IsOK[{Arg, Abs}]
  248. LeafConjugate[f_, False]:= AbsExpr[f, False] * 
  249.     Exp[-I ArgExpr[f, False]] /; IsOK[{Arg, Abs},f]
  250.  
  251. (* other cases *)
  252. LeafConjugate[f_, False]:=
  253.     Which[  IsOK[Arg, f], f*E^(-2*I*ArgExpr[f, False]),
  254.         IsOK[Im, f], f-2*I*ReImExpr[f, False][[2]],
  255.         IsOK[Re, f], 2*ReImExpr[f, False][[1]]-f,
  256.         IsOK[Abs, f], AbsExpr[f, False]^2/f,
  257.         IsOK[Sign, f], f*SignExpr[f, False]^(-2),
  258.         True, ReImFail]
  259.  
  260. LeafConjugate[f_, True]:=
  261.     Which[  IsOK[Arg, f], f*E^(-2*I*Arg[f]),
  262.         IsOK[Im, f], f-2*I*Im[f],
  263.         IsOK[Re, f], 2*Re[f]-f,
  264.         IsOK[Abs, f], Abs[f]^2/f,
  265.         IsOK[Sign, f], f*Sign[f]^(-2),
  266.         True, ReImFail]
  267.  
  268.  
  269. LeafReIm[f_?NumberQ, _] := {Re[f], Im[f]}
  270. LeafReIm[f_?ConstSymbQ, _] := {Re[f], Im[f]}
  271.  
  272. LeafReIm[f_, True]:= {ReImFail, ReImFail} /; IsOK[{Re, Im},f]
  273. LeafReIm[f_, True]:= {ReImFail, -I*(f-Re[f])} /; IsOK[Re,f]
  274. LeafReIm[f_, True]:= {f-I*Im[f], ReImFail} /; IsOK[Im,f]
  275.  
  276. (* cartesian *)
  277.  
  278. LeafReIm[f_?RealValQ, False]:= {f, 0}
  279. LeafReIm[f_?((I * RealValQ) || (RealValQ * I)), False]:= {0, f}
  280.  
  281. LeafReIm[f_, False]:= {Re[f], Im[f]} /; IsOK[{Re, Im},f]
  282.  
  283. (* mixed mode, evaluate only once *)
  284. LeafReIm[f_, False]:= With[{re = Re[f]}, {re, -I*(f-re)}] /;
  285. IsOK[Re,f]
  286. LeafReIm[f_, False]:= With[{im = Im[f]}, {f-I*im, im}] /;
  287. IsOK[Im,f]
  288.  
  289. (* polar *)
  290. LeafReIm[f_, False]:=
  291.     With[{arg = ArgExpr[f, False], abs = AbsExpr[f, False]},
  292.         {abs Cos[arg], abs Sin[arg]}
  293.     ] /; IsOK[{Arg, Abs},f]
  294. LeafReIm[f_, True]:=
  295.     With[{arg = Arg[f], abs = Abs[f]},
  296.         {abs Cos[arg], abs Sin[arg]}
  297.     ] /; IsOK[{Arg, Abs},f]
  298.  
  299. (* other cases *)
  300. LeafReIm[f_, False]:=
  301.     Which[     IsOK[Conjugate, f], With[{conj = ConjugateExpr[f, False]},
  302.               {(f+conj)/2, (f-conj)/(2*I)}],
  303.         IsOK[Arg, f], With[{arg = ArgExpr[f, False]},
  304.               {f/2+f*E^(-2*I*arg)/2,
  305.                f/(2*I)-f*E^(-2*I*arg)/(2*I)}],
  306.         IsOK[Abs, f], With[{abs = AbsExpr[f, False]},
  307.               {f/2+abs^2/(2*f), f/(2*I)-abs^2/(2*I*f)}],
  308.         IsOK[Sign, f], With[{sign = SignExpr[f, False]},
  309.               {f/2+f*sign^(-2)/2,
  310.                f/(2*I)-f*sign^(-2)/(2*I)}],
  311.         True, ReImFail]
  312.  
  313. LeafReIm[f_, True]:=
  314.     Which[  IsOK[Conjugate, f], With[{con = Conjugate[f]},
  315.               {(f+con)/2, (f-con)/(2*I)}],
  316.         IsOK[Arg, f], With[{arg = Arg[f]},
  317.               {f/2+f*E^(-2*I*Arg[f])/2,
  318.                f/(2*I)-f*E^(-2*I*Arg[f])/(2*I)}],
  319.         IsOK[Abs, f], With[{abs = Abs[f]},
  320.               {f/2+abs^2/(2*f), f/(2*I)-abs^2/(2*I*f)}],
  321.         IsOK[Sign, f], With[{sign = Sign[f]},
  322.               {f/2+f*sign^(-2)/2,
  323.                f/(2*I)-f*sign^(-2)/(2*I)}],
  324.         True, ReImFail]
  325.  
  326.  
  327. LeafSign[f_?NumberQ, True] := Sign[f]
  328. (*LeafSign[f_?ConstSymbQ, True] := Sign[f]*)
  329.  
  330. LeafSign[f_?ConstantQ, toplevel_] := ConstantSign[f, toplevel]
  331.  
  332. LeafSign[f_, True]:= ReImFail /; IsOK[Sign,f]
  333. LeafSign[f_, False]:= Sign[f] /; IsOK[Sign,f]
  334.  
  335. (* cartesian *)
  336. LeafSign[f_, True]:= (Re[f] + I Im[f])/Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im}]
  337. LeafSign[f_, False]:= Module[{x,y,norm,normsqrt}, 
  338.             {x,y}=ReImExpr[f, False];
  339.             norm=GetNorm[x,y];
  340.             If[norm===1,normsqrt=1,normsqrt=Sqrt[norm]];
  341.             (x + I y)/normsqrt
  342.         ] /; IsOK[{Re, Im},f]
  343.  
  344. (* other cases *)
  345. LeafSign[f_, False]:=
  346.     Which[  IsOK[Arg, f], Exp[I ArgExpr[f, False]],
  347.         IsOK[Abs, f], f/AbsExpr[f, False],
  348.         IsOK[Conjugate, f], 
  349.             f/Sqrt[f*ConjugateExpr[f, False]],
  350.         IsOK[Re, f], f/Sqrt[2*f*ReImExpr[f, False][[1]]-f^2],
  351.         IsOK[Im, f], f/Sqrt[f^2-2*I*f*ReImExpr[f, False][[2]]],
  352.         True, ReImFail] 
  353.  
  354. LeafSign[f_, True]:=
  355.     Which[  IsOK[Arg, f], Exp[I Arg[f]],
  356.         IsOK[Abs, f], f/Abs[f],
  357.         IsOK[Conjugate, f], f/Sqrt[f*Conjugate[f]],
  358.         IsOK[Re, f], f/Sqrt[2*f*Re[f]-f^2],
  359.         IsOK[Im, f], f/Sqrt[f^2-2*I*f*Im[f]],
  360.         True, ReImFail]
  361.  
  362.  
  363. (* End of leaf code ********************************************** *)
  364.  
  365. (* ****************************************************************
  366. *
  367. *     Arg
  368. *
  369. ***************************************************************** *)
  370.  
  371. ArgExpr[x_?NonNegValQ, False] := 0
  372. ArgExpr[Abs[_], _] := 0
  373. ArgExpr[Arg[x_?RealValQ], False] := 0
  374. ArgExpr[Arg[Arg[_]], False] := 0
  375. ArgExpr[Sign[x_], False] := ArgExpr[x, False]
  376.  
  377. ArgExpr[Re[x_], False] := With[{re=ReImExpr[x, False][[1]]},
  378.     Which[re===Re[x], ArgRecursive[re, False],
  379.     True, ArgExpr[re, False] ] ]
  380.  
  381. ArgExpr[Im[x_], False] := With[{im=ReImExpr[x, False][[2]]},
  382.     Which[im===Im[x], ArgRecursive[im, False],
  383.     True, ArgExpr[im, False] ] ]
  384.  
  385. ArgExpr[Arg[x_], False] := With[{arg=ArgExpr[x, False]},
  386.     Which[arg===Arg[x], ArgRecursive[arg, False],
  387.     True, ArgExpr[arg, False] ] ]
  388.  
  389. ArgExpr[Conjugate[x_], False] := With[{conj=ConjugateExpr[x, False]},
  390.     Which[conj===Conjugate[x], ArgRecursive[conj, False],
  391.     True, ArgExpr[conj, False] ] ]
  392.  
  393. ArgExpr[f_?ConstantQ, False] := ConstantArg[f, False]
  394.  
  395. ArgExpr[e_, True] := (ArgRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
  396.     IsRecursive[Arg] || ConstantQ[e]
  397. ArgExpr[e_, False] := ArgRecursive[e, False] /;
  398.     IsRecursive[Arg] || ConstantQ[e]
  399.  
  400. ArgExpr[e_, True] := LeafArg[e, True] /. ArcTanTemp -> ArcTan
  401. ArgExpr[e_, False] := LeafArg[e, False]
  402.  
  403. ArgRecursive[e_ ^ r_Rational, toplevel_] :=
  404.     r ArgExpr[e,toplevel] /; r < 1 (* principal root *)
  405.  
  406. ArgRecursive[e_ ^ (_Rational|_Integer), False] := 0 /;
  407.      ArgExpr[e, toplevel] == 0
  408. ArgRecursive[e_ ^ (_Rational|_Integer), True] := 0 /;
  409.      Arg[e] == 0
  410.  
  411. ArgRecursive[e_?RealValQ ^ p_Integer, False] := 0 /; EvenQ[p]
  412.  
  413. (* internal code makes -1 < Im[c] <= 1 *)
  414. ArgRecursive[ Exp[c_Complex Pi], _ ] := Im[c] Pi
  415. ArgRecursive[ Exp[e_], _ ] := 0 /; Im[e] === 0 (* Exp[e] positive for real e *)
  416.  
  417. ArgRecursive[e_ ^ p_, toplevel_:False] := 
  418.     Module[{abs,arg,x,y,theta},
  419.         abs=AbsExpr[e,toplevel];
  420.         arg=ArgExpr[e,toplevel];
  421.         {x,y}=ReImExpr[p,toplevel];
  422.         theta = y * RealLog[abs] + x * arg;
  423.         MyArcTan[Cos[theta], Sin[theta], toplevel]
  424.     ]
  425.  
  426. (* special exponentials *)
  427.  
  428. ArgRecursive[e_, toplevel_] := LeafArg[e, toplevel]
  429.  
  430. (* ****************************************************************
  431. *
  432. *     Abs
  433. *
  434. ***************************************************************** *)
  435.  
  436. AbsExpr[t:Abs[_], True] := t
  437.  
  438. AbsExpr[Abs[x_], False] := AbsExpr[x,False]
  439.  
  440. AbsExpr[Re[x_], False] := With[{re=ReImExpr[x, False][[1]]},
  441.     Which[re===Re[x], AbsRecursive[re, False],
  442.     True, AbsExpr[re, False] ] ]
  443.  
  444. AbsExpr[Im[x_], False] := With[{im=ReImExpr[x, False][[2]]},
  445.     Which[im===Im[x], AbsRecursive[im, False],
  446.     True, AbsExpr[im, False] ] ]
  447.  
  448. AbsExpr[Sign[e_], False] := Which[
  449.     ZeroQ[e], 0,
  450.     PosValQ[e] || NegValQ[e], 1,
  451.     True, With[{sign=SignExpr[e, False]},
  452.         If[sign===Sign[e],AbsRecursive[sign, False],
  453.             AbsExpr[sign, False]] ] ]    
  454.  
  455. AbsExpr[x_?NonNegValQ, False] := x
  456.  
  457. AbsExpr[Conjugate[e_], toplevel_] := AbsExpr[e, toplevel]
  458. AbsExpr[Literal[Arg[e_]], False] := Sqrt[ArgExpr[e, False]^2]
  459.  
  460. AbsExpr[e_, True] := (AbsRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
  461.     IsRecursive[Abs] || ConstantQ[e]
  462. AbsExpr[e_, False] := AbsRecursive[e, False] /;
  463.     IsRecursive[Abs] || ConstantQ[e]
  464.  
  465. AbsExpr[e_, True] := LeafAbs[e, True] /. ArcTanTemp -> ArcTan
  466. AbsExpr[e_, False] := LeafAbs[e, False]
  467.  
  468. AbsRecursive[e_Plus, toplevel_] := AbsPlus[e, toplevel]
  469. AbsRecursive[e_Times, toplevel_] :=  AbsExpr[#, False]& /@ e
  470.  
  471. AbsRecursive[g_ ^ n_Integer, toplevel_] := 
  472.     Module[{x, y},
  473.         {x,y} = ReImExpr[g, toplevel];
  474.         Which[    y === 0, If[EvenQ[n], g^n, AbsExpr[g, toplevel]^n],
  475.             x === 0, If[Mod[n,4] === 0, g^n,
  476.                 AbsExpr[g, toplevel]^n],
  477.             True, AbsExpr[g, toplevel]^n]
  478.            ]
  479.  
  480. AbsRecursive[t:g_ ^ h_, toplevel_] := 
  481.     Module[{x, y}, 
  482.         {x, y} = ReImExpr[h, toplevel];
  483.         AbsExpr[g, toplevel]^x E^(-y*ArgExpr[g, toplevel])
  484.            ]
  485.  
  486. AbsRecursive[e_, toplevel_] := LeafAbs[e, toplevel]
  487.  
  488. AbsPlus[f_Plus, toplevel_] := AbsExpr[#, False]& /@ f /;
  489.     SameSignsQ[f, toplevel] =!= ReImFail
  490.  
  491. AbsPlus[e_, toplevel_] := LeafAbs[e, toplevel]
  492.  
  493. (* ****************************************************************
  494. *
  495. *     ConjugateExpr
  496. *
  497. ***************************************************************** *)
  498.  
  499. (* commute with Conjugate -> real valued for real arg *)
  500.  
  501. ConjugateFunctions={
  502.   ArcCot,  ArcCsch, ArcSinh, ArcTan, 
  503.   Cos, Cosh, Cot, Coth, Csc, Csch,
  504.   Sec, Sech, Sin, Sinh, Tan, Tanh,
  505.   Plus, Times, Sign(* Sign --Handle this a bit differently in ReImExpr *)
  506. }
  507.  
  508. (* are real valued always: *)
  509.  
  510. RealFunctions = { Abs, Arg, Im, Re }
  511.  
  512. ConjugateExpr[z_?(!ConstantQ), toplevel_] := Module[{x,y}, 
  513.     {x,y} =  ReImExpr[z, toplevel]; 
  514.     Which[ (x===ReImFail) || (y===ReImFail), ReImFail,
  515.     True, x - I y] ] /; PolarQ[z] || CartQ[z] || ! IsOK[Conjugate, z]
  516.  
  517. ConjugateExpr[Conjugate[f_], _] := f
  518.  
  519. ConjugateExpr[t:f_Symbol[_], True] := t /; MemberQ[RealFunctions,f]
  520. ConjugateExpr[t:f_Symbol[_], False] := ReImExpr[t, False][[1]] /;
  521.     MemberQ[RealFunctions,f]
  522.  
  523. ConjugateExpr[(-1) ^ y_?RealValQ,_] := (-1)^(-y)
  524. ConjugateExpr[x_?NegValQ ^ y_?RealValQ,_] := (-1)^(-y) * (-x)^y
  525. ConjugateExpr[x_?PosValQ ^ Complex[0,z_], _] := x^(-Complex[0,z])
  526. ConjugateExpr[x_?PosValQ ^ (y_?RealValQ * Complex[0,z_]), _] :=
  527.     x^(-y*Complex[0,z])
  528.  
  529. ConjugateExpr[e_, True] :=
  530.     (ConjugateRecursive[e, True] /.ArcTanTemp -> ArcTan) /;
  531.         IsRecursive[Conjugate] || ConstantQ[e]
  532. ConjugateExpr[e_, False] := ConjugateRecursive[e, False] /;
  533.     IsRecursive[Conjugate] || ConstantQ[e]
  534.     
  535. ConjugateExpr[e_, True] := LeafConjugate[e, True] /. ArcTanTemp -> ArcTan
  536. ConjugateExpr[e_, False] := LeafConjugate[e, False]
  537.  
  538. ConjugateRecursive[t:f_Symbol[__], True] := Conjugate /@ t /; 
  539.         MemberQ[ConjugateFunctions,f] && ConstantQ[t]
  540.  
  541. ConjugateRecursive[t:f_Symbol[__],toplevel_] := 
  542.     ConjugateExpr[#, False]& /@ t /;
  543.         MemberQ[ConjugateFunctions,f] && f =!= Sign && !ConstantQ[t]
  544.  
  545. (* real powers *)
  546. ConjugateRecursive[e_?PosValQ ^ (r_Rational|r_Integer), _] := e^r
  547.  
  548. ConjugateRecursive[e_?RealValQ, toplevel_] := e
  549.  
  550. ConjugateRecursive[e_, toplevel_] := LeafConjugate[e, toplevel]
  551.  
  552.  
  553. (* ****************************************************************
  554. *
  555. *     SignExpr
  556. *
  557. ***************************************************************** *)
  558.  
  559. SignExpr[t:Sign[_], True] := t
  560. SignExpr[Sign[x_], False] := SignExpr[x,False]
  561.  
  562. SignExpr[t:f_Symbol[_], False] := With[{x=ReImExpr[t,
  563. False][[1]]},
  564.     Which[ZeroQ[x], 0,
  565.     PosValQ[x], 1,
  566.     NegValQ[x], -1,
  567.     x===t, SignRecursive[x, False],
  568.     True, SignExpr[x, False]
  569.      ] ] /; MemberQ[RealFunctions,f]
  570.  
  571. SignExpr[e_, True] := (SignRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
  572.     IsRecursive[Sign] || ConstantQ[e]
  573.  
  574. SignExpr[e_, False] := SignRecursive[e, False] /;
  575.     IsRecursive[Sign] || ConstantQ[e]
  576.  
  577. SignExpr[e_, True] := LeafSign[e, True] /. ArcTanTemp -> ArcTan
  578. SignExpr[e_, False] := LeafSign[e, False]
  579.  
  580. SignRecursive[e_Times, False] := SignExpr[#, False]& /@ e
  581. SignRecursive[e_Times, True] := Sign /@ e
  582.  
  583. SignRecursive[e_Plus, toplevel_] :=
  584.     Module[{answer = SameSignsQ[e, toplevel]},
  585.         If[ answer =!= ReImFail, answer, LeafSign[e, toplevel] ]
  586.     ] 
  587.  
  588. SignRecursive[g_ ^ n_Integer, True]:=
  589.     Module[{x,y},
  590.         {x,y}=ReImExpr[g, False];
  591.         Which[    y === 0, If[EvenQ[n], Sign[g]^2, Sign[g]],
  592.             x === 0, Switch[Mod[n,4],
  593.                               0, Sign[g]^2,
  594.                       1, I*Sign[g],
  595.                       2, -(Sign[g]^2),
  596.                       _, -I*Sign[g]],
  597.             True, Sign[g]^n]
  598.     ]
  599.     
  600. SignRecursive[g_ ^ n_Integer, False]:=
  601.     Module[{x,y},
  602.         {x,y}=ReImExpr[g,False];
  603.         Which[    y === 0, If [EvenQ[n],
  604.             SignExpr[g, False]^2, SignExpr[g, False]],
  605.             x === 0, Switch[Mod[n,4],
  606.                               0, SignExpr[g, False]^2,
  607.                       1, I*SignExpr[g, False],
  608.                       2, -(SignExpr[g, False]^2),
  609.                       _, -I*SignExpr[g, False]], 
  610.             True, SignExpr[g, False] ^ n]
  611.     ]
  612.  
  613. (* number roots *)
  614.  
  615. SignRecursive[g_ ^ r_?RealValQ, toplevel_] :=
  616.     Which[toplevel, Sign[g]^r,
  617.     True, SignExpr[g, toplevel] ^ r]
  618.     
  619. (* exponentials *)
  620.  
  621. SignRecursive[ Exp[e_], _ ] := Exp[ I Im[e] ]
  622.  
  623. SignRecursive[e_, toplevel_] := LeafSign[e, toplevel]
  624.  
  625. SameSignsQ[f_, toplevel_]:=
  626.   (* If all args have same sign (skipping 0), then return sign. *)
  627.     Module[{answer = 0, sign, i = 1, n = Length[f]},
  628.         While[i<=n,
  629.             (*sign=SignExpr[f[[i]], toplevel];*)
  630.             sign=Sign[f[[i]]];
  631.             Which[    SameQ[answer,0], answer=sign,
  632.                 SameQ[answer, ReImFail], Break[],
  633.                 SameQ[sign,0], _,
  634.                 SameQ[answer,sign], _,
  635.                 True, answer=ReImFail; Break[] ];
  636.             i=i+1];
  637.         answer]
  638.  
  639. SignRe[n_]:=
  640.   (* Exact sign of real part of number. *)
  641.   Module[{sign},
  642.     sign=Sign[Re[n]];
  643.     sign=Which[sign<0, -1,
  644.                 sign==0, 0,
  645.         True, 1];
  646.     sign]
  647.  
  648. SignIm[n_]:=
  649.   (* Exact sign of imaginary part of number. *)
  650.   Module[{sign},
  651.     sign=Sign[Im[n]];
  652.     sign=Which[sign<0, -1,
  653.                 sign==0, 0,
  654.         True, 1];
  655.     sign]
  656.  
  657. (* ****************************************************************
  658. *
  659. *     Re and Im 
  660. *
  661. ***************************************************************** *)
  662.  
  663. ReImExpr[t:f_Symbol[_], True] := {t, 0} /; MemberQ[RealFunctions,f]
  664.  
  665. ReImExpr[Literal[Conjugate[g_]], True] := {1, -1} {Re[g],Im[g]} /.
  666.     ArcTanTemp -> ArcTan
  667.  
  668. ReImExpr[Literal[Abs[g_]], False] := {AbsExpr[g, False],0}
  669. ReImExpr[Literal[Arg[g_]], False] := {ArgExpr[g, False],0}
  670.  
  671. ReImExpr[Re[g_], False] := {(ReImExpr[g,False])[[1]],0}
  672. ReImExpr[Im[g_], False] := {(ReImExpr[g,False])[[2]],0}
  673.  
  674. ReImExpr[Literal[Sign[g_?RealValQ]], False] := {SignExpr[g, False],0}(*try*)
  675. ReImExpr[Literal[Sign[g_]], False] := 
  676.     ReImRecursive[SignExpr[g, False], False] /; IsOK[Sign, g]
  677.  
  678. ReImExpr[Literal[Sign[g_]], False] := Module[{x, y},
  679.     {x, y} = ReImExpr[g, False];
  680.     {x, y} / Sqrt[x^2 + y^2]        
  681.      ] /; ! IsOK[Sign, g] 
  682.  
  683. ReImExpr[Literal[Conjugate[g_]], False] := 
  684.     {1, -1} ReImExpr[g, False];
  685.  
  686. ReImExpr[e_, True] := (ReImRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
  687.     IsRecursive[Re] || IsRecursive[Im] || ConstantQ[e]
  688. ReImExpr[e_, toplevel_:False] := ReImRecursive[e, toplevel] /;
  689.     IsRecursive[Re] || IsRecursive[Im] || ConstantQ[e]
  690.     
  691. ReImExpr[e_, True] := LeafReIm[e, True] /. ArcTanTemp -> ArcTan
  692. ReImExpr[e_, False] := LeafReIm[e, False]
  693.  
  694.  
  695. ReImRecursive[t:f_?IsConjugate[g_?RealValQ], True] := {t, 0}
  696. ReImRecursive[t:f_?IsConjugate [g_?(ConstantQ&&RealValQ)], False] := 
  697.     {t, 0} (*/; MemberQ[ConjugateFunctions, f] && ConstantQ[g]*)
  698.  
  699. ReImRecursive[ArcCos[g_], tl_] :=
  700.     ReImExpr[ Pi/2+I*Log[I*g+Sqrt[1-g^2]] ,tl]
  701. ReImRecursive[ArcCosh[g_?(ZeroQ[# + 1]&)], tl_] := {0, Pi}
  702. ReImRecursive[ArcCosh[g_?(NotZeroQ[# + 1]&)], tl_] := 
  703.     ReImExpr[Log[g+(g+1)*Sqrt[(g-1)/(g+1)]], tl]
  704. ReImRecursive[ArcCot[g_?ZeroQ], tl_] := {Pi/2, 0}
  705. ReImRecursive[ArcCot[g_?NotZeroQ], tl_] :=
  706.     ReImExpr[ -I*Log[(1+I/g)/Sqrt[1+g^(-2)]], tl]
  707. ReImRecursive[ArcCoth[g_?ZeroQ], tl_] := {0, Pi/2}
  708. ReImRecursive[ArcCoth[g_?NotZeroQ], tl_] :=
  709.     ReImExpr[ Log[(1+1/g)/Sqrt[1-g^(-2)]], tl]
  710. ReImRecursive[ArcCsc[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
  711. ReImRecursive[ArcCsc[g_?NotZeroQ], tl_] :=
  712.     ReImExpr[ -I*Log[I/g+Sqrt[1-g^(-2)]], tl]
  713. ReImRecursive[ArcCsch[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
  714. ReImRecursive[ArcCsch[g_?NotZeroQ], tl_] :=
  715.     ReImExpr[ Log[1/g+Sqrt[1+g^(-2)]], tl]
  716. ReImRecursive[ArcSec[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
  717. ReImRecursive[ArcSec[g_?NotZeroQ], tl_] :=
  718.     ReImExpr[ ArcCos[1/g], tl]
  719. ReImRecursive[ArcSech[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
  720. ReImRecursive[ArcSech[g_?NotZeroQ], tl_] :=
  721.     ReImExpr[Log[1/g+(1/g+1)*Sqrt[(1-g)/(1+g)]], tl]
  722. ReImRecursive[ArcSin[g_], tl_] :=
  723.     ReImExpr[ -I*Log[I*g+Sqrt[1-g^2]], tl]
  724. ReImRecursive[ArcSinh[g_], tl_] := ReImExpr[ Log[g+Sqrt[1+g^2]], tl]
  725. ReImRecursive[ArcTan[g_?(ZeroQ[#^2 + 1]&)], tl_] := 
  726.     {Indeterminate, Indeterminate}
  727. ReImRecursive[ArcTan[g_?(NotZeroQ[#^2 + 1]&)], tl_] :=
  728.     ReImExpr[ -I*Log[(1+I*g)/Sqrt[1+g^2]], tl]
  729. ReImRecursive[ArcTan[g_,h_], tl_] := Module[{norm}, norm=GetNorm[g,h];
  730.     ReImExpr[-I*Log[(g+I*h)/(Which[norm===1,1,
  731.          True, Sqrt[norm]])],tl]] /;  
  732.         !TrueQ[RealValQ[g]] || !TrueQ[RealValQ[h]]    
  733.  
  734. ReImRecursive[ArcTanh[g_?(ZeroQ[#^2 - 1]&)], tl_] :=
  735.     {Indeterminate, Indeterminate}
  736. ReImRecursive[ArcTanh[g_?(NotZeroQ[#^2 - 1]&)], tl_] :=
  737.     ReImExpr[ Log[(1+g)/Sqrt[1-g^2]], tl]
  738.  
  739. ReImRecursive[Cos[g_], tl_] :=
  740.   Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Cos[x]*Cosh[y],-Sin[x]*Sinh[y]} ]
  741. ReImRecursive[Cosh[g_], tl_] :=
  742.   Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Cosh[x]*Cos[y],Sinh[x]*Sin[y]} ]
  743. ReImRecursive[Cot[g_], tl_] :=
  744.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  745.       {-Sin[2*x], Sinh[2*y]}/(Cos[2*x]-Cosh[2*y]) ]
  746. ReImRecursive[Coth[g_], tl_] :=
  747.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  748.       {-Sinh[2*x], Sin[2*y]}/(Cos[2*y]-Cosh[2*x]) ]
  749. ReImRecursive[Csc[g_], tl_] :=
  750.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  751.       {-2*Sin[x]*Cosh[y], 2*Cos[x]*Sinh[y]}/(Cos[2*x]-Cosh[2*y]) ]
  752. ReImRecursive[Csch[g_], tl_] :=
  753.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  754.       {-2*Cos[y]*Sinh[x], 2*Sin[y]*Cosh[x]}/(Cos[2*y]-Cosh[2*x]) ]
  755. ReImRecursive[Sec[g_], tl_] :=
  756.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  757.       {2*Cos[x]*Cosh[y], 2*Sin[x]*Sinh[y]}/(Cos[2*x]+Cosh[2*y]) ]
  758. ReImRecursive[Sech[g_], tl_] :=
  759.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  760.       {2*Cos[y]*Cosh[x], -2*Sin[y]*Sinh[x]}/(Cos[2*y]+Cosh[2*x]) ]
  761. ReImRecursive[Sin[g_], tl_] :=
  762.   Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Sin[x]*Cosh[y],Cos[x]*Sinh[y]} ]
  763. ReImRecursive[Sinh[g_], tl_] :=
  764.   Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Sinh[x]*Cos[y],Cosh[x]*Sin[y]} ]
  765. ReImRecursive[Tan[g_], tl_] :=
  766.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  767.       {Sin[2*x], Sinh[2*y]}/(Cos[2*x]+Cosh[2*y]) ]
  768. ReImRecursive[Tanh[g_], tl_] :=
  769.   Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  770.       {Sinh[2*x], Sin[2*y]}/(Cos[2*y]+Cosh[2*x]) ]
  771.  
  772.  
  773. ReImRecursive[Log[g_?ConstantQ], toplevel_] := ConstantLog[g, toplevel]
  774. ReImRecursive[Log[g_], toplevel_] := 
  775.     {RealLog[AbsExpr[g,toplevel]], ArgExpr[g, toplevel]}
  776.  
  777. ReImRecursive[e_Plus, toplevel_] :=  ReImExpr[#, False]& /@ e
  778.  
  779. ReImRecursive[f_Times, toplevel_] :=
  780.     Module[{u,v,x,y,k},
  781.         {u,v}=ReImExpr[f[[1]], toplevel];
  782.         Do[    ({x,y}=ReImExpr[f[[k]], toplevel];
  783.             {u,v}={u*x-v*y,u*y+v*x}),
  784.            {k,2,Length[f]}];
  785.         {u,v} ]
  786.  
  787. ReImRecursive[g_ ^ n_Integer, toplevel_] :=
  788.   Module[{x,y,answer},
  789.     {x,y}=ReImExpr[g, toplevel];
  790.     answer=Which[SameQ[y,0], {x^n,0},
  791.                   SameQ[x,0], If[EvenQ[n],
  792.                           {((-1)^(n/2))*(y^n),0},
  793.                   {0,((-1)^((n-1)/2))*(y^n)}],
  794.           True, Module[{sign,m,u,v},
  795.                   (* Note: N==0 should not occur. *)
  796.                   (* If N<0, use identity (X+I*Y)^N ==
  797.                  (X-I*Y)^(-N)/(X^2+Y^2)^(-N) *)
  798.               sign=Sign[n];
  799.               m=Abs[n];
  800.                   If[sign<0,y=-y];
  801.                           (* Apply Binomial Theorem for (X+I*Y)^N. *)
  802.                           u=Sum[((-1)^(i/2))*Binomial[m,i]*(x^(m-i))*(y^i),
  803.                                  {i,0,m,2}];
  804.                           v=Sum[((-1)^((i-1)/2))*Binomial[m,i]*(x^(m-i))*(y^i),
  805.                                  {i,1,m,2}];
  806.               If[sign<0,
  807.                  Module[{den},
  808.                    den=AbsExpr[g, toplevel]^(2*m);
  809.                    u=u/den;
  810.                    v=v/den]];
  811.               {u,v}]];
  812.     answer]
  813.  
  814. ReImRecursive[g_ ^ h_, toplevel_] :=
  815.     Module[{abs},
  816.         abs=AbsExpr[g, toplevel];
  817.         If[ZeroQ[abs] && PosValQ[h], {0, 0}, Module[{arg,x,y,theta,power},
  818.             arg=ArgExpr[g, toplevel];
  819.             {x,y}=ReImExpr[h, toplevel];
  820.             theta=y*RealLog[abs]+x*arg;
  821.             power=abs^x*E^(-y*arg);
  822.             {power*Cos[theta],power*Sin[theta]}  ]
  823.         ]
  824.     ]
  825.  
  826. ReImRecursive[f_?NumberQ, _] := {Re[f], Im[f]}
  827.  
  828. (* top level *)
  829.  
  830. ReImRecursive[e_?RealValQ, True] := {e, 0}
  831.  
  832. ReImRecursive[e_, True] := {ReImFail, ReImFail}
  833.  
  834. (* leaves *)
  835.  
  836. ReImRecursive[e_, toplevel_] := LeafReIm[e,toplevel]
  837.  
  838. ReImInternalTimes[{x1_,y1_},{x2_,y2_}]:=
  839.   {x1*x2-y1*y2,x1*y2+x2*y1}
  840.  
  841. ReImInternalDivide[{x1_,y1_},{x2_,y2_}]:=
  842.   Module[{norm},
  843.     norm=x2^2+y2^2;
  844.     {(x1*x2+y1*y2)/norm,(x2*y1-x1*y2)/norm}]
  845.  
  846. (* ****************************************************************
  847. *
  848. *     Constant Functions
  849. *
  850. *     These functions implement Log and inverse trigonometric functions 
  851. *on constant args.  These are slightly clever about paying attention
  852. *to branch cuts.
  853. ***************************************************************** *)
  854.  
  855. ConstantSymbols={Catalan,Degree,E,EulerGamma,GoldenRatio,Pi}
  856.  
  857. ConstantAbs[g_, toplevel_:False]:=
  858.   Module[{x,y,answer,norm},
  859.     {x,y}=ReImExpr[g, toplevel]; 
  860.     answer=Which[SameQ[y,0], Switch[SignRe[N[x]],
  861.                                       1, x,
  862.                       0, 0,
  863.                       _, Foil[-1,x]],
  864.                   SameQ[x,0], Switch[SignRe[N[y]],
  865.                               1, y,
  866.                       0, 0,
  867.                       _, Foil[-1,y]],
  868.                   True, If[(norm=GetNorm[x,y])===1,1,Sqrt[norm]]];
  869.     answer]
  870.  
  871. ConstantSign[g_, toplevel_]:= Module[{x,y},
  872.     {x,y}=ReImExpr[g, toplevel];
  873.     Which[SameQ[y,0], SignRe[N[x]],
  874.                   SameQ[x,0], I*SignRe[N[y]],
  875.                   True, If[(norm=GetNorm[x,y])===1, g, g/Sqrt[norm]] ] ]
  876.  
  877. ConstantArg[g_, toplevel_:False]:= (*ConstantLog[g, toplevel][[2]] ]*)
  878.     Module[{x,y},
  879.     {x,y}=ReImExpr[g, toplevel];
  880.     Which [SameQ[y,0], Switch[SignRe[N[x]],
  881.                                       1, 0,
  882.                       0, 0,
  883.                       _, Pi],
  884.          SameQ[x,0], Switch[SignRe[N[y]],
  885.                               1, Pi/2,
  886.                       0, 0,
  887.                       _, -Pi/2],
  888.          True, MyArcTan[x, y, toplevel] ] ]
  889.    
  890. ConstantLog[g_, toplevel_:False]:= Module[{x,y}, 
  891.     {x,y} = ReImExpr[g, toplevel];
  892.     Which[SameQ[y,0], If[N[x]>0,
  893.                                   {Log[x],0},
  894.                   {Log[-x],Pi}],
  895.                   SameQ[x,0], If[N[y]>0,
  896.                           {Log[y],Pi/2},
  897.                   {Log[-y],-Pi/2}],
  898.                   True, Module[{norm,q,log,arctan,r},
  899.             norm = GetNorm[x,y];
  900.                         log=Which[ norm===1, 0,
  901.                   True, RealLog[norm] ];
  902.             arctan = MyArcTan[x, y, toplevel];
  903.             {1/2*log,arctan}] ] ]
  904.    
  905. ConstantArcCot[g_]:=
  906.   Module[{x,y,u,v,k,answer},
  907.     {x,y}=ReImExpr[g];
  908.     {u,v}=If[SameQ[y,0],
  909.                {ArcCot[g],0},
  910.            Module[{arccot,log},
  911.              arccot=ArcCot[2*x/(-1+x^2+y^2)];
  912.          log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
  913.          (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
  914.          {3/4*Pi-1/2*arccot,-1/4*log}]];
  915.     (* Place on proper branch. *)
  916.     k=Round[Re[N[(ArcCot[g]-u)/(Pi/2)]]];
  917.     u=u+k*Pi/2;
  918.     answer={u,v};
  919.     answer]
  920.  
  921. ConstantArcTan[g_]:=
  922.   Module[{x,y,u,v,k,answer},
  923.     {x,y}=ReImExpr[g];
  924.     {u,v}=If[SameQ[y,0],
  925.                {ArcTan[g],0},
  926.            Module[{arctan,log},
  927.              arctan=ArcTan[2*x/(-1+x^2+y^2)];
  928.          log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
  929.          (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
  930.          {-1/2*arctan,1/4*log}]];
  931.     (* Place on proper branch. *)
  932.     k=Round[Re[N[(ArcTan[g]-u)/(Pi/2)]]];
  933.     u=u+k*Pi/2;
  934.     answer={u,v};
  935.     answer]
  936.  
  937. (* *****************************
  938.  
  939. Utility Functions
  940.  
  941. *******************************)
  942.     
  943. GetNorm[(Cos|-Cos)[(theta_|-theta_)], (Sin|-Sin)[(theta_|-theta_)]]:= 1
  944. GetNorm[(Sin|-Sin)[(theta_|-theta_)], (Cos|-Cos)[(theta_|-theta_)]]:= 1
  945. GetNorm[r_ (Cos|-Cos)[(theta_|-theta_)], r_ (Sin|-Sin)[(theta_|-theta_)]]:= r^2
  946. GetNorm[r_ (Sin|-Sin)[(theta_|-theta_)], r_ (Cos|-Cos)[(theta_|-theta_)]]:= r^2
  947. GetNorm[re_,im_]:= re^2 + im^2
  948.  
  949. MyArcTan[x_, y_, toplevel_] :=
  950.     Which[
  951.     NonNegValQ[x] && ZeroQ[y], 0,
  952.     ConstantQ[x] && ConstantQ[y], MyArcTan2[x, y, toplevel],
  953.     True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]
  954.     ]
  955.     
  956. MyArcTan2[x_, y_, toplevel_] :=
  957.     Which[
  958.         PosValQ[x], Which[y===x, Pi/4, y===-x, -Pi/4,
  959.             y/x === Sqrt[3], Pi/3, x/y === Sqrt[3], Pi/6,
  960.             y/x === -Sqrt[3], -Pi/3, x/y === -Sqrt[3], -Pi/6,
  961.             True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]],
  962.          NegValQ[x], Which[ZeroQ[y], Pi, y===x, -3 Pi/4, y===-x, 3 Pi/4,
  963.             y/x === Sqrt[3], -2 Pi/3, x/y === Sqrt[3], -5 Pi/6,
  964.             y/x === -Sqrt[3], 2 Pi/3, x/y === -Sqrt[3], 5 Pi/6,
  965.             True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]],
  966.         PosValQ[y] && ZeroQ[x], Pi/2,
  967.         NegValQ[y] && ZeroQ[x], -Pi/2,
  968.          True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]] ]
  969.  
  970. CESimplify[w_] := w
  971.     
  972. (* ****************************************************************
  973. *
  974. *     Real Functions
  975. *
  976. ***************************************************************** *)
  977.  
  978. (*RealAbs is currently not used by ComplexExpand*)
  979. RealAbs[x_]:=
  980.   Which[ConstantQ[x], If[Re[N[x]]>=0,x,Foil[-1,x]],
  981.         SameQ[SyntacticSign[x],1], Abs[x],
  982.         True, Abs[Foil[-1,x]]]
  983.  
  984. RealLog[x_?NonNegValQ ^ y_?NumberQ] := y RealLog[x] /;
  985.     Im[y] == 0 (*&& AbsExpr[x, toplevel]==x*)
  986.  
  987. RealLog[x_] := 
  988.   If [ Head[x]===Power && x[[1]]===E,  x[[2]], Log[x]]
  989.     
  990. RealSign[x_]:=
  991.   Which[ConstantQ[x], SignRe[N[x]],
  992.         SameQ[SyntacticSign[x],1], Sign[x],
  993.         True, -Sign[Foil[-1,x]]]
  994.  
  995. SyntacticSign[expr_]:=
  996.   (* If the first character of EXPR that will print is "-" then -1,
  997.      if SameQ[EXPR,0] then 0,
  998.      otherwise 1. *)
  999.   Module[{answer},
  1000.     answer=Switch[Head[expr],
  1001.                    Integer, Sign[expr],
  1002.            Rational, Sign[expr],
  1003.            Real, Sign[expr],
  1004.                    Complex, If[SameQ[Re[expr],0],
  1005.                        Sign[Im[expr]],
  1006.                    Sign[Re[expr]]],
  1007.            Plus, SyntacticSign[expr[[1]]],
  1008.            Times, SyntacticSign[expr[[1]]],
  1009.            _, 1];
  1010.     answer]
  1011.  
  1012. Foil[expr1_,expr2_]:=
  1013.   (* Distributive product. *)
  1014.   Module[{i,j,answer},
  1015.     answer=If[Not[SameQ[Head[expr1],Plus]],
  1016.                If[Not[SameQ[Head[expr2],Plus]],
  1017.               expr1*expr2,
  1018.               Sum[expr1*expr2[[j]],
  1019.               {j,Length[expr2]}]],
  1020.            If[Not[SameQ[Head[expr2],Plus]],
  1021.               Sum[expr1[[i]]*expr2,
  1022.               {i,Length[expr2]}],
  1023.               Sum[Sum[expr1[[i]]*expr2[[j]],
  1024.                   {j,Length[expr2]}],
  1025.               {i,Length[expr1]}]]];
  1026.     answer]
  1027.  
  1028. (* ****************************************************************
  1029. *
  1030. *     Install our definitions of Abs, Arg, Re and Im
  1031. *
  1032. *     Info about the C Routines:
  1033. *(1) Abs -- does numbers.
  1034. *(2) Arg -- does infinities, real and imaginary exact numbers, and floats.
  1035. *(3) Conjugate -- does numbers.
  1036. *(4) Im -- does numbers.
  1037. *(5) Re -- does numbers.
  1038. *(6) Sign -- does infinities, numbers, constant symbols, integer Power's
  1039. *recursively, Plus recursively if all args have same sign, and 
  1040. *Times recursively.
  1041. ***************************************************************** *)
  1042.  
  1043. ConstSymbQ[s_Symbol] := MemberQ[ConstantSymbols, s]
  1044.  
  1045. ConstantQ[s_Symbol] := MemberQ[ConstantSymbols, s]
  1046. ConstantQ[s_Symbol[x_]] := MemberQ[AllowedSix,s] && ConstantQ[x]
  1047. ConstantQ[_?NumberQ] := True
  1048. ConstantQ[_?ZeroQ] := True
  1049. ConstantQ[a_ + b_] := ConstantQ[a] && ConstantQ[b]
  1050. ConstantQ[a_ * b_] := ConstantQ[a] && ConstantQ[b]
  1051. ConstantQ[a_ ^ b_] := ConstantQ[a] && ConstantQ[b]
  1052. ConstantQ[_] := False
  1053.  
  1054. RealValQ[s_Symbol] := MemberQ[ConstantSymbols, s] || 
  1055.     (AtomQ[s] && !IsComplex[s])
  1056.  
  1057. RealValQ[s_Symbol[x_]] := MemberQ[RealFunctions,s] ||
  1058.     (MemberQ[ConjugateFunctions,s] && RealValQ[x])
  1059.  
  1060. RealValQ[Log[_?PosValQ]]:= True
  1061. RealValQ[ArcTan[a_,b_]]:= RealValQ[a] && RealValQ[b]
  1062. RealValQ[ArcTanTemp[_,_]]:= True
  1063. RealValQ[(_Real | _Rational | _Integer)] := True
  1064. RealValQ[x_ - I*Im[x_]] := True
  1065. RealValQ[I*Im[x_] - x_] := True
  1066. RealValQ[I*z_ + Im[z_]] := True
  1067. RealValQ[I*z_ - I*Re[z_]] := True
  1068. RealValQ[I*Re[z_] - I*z_] := True
  1069. RealValQ[a_ - Im[a_]] := True
  1070. RealValQ[-a_ + Im[a_]] := True
  1071. RealValQ[a_ + Literal[Conjugate[a_]]] := True
  1072. RealValQ[a_ + b_] := (RealValQ[a] && RealValQ[b]) (*|| a === -Im[b] ||
  1073.     -a === Im[b] || b===Conjugate[a]*)
  1074. RealValQ[a_ * Literal[Conjugate[a_]]] := True
  1075. RealValQ[a_ * b_] := (RealValQ[a] && RealValQ[b]) (*|| b===Conjugate[a]*) ||
  1076.     (ImagValQ[a] && ImagValQ[b])
  1077. RealValQ[a_?PosValQ ^ b_] := RealValQ[b]
  1078. RealValQ[a_^b_?IntegerQ] := RealValQ[a] || (ImagValQ[a] && EvenQ[b])
  1079. RealValQ[_] := False
  1080.  
  1081. ImagValQ[x_ - Re[x_]] := True
  1082. ImagValQ[Re[x_] - x_] := True
  1083. ImagValQ[a_ - Literal[Conjugate[a_]] ] := True
  1084. ImagValQ[-a_ + Literal[Conjugate[a_]] ] := True
  1085. ImagValQ[a_ + Literal[Conjugate[-a_]] ] := True
  1086. ImagValQ[a_ - Re[a]] := True
  1087. ImagValQ[-a_ + Re[a]] := True
  1088. ImagValQ[a_ + Re[-a]] := True
  1089. ImagValQ[a_ + b_ ] := (ImagValQ[a] && ImagValQ[b]) (*|| b === -Conjugate[a] ||
  1090.     -b === Conjugate[a] || b === -Re[a] || -b === Re[a] ||
  1091.     a === -Re[b] || -a === Re[b]*)
  1092. ImagValQ[a_ * b_] := (ImagValQ[a] && RealValQ[b]) ||
  1093.     (ImagValQ[b]&&RealValQ[a])
  1094. ImagValQ[_] := False
  1095.  
  1096. NonNegValQ[Literal[Conjugate[a_]]] := NonNegValQ[a]
  1097. NonNegValQ[Literal[Arg[_?NonNegValQ]]] := True
  1098. NonNegValQ[Literal[Re[_?ImagValQ]]] := True
  1099. NonNegValQ[Literal[Abs[_?ZeroQ]]] := True
  1100. NonNegValQ[Literal[Im[_?RealValQ]]] := True
  1101. NonNegValQ[Literal[Sign[a_]]] := NonNegValQ[a]
  1102. NonNegValQ[ArcTanTemp[_,_?ZeroQ]]:= True
  1103. NonNegValQ[a_ + b_] := NonNegValQ[a] && NonNegValQ[b]
  1104. NonNegValQ[a_ * Literal[Conjugate[a_]]] := True
  1105. NonNegValQ[a_ * b_] := (NonNegValQ[a] && NonNegValQ[b]) ||
  1106.     (NonPosValQ[a] && NonPosValQ[b]) (*|| b===Conjugate[a]*)
  1107. NonNegValQ[a_ ^ b_] := NonNegValQ[a] && RealValQ[b]
  1108. NonNegValQ[a_?RealValQ ^ b_?EvenQ] := True
  1109. NonNegValQ[Literal[Abs[_]]] := True
  1110. NonNegValQ[Literal[Im[_?RealValQ]]] := True
  1111.  
  1112. NonNegValQ[a_Real | a_Rational | a_Integer] := NonNegative[a]
  1113. NonNegValQ[a_?ConstSymbQ] := True
  1114.  
  1115. NonNegValQ[_] := False
  1116.  
  1117. NonPosValQ[a_Real | a_Rational | a_Integer] := Not[Positive[a]]
  1118. NonPosValQ[Literal[Arg[_?NonNegValQ]]] := True
  1119. NonPosValQ[Literal[Conjugate[a_]]] := NonPosValQ[a]
  1120. NonPosValQ[Literal[Sign[a_]]] := NonPosValQ[a]
  1121. NonPosValQ[Literal[Re[_?ImagValQ]]] := True
  1122. NonPosValQ[Literal[Abs[_?ZeroQ]]] := True
  1123. NonPosValQ[Literal[Im[_?RealValQ]]] := True
  1124. NonPosValQ[a_ + b_] := NonPosValQ[a] && NonPosValQ[b]
  1125. NonPosValQ[a_ * Literal[Conjugate[-a_]]] := True
  1126. NonPosValQ[-a_ * Literal[Conjugate[a_]]] := True
  1127. NonPosValQ[a_ * (-Literal[Conjugate[a_]])] := True
  1128. NonPosValQ[a_ * b_] := (NonPosValQ[a] && NonNegValQ[b]) ||
  1129.     (NonNegValQ[a] && NonPosValQ[b])
  1130. NonPosValQ[_] := False
  1131.  
  1132. ZeroQ[Literal[Arg[_?NonNegValQ]]] := True
  1133. ZeroQ[Literal[Im[_?RealValQ]]] := True
  1134. ZeroQ[Literal[Re[_?ImagValQ]]] := True
  1135. ZeroQ[Literal[(Abs | Conjugate | Sign)[_?ZeroQ]]] := True
  1136. ZeroQ[ArcTanTemp[_?ZeroQ,_?ZeroQ]]:= True
  1137. ZeroQ[x_] := NonNegValQ[x] && NonPosValQ[x]
  1138. ZeroQ[_] := False
  1139. NotZeroQ[x_] := ! ZeroQ[x]
  1140.  
  1141. Epsilon = 10^(-30)
  1142.  
  1143. PosValQ[a_?RealValQ] := NonNegValQ[N[a]-Epsilon]
  1144. NegValQ[a_?RealValQ] := NonPosValQ[N[a]+Epsilon]
  1145. PosValQ[_] := False
  1146. NegValQ[_] := False
  1147.  
  1148. AllowedSix={Re, Im, Abs, Arg, Conjugate, Sign}
  1149. AllowedFunctions=AllowedSix
  1150. RecursiveFunctions={Sign} (* for top-level *)
  1151.  
  1152. AssumedComplex={ _ }  (* for top-level *)
  1153.  
  1154. IsComplex[e_] := Or @@ (MatchQ[e, #]& /@ AssumedComplex)
  1155.  
  1156. PolarQ[x_] := IsOK[{Arg, Abs},x]
  1157. CartQ[x_]  := IsOK[{Re, Im},x]
  1158.  
  1159. IsOK[f_List,x_] := And @@ (IsOK[#, x]& /@ f )
  1160. IsOK[f_,x_] := MemberQ[AllowedFunctions, f]
  1161.  
  1162. IsAnyOK[f_List,x_] := Or @@ (IsOK[#, x]& /@ f )
  1163. IsRecursive[f_] := MemberQ[RecursiveFunctions, f]
  1164. IsConjugate[f_] := MemberQ[ConjugateFunctions, f]
  1165.  
  1166. (* symbols ComplexExpand, TargetFunctions
  1167.    have been created by internal code *)
  1168.  
  1169. Begin["`Private`"]
  1170.  
  1171. Unprotect[ComplexExpand]
  1172.  
  1173. Options[ComplexExpand] = {TargetFunctions -> AllowedSix}
  1174.  
  1175. ComplexExpand::exf = "Value of option TargetFunctions -> `1` does not
  1176. contain any of the allowed functions `2`."
  1177.  
  1178. ComplexExpand::rlst="Warning: ComplexExpand does not support `1`. It only supports the rule TargetFunctions -> <function or list of functions>."
  1179.  
  1180. ComplexExpand[e_List, args___] := ComplexExpand[#, args]& /@ e
  1181.  
  1182. ComplexExpand[expr_, opts___Rule] := ComplexExpand[expr, {}, opts]
  1183.  
  1184. ComplexExpand[expr_, cvars_List, opts___Rule] :=
  1185.     Block[{AssumedComplex, AllowedFunctions, RecursiveFunctions,
  1186.            badrulelist, x, `y, `funs},
  1187.         AssumedComplex = cvars;
  1188.         RecursiveFunctions = AllowedSix;
  1189.         badrulelist=Select[{opts},
  1190.             Function[!MatchQ[#,TargetFunctions->_]]];
  1191.         Do[Message[ComplexExpand::rlst,badrulelist[[i]]],
  1192.             {i, Length[badrulelist]}];
  1193.         funs = TargetFunctions /. {opts} /. Options[ComplexExpand];
  1194.         If[Head[funs] =!= List, funs = List[funs]];
  1195.         AllowedFunctions = Select[funs, MemberQ[AllowedSix, #]&];
  1196.         If[ Length[AllowedFunctions] < 1,
  1197.             Message[ComplexExpand::exf, funs, AllowedSix];
  1198.             Return[expr] ];
  1199.         {x, y} = ReImExpr[expr, False];
  1200.         If[ x === ReImFail || y === ReImFail, expr,
  1201.          Block[{BadFunctions, Badpatterns, pos, return = x + I y},
  1202.             BadFunctions = Select[AllowedSix, 
  1203.             ! MemberQ[AllowedFunctions, #]& ];
  1204.             BadPatterns = Apply[Alternatives,
  1205.             Table[BadFunctions[[j]] ,
  1206.             {j, Length[BadFunctions]}]];
  1207.             pos =  Position[return, BadPatterns[_]];
  1208.         RecursiveFunctions={Sign};(* for top-level *)
  1209.         AllowedFunctions=AllowedSix;
  1210.         return = MapAt[ComplexExpand[#, cvars, opts]&,
  1211.                 return, pos] /. {ArcTanTemp->ArcTan};
  1212.              CESimplify[return]
  1213.                  ] 
  1214.         ]
  1215.     ]
  1216.  
  1217.  
  1218. ComplexExpand[expr_, var_, opts___Rule] := ComplexExpand[expr, {var}, opts]
  1219.  
  1220. ComplexExpand[args___] := $Failed /;
  1221.     Message[General::argct, ComplexExpand, Length[{args}]]
  1222.     
  1223.  
  1224. Protect[ComplexExpand]
  1225.  
  1226. End[]
  1227.  
  1228. On[ General::shdw]
  1229.  
  1230. EndPackage[]
  1231.  
  1232. (* do not leave context on path *)
  1233.  
  1234. $ContextPath = DeleteCases[$Packages, "System`ComplexExpand`"]
  1235.  
  1236. Unprotect[$Packages]
  1237. $Packages = DeleteCases[$Packages, "System`ComplexExpand`"]
  1238. Protect[$Packages]
  1239.  
  1240. Null
  1241.