home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e031 / 3.ddi / MATHZIP2 / STARTUP / COMPLEXE.M < prev    next >
Encoding:
Text File  |  1991-09-23  |  25.1 KB  |  837 lines
  1.  
  2. (* ****************************************************************
  3. *
  4. *     ComplexExpand, Kelly Roach and Roman Maeder, V2.0
  5. *
  6. ***************************************************************** *)
  7.  
  8. BeginPackage["System`ComplexExpand`"]
  9.  
  10. (*
  11. This is not too dangerous since the System`ComplexExpand` context is removed at the end
  12. *)
  13.  
  14. Off[ General::shdw]
  15.  
  16. Unprotect[ ComplexExpand, TargetFunctions, AbsExpr, ArgExpr, ConjugateExpr,
  17.            ReImExpr, SignExpr]
  18.  
  19. Map[ Clear, {ComplexExpand, TargetFunctions, AbsExpr, ArgExpr, ConjugateExpr,
  20.            ReImExpr, SignExpr}]
  21.  
  22.  
  23. (* ****************************************************************
  24. *
  25. *     Rigorous Definitions Of The Elementary Inverse Functions
  26. *
  27. *Definitions:
  28. *
  29. *(1) Log[z_]:=If[And[SameQ[Im[z],0],Re[z]<0],Log[-z]+Pi*I,Log[z]]
  30. *(2) Power[u_,v_]:=E^(v*Log[u])
  31. *(3) Sqrt[z_]:=Power[z,1/2]
  32. *(4) ArcSin[z_]:=-I*Log[I*z+Sqrt[1-z^2]]
  33. *(5) ArcCos[z_]:=Pi/2-ArcSin[z]
  34. *(6) ArcTan[z_]:=-I*Log[(1+I*z)/Sqrt[1+z^2]]
  35. *    ArcTan[u_,v_]:=-I*Log[(u+I*v)/Sqrt[u^2+v^2]]
  36. *(7) ArcCsc[z_]:=ArcSin[1/z]
  37. *(8) ArcSec[z_]:=ArcCos[1/z]
  38. *(9) ArcCot[z_]:=If[SameQ[z,0],Pi/2,ArcTan[1/z]]
  39. *(10) ArcSinh[z_]:=Log[z+Sqrt[1+z^2]]
  40. *(11) ArcCosh[z_]:=Log[z+(z+1)*Sqrt[(z-1)/(z+1)]]
  41. *(12) ArcTanh[z_]:=Log[(1+z)/Sqrt[1-z^2]]
  42. *(13) ArcCsch[z_]:=ArcSinh[1/z]
  43. *(14) ArcSech[z_]:=ArcCosh[1/z]
  44. *(15) ArcCoth[z_]:=If[SameQ[z,0],Pi*I/2,ArcTanh[1/z]]
  45. *
  46. *Branch Cuts:
  47. *     
  48. *(1) Log
  49. *     (-Infinity,0) attaches to quadrant II.
  50. *     Im[Log[z]] == Pi     for z in (-Infinity,0)
  51. *(2) Power  (Considering Power for constant v.)
  52. *     (-Infinity,0) attaches to quadrant II, for u not an integer.
  53. *     Power[u,v] == (-u)^v*(Cos[v]+I*Sin[v])     for v in (-Infinity,0)
  54. *(3) Sqrt
  55. *     (-Infinity,0) attaches to quadrant II.
  56. *     Sign[Im[Sqrt[z]]] > 0     for z in (-Infinity,0)
  57. *(4) ArcSin
  58. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  59. *     Re[ArcSin[z]] == -Pi/2     for z in (-Infinity,-1)
  60. *     Re[ArcSin[z]] == Pi/2      for z in (1,Infinity)
  61. *(5) ArcCos
  62. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  63. *     Re[ArcCos[z]] == Pi     for z in (-Infinity,-1)
  64. *     Re[ArcCos[z]] == 0      for z in (1,Infinity)
  65. *(6) ArcTan
  66. *     (-I*Infinity,-I) attaches to quadrant IV and (I,I*Infinity) to quadrant II.
  67. *     Re[ArcTan[z]] == -Pi/2     for z in (-I*Infinity,-I)
  68. *     Re[ArcTan[z]] == Pi/2      for z in (I,I*Infinity)
  69. *(7) ArcCsc
  70. *     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
  71. *     Re[ArcCsc[z]] == -Pi/2     for z in (-1,0).
  72. *     Re[ArcCsc[z]] == Pi/2      for z in (0,1).
  73. *(8) ArcSec
  74. *     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
  75. *     Re[ArcSec[z]] == Pi     for z in (-1,0).
  76. *     Re[ArcSec[z]] == 0      for z in (0,1).
  77. *(9) ArcCot
  78. *     (-I,0) attaches to quadrant III and (0,I) to quadrant I.
  79. *     Re[ArcCot[z]] == Pi/2     for z in (-I,0).
  80. *     Re[ArcCot[z]] == -Pi/2    for z in (0,I).
  81. *(10) ArcSinh
  82. *     (-I*Infinity,-I) attaches to quadrant I and (I,I*Infinity) to quadrant III.
  83. *     Im[ArcSinh[z]] == -Pi/2   for z in (-I*Infinity,-I)
  84. *     Im[ArcSinh[z]] == Pi/2    for z in (I,I*Infinity)
  85. *(11) ArcCosh
  86. *     (-Infinity,1) attaches to quadrants I and II.
  87. *     Sign[Im[ArcCosh[z]]] > 0     for z in (-Infinity,1).
  88. *(12) ArcTanh
  89. *     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
  90. *     Im[ArcTanh[z]]] == Pi/2      for z in (-Infinity,-1).
  91. *     Im[ArcTanh[z]]] == -Pi/2     for z in (1,Infinity).
  92. *(13) ArcCsch
  93. *     (-I,0) attaches to quadrant IV and (0,I) to quadrant I.
  94. *     Sign[Re[ArcCsch[z]]] > 0     for z in (-I,0)
  95. *     Sign[Re[ArcCsch[z]]] < 0     for z in (0,I)
  96. *(14) ArcSech
  97. *     (-Infinity,0) attaches to quadrant III and (1,Infinity) to quadrant IV.
  98. *     Sign[Im[ArcSech[z]]] > 0     for z in (-Infinity,0)
  99. *     Sign[Im[ArcSech[z]]] > 0     for z in (1,Infinity)
  100. *(15) ArcCoth
  101. *     (-1,0] attaches to quadrant III and (0,1) to quadrant I.
  102. *     Im[ArcCoth[z]] == Pi/2        for z in (-1,0]
  103. *     Im[ArcCoth[z]] == -Pi/2      for z in (0,1)
  104. ***************************************************************** *)
  105.  
  106. (* ****************************************************************
  107. *
  108. *     Useful Trigonometric Identities
  109. *
  110. *    Sinh[X] == -I*Sin[I*X] == -Sinh[-X] == I*Sin[-I*X]
  111. *    Cosh[X] == Cos[I*X] == Cosh[-X] == Cos[-I*X]
  112. *    Tanh[X] == -I*Tan[I*X] == -Tanh[-X] == I*Tan[-I*X]
  113. *    ArcSinh[X] == I*ArcSin[-I*X] == Log[X+Sqrt[1+X^2]]
  114. *    Sin[X] == -I*Sinh[I*X] == -Sin[-X] == I*Sinh[-I*X]
  115. *    Cos[X] == Cosh[I*X] == Cos[-X] == Cosh[-I*X]    
  116. *    Tan[X] == -I*Tanh[I*X] == -Tan[-X] == I*Tanh[-I*X]
  117. *    ArcSin[X] == -I*ArcSinh[I*X] == -I*Log[I*X+Sqrt[1-X^2]]
  118. *
  119. *Note: Be very careful with ArcSinh and ArcSin.  The above equations
  120. *are correct.  Cannot assume X is real and must pay attention to branch
  121. *cuts.
  122. ***************************************************************** *)
  123.  
  124. (* ****************************************************************
  125. *
  126. *     Abs, Arg, Conjugate, Im, Re, Sign Identities
  127. *
  128. *For X!=0:
  129. *
  130. *     Abs[X] == Sqrt[X*Conjugate[X]]
  131. *     Arg[X] == -I*Log[X/Sqrt[X*Conjugate[X]]]
  132. *     Im[X] == (X-Conjugate[X])/2
  133. *     Re[X] == (X+Conjugate[X])/2
  134. *     Sign[X] == X/Sqrt[X*Conjugate[X]]
  135. *
  136. *All other relations follow from these relations.  Examples:
  137. *
  138. *     Conjugate[X]=Abs[X]^2/X
  139. *     Re[X]=X/2+X/(2*Sign[X]^2)
  140. *     Im[X]=X*(1-E^(-2*I*Arg[X]))/2
  141. *
  142. *Etc.
  143. ***************************************************************** *)
  144.  
  145. (* ****************************************************************
  146. *
  147. *     Leaf Functions
  148. *
  149. ***************************************************************** *)
  150.  
  151. LeafAbs[f_?NumberQ, _] := Abs[f]
  152.  
  153. LeafAbs[_, True]:= ReImFail /; IsOK[Abs]
  154. LeafAbs[f_, _]:= Abs[f] /; IsOK[Abs]
  155.  
  156. (* cartesian *)
  157. LeafAbs[f_, _]:= Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im}]
  158.  
  159. (* polar: n/a *)
  160.  
  161. (* other cases *)
  162. LeafAbs[f_, toplevel_]:=
  163.     Switch[First[AllowedFunctions],
  164.         Sign, f/Sign[f],
  165.         Arg, f/E^(I*Arg[f]),
  166.         Conjugate, Sqrt[f*Conjugate[f]],
  167.         Im, Sqrt[2*f*Re[f]-f^2],
  168.         Re, Sqrt[f^2-2*I*f*Im[f]],
  169.         _, ReImFail]
  170.  
  171.  
  172. LeafArg[f_?NumberQ, _] := Arg[f]
  173.  
  174. LeafArg[_, True]:= ReImFail /; IsOK[Arg]
  175. LeafArg[f_, _]:= Arg[f] /; IsOK[Arg]
  176.  
  177. (* cartesian *)
  178. LeafArg[f_, _]:= ArcTan[Re[f],Im[f]] /; IsOK[{Re, Im}]
  179.  
  180. (* polar: n/a *)
  181.  
  182. (* other cases *)
  183. LeafArg[f_, toplevel_]:=
  184.     Switch[First[AllowedFunctions],
  185.         Abs, -I*Log[f/Abs[f]],
  186.         Conjugate, -I*Log[f/Sqrt[f*Conjugate[f]]],
  187.         Im, ArcTan[f-I*Im[f],Im[f]],
  188.         Re, ArcTan[Re[f],-I*(f-Re[f])],
  189.         Sign, -I*Log[Sign[f]],
  190.         _, ReImFail]
  191.  
  192.  
  193. LeafConjugate[f_?NumberQ, _] := Conjugate[f]
  194.  
  195. LeafConjugate[f_, _] := f /; !IsComplex[f]
  196.  
  197. LeafConjugate[_, True]:= ReImFail /; IsOK[Conjugate]
  198. LeafConjugate[f_, _]:= Conjugate[f] /; IsOK[Conjugate]
  199.  
  200. (* cartesian *)
  201. LeafConjugate[f_, _]:= Re[f]-I*Im[f] /; IsOK[{Re, Im}]
  202.  
  203. (* polar *)
  204. LeafConjugate[f_, _]:= Abs[f] Exp[-I Arg[f]] /; IsOK[{Arg, Abs}]
  205.  
  206. (* other cases *)
  207. LeafConjugate[f_, toplevel_]:=
  208.     Switch[First[AllowedFunctions],
  209.         Abs, Abs[f]^2/f,
  210.         Arg, f*E^(-2*I*Arg[f]),
  211.         Im, f-2*I*Im[f],
  212.         Re, 2*Re[f]-f,
  213.         Sign, f*Sign[f]^(-2),
  214.         _, ReImFail]
  215.  
  216.  
  217. LeafReIm[f_?NumberQ, _] := {Re[f], Im[f]}
  218.  
  219. LeafReIm[_, True]:= {ReImFail, ReImFail} /; IsOK[{Re, Im}]
  220. LeafReIm[f_, True]:= {ReImFail, -I*(f-Re[f])} /; IsOK[Re]
  221. LeafReIm[f_, True]:= {f-I*Im[f], ReImFail} /; IsOK[Im]
  222.  
  223. (* cartesian *)
  224. LeafReIm[f_, _]:= {Re[f], Im[f]} /; IsOK[{Re, Im}]
  225.  
  226. (* mixed mode, evaluate only once *)
  227. LeafReIm[f_, _]:= With[{re = Re[f]}, {re, -I*(f-re)}] /; IsOK[Re]
  228. LeafReIm[f_, _]:= With[{im = Im[f]}, {f-I*im, im}]    /; IsOK[Im]
  229.  
  230. (* polar *)
  231. LeafReIm[f_, _]:=
  232.     With[{arg = Arg[f], abs = Abs[f]},
  233.         {abs Cos[arg], abs Sin[arg]}
  234.     ] /; IsOK[{Arg, Abs}]
  235.  
  236. (* other cases *)
  237. LeafReIm[f_, toplevel_]:=
  238.     Switch[First[AllowedFunctions],
  239.         Abs, With[{abs = Abs[f]},
  240.               {f/2+abs^2/(2*f), f/(2*I)-abs^2/(2*I*f)}],
  241.         Arg, With[{arg = Arg[f]},
  242.               {f/2+f*E^(-2*I*Arg[f])/2,
  243.                f/(2*I)-f*E^(-2*I*Arg[f])/(2*I)}],
  244.         Conjugate, With[{con = Conjugate[f]},
  245.               {(f+con)/2, (f-con)/(2*I)}],
  246.         Sign, With[{sign = Sign[f]},
  247.               {f/2+f*sign^(-2)/2,
  248.                f/(2*I)-f*sign^(-2)/(2*I)}],
  249.         _, ReImFail]
  250.  
  251.  
  252. LeafSign[f_?NumberQ, _] := Sign[f]
  253.  
  254. LeafSign[_, True]:= ReImFail /; IsOK[Sign]
  255. LeafSign[f_, _]:= Sign[f] /; IsOK[Sign]
  256.  
  257. (* cartesian *)
  258. LeafSign[f_, _]:= (Re[f] + I Im[f])/Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im}]
  259.  
  260. (* polar *)
  261. LeafSign[f_, _]:= f/Abs[f] /; IsOK[{Arg, Abs}]
  262.  
  263. (* other cases *)
  264. LeafSign[f_, toplevel_]:=
  265.     Switch[First[AllowedFunctions],
  266.         Abs, f/Abs[f],
  267.         Arg, Exp[I Arg[f]],
  268.         Conjugate, f/Sqrt[f*Conjugate[f]],
  269.         Re, f/Sqrt[2*f*Re[f]-f^2],
  270.         Im, f/Sqrt[f^2-2*I*f*Im[f]],
  271.         _, ReImFail]
  272.  
  273.  
  274. (* ****************************************************************
  275. *
  276. *     Arg
  277. *
  278. ***************************************************************** *)
  279.  
  280. ArgExpr[Abs[_], _] := 0
  281.  
  282. ArgExpr[e_, toplevel_:False] := ArgRecursive[e, toplevel] /;
  283.     IsRecursive[Arg] || ConstantQ[e]
  284.  
  285. ArgExpr[e_, toplevel_:False] := LeafArg[e, toplevel]
  286.  
  287. ArgRecursive[e_ ^ r_Rational, _] :=
  288.     r Arg[e] /; r < 1 (* principal root *)
  289. ArgRecursive[e_ ^ (_Rational|_Integer), _] := 0 /; Arg[e] == 0 (* positive *)
  290.  
  291. (* special exponentials *)
  292.  
  293. (* internal code makes -1 < Im[c] <= 1 *)
  294. ArgRecursive[ Exp[c_Complex Pi], _ ] := Im[c] Pi
  295. ArgRecursive[ Exp[e_], _ ] := 0 /; Im[e] === 0 (* positive for real e *)
  296.  
  297. ArgRecursive[e_, toplevel_] := LeafArg[e, toplevel]
  298.  
  299. (* ****************************************************************
  300. *
  301. *     Abs
  302. *
  303. ***************************************************************** *)
  304.  
  305. AbsExpr[t:Abs[_], _] := t
  306. (* AbsExpr[Sign[_], _] := 1 : not true for 0 *)
  307.  
  308. AbsExpr[e_, toplevel_:False] := AbsRecursive[e, toplevel] /;
  309.     IsRecursive[Abs] || ConstantQ[e]
  310.  
  311. AbsExpr[e_, toplevel_:False] := LeafAbs[e, toplevel]
  312.  
  313. AbsRecursive[Literal[Conjugate[e_]], _] := Abs[e]
  314. AbsRecursive[e_Plus, toplevel_] := AbsPlus[e, toplevel]
  315. AbsRecursive[e_Times, toplevel_] := Abs /@ e
  316.  
  317. AbsRecursive[g_ ^ n_Integer, _] := 
  318.     Block[{x, y},
  319.         {x,y} = ReImExpr[g];
  320.         Which[    y === 0, If[EvenQ[n], g^n, Abs[g]^n],
  321.             x === 0, If[Mod[n,4] === 0, g^n, Abs[g]^n],
  322.             True, Abs[g]^n]
  323.            ]
  324.  
  325. AbsRecursive[g_ ^ h_, _] := 
  326.     Block[{x, y},
  327.         {x, y} = ReImExpr[h];
  328.         Abs[g]^x E^(-y*Arg[g])
  329.            ]
  330.  
  331. AbsRecursive[e_, toplevel_] := LeafAbs[e, toplevel]
  332.  
  333. AbsPlus[f_Plus, _] := Abs /@ f /; SameSignsQ[f] =!= ReImFail
  334. AbsPlus[e_, toplevel_] := LeafAbs[e, toplevel]
  335.  
  336.  
  337. (* ****************************************************************
  338. *
  339. *     ConjugateExpr
  340. *
  341. ***************************************************************** *)
  342.  
  343. (* commute with Conjugate -> real valued for real arg *)
  344.  
  345. ConjugateFunctions={
  346.   ArcCot, ArcCoth, ArcCsc, ArcCsch, 
  347.   ArcSec, ArcSech, ArcSinh, ArcTan, 
  348.   Cos, Cosh, Cot, Coth, Csc, Csch,
  349.   Sec, Sech, Sin, Sinh, Tan, Tanh,
  350.   Plus, Times, Sign
  351. }
  352.  
  353. (* are real valued always: *)
  354.  
  355. RealFunctions = { Abs, Arg, Im, Re }
  356.  
  357. ConjugateExpr[Literal[Conjugate[f_]], _] := f
  358. ConjugateExpr[t:f_Symbol[_], _] := t /; MemberQ[RealFunctions,f]
  359.  
  360. ConjugateExpr[e_, toplevel_:False] := ConjugateRecursive[e, toplevel] /;
  361.     IsRecursive[Conjugate] || ConstantQ[e]
  362.  
  363. ConjugateExpr[e_, toplevel_:False] := LeafConjugate[e,toplevel]
  364.  
  365. ConjugateRecursive[t:f_Symbol[__], _] := Conjugate /@ t /;
  366.         MemberQ[ConjugateFunctions,f]
  367.  
  368. (* real powers *)
  369. ConjugateRecursive[e_?Positive ^ (r_Rational|r_Integer), _] := e^r
  370.  
  371. ConjugateRecursive[e_, toplevel_] := LeafConjugate[e, toplevel]
  372.  
  373.  
  374. (* ****************************************************************
  375. *
  376. *     SignExpr
  377. *
  378. ***************************************************************** *)
  379.  
  380. SignExpr[t:Sign[_], _] := t
  381.  
  382. SignExpr[e_, toplevel_:False] := SignRecursive[e, toplevel] /;
  383.     IsRecursive[Sign] || ConstantQ[e]
  384.  
  385. SignExpr[e_, toplevel_:False] := LeafSign[e, toplevel]
  386.  
  387. SignRecursive[e_Times, _] := Sign /@ e
  388.  
  389. SignRecursive[e_Plus, toplevel_] :=
  390.     Block[{answer = SameSignsQ[e]},
  391.         If[ answer =!= ReImFail, answer, LeafSign[e, toplevel] ]
  392.     ]
  393.  
  394.  
  395. SignRecursive[g_ ^ n_Integer, _]:=
  396.     Block[{x,y},
  397.         {x,y}=ReImExpr[g];
  398.         Which[    y === 0, If[EvenQ[n],1,Sign[g]],
  399.             x === 0, Switch[Mod[n,4],
  400.                               0, 1,
  401.                       1, I*Sign[g],
  402.                       2, -1,
  403.                       _, -I*Sign[g]],
  404.             True, Sign[g] ^ n]
  405.     ]
  406.  
  407. (* number roots *)
  408.  
  409. SignRecursive[g_?NumberQ ^ r_Rational, _]:= g^r/Abs[g]^r
  410.  
  411. (* exponentials *)
  412.  
  413. SignRecursive[ Exp[e_], _ ] := Exp[ I Im[e] ]
  414.  
  415. (*
  416. SignRecursive[g_ ^ h_, _]:=
  417.     Block[{x, y},
  418.         {x,y}=ReImExpr[h];
  419.         g^h / (Abs[g]^x*E^(-y*Arg[g]))
  420.     ]
  421. *)
  422.  
  423. SignRecursive[e_, toplevel_] := LeafSign[e, toplevel]
  424.  
  425. SameSignsQ[f_]:=
  426.   (* If all args have same sign (skipping 0), then return sign. *)
  427.     Block[{answer = 0, sign, i = 1, n = Length[f]},
  428.         While[i<=n,
  429.             sign=Sign[f[[i]]];
  430.             Which[    SameQ[answer,0], answer=sign,
  431.                 SameQ[sign,0], _,
  432.                 SameQ[answer,sign], _,
  433.                 True, answer=ReImFail; Break[] ];
  434.         i=i+1];
  435.         answer]
  436.  
  437.  
  438. SignRe[n_]:=
  439.   (* Exact sign of real part of number. *)
  440.   Block[{sign},
  441.     sign=Sign[Re[n]];
  442.     sign=Which[sign<0, -1,
  443.                 sign==0, 0,
  444.         True, 1];
  445.     sign]
  446.  
  447. SignIm[n_]:=
  448.   (* Exact sign of imaginary part of number. *)
  449.   Block[{sign},
  450.     sign=Sign[Im[n]];
  451.     sign=Which[sign<0, -1,
  452.                 sign==0, 0,
  453.         True, 1];
  454.     sign]
  455.  
  456. (* ****************************************************************
  457. *
  458. *     Re and Im 
  459. *
  460. ***************************************************************** *)
  461.  
  462. ReImExpr[t:f_Symbol[_], top_:False] := {t, 0} /; MemberQ[RealFunctions,f]
  463.  
  464. ReImExpr[e_, toplevel_:False] := ReImRecursive[e, toplevel] /;
  465.     IsRecursive[Re] || IsRecursive[Im] || ConstantQ[e]
  466.  
  467. ReImExpr[Literal[Conjugate[f_]], True] := {1, -1} ReImExpr[f, False]
  468.  
  469. ReImExpr[e_, toplevel_:False] := LeafReIm[e,toplevel]
  470.  
  471. ReImRecursive[t:f_Symbol[g_], _] := {t, 0} /;
  472.         MemberQ[ConjugateFunctions, f] && Im[g] == 0
  473.  
  474. ReImRecursive[ArcCos[g_], _] := ReImExpr[ Pi/2+I*Log[I*g+Sqrt[1-g^2]] ]
  475. ReImRecursive[ArcCosh[g_], _] := ReImExpr[ Log[g+(g+1)*Sqrt[(g-1)/(g+1)]] ]
  476. ReImRecursive[ArcCot[g_], _] := ReImExpr[ -I*Log[(1+I/g)/Sqrt[1+g^(-2)]] ]
  477. ReImRecursive[ArcCoth[g_], _] := ReImExpr[ Log[(1+1/g)/Sqrt[1-g^(-2)]] ]
  478. ReImRecursive[ArcCsc[g_], _] := ReImExpr[ -I*Log[I/g+Sqrt[1-g^(-2)]] ]
  479. ReImRecursive[ArcCsch[g_], _] := ReImExpr[ Log[1/g+Sqrt[1+g^(-2)]] ]
  480. ReImRecursive[ArcSec[g_], _] := ReImExpr[ ArcCos[1/g] ]
  481. ReImRecursive[ArcSech[g_], _] := ReImExpr[ Log[1/g+(1/g+1)*Sqrt[(1-g)/(1+g)]] ]
  482. ReImRecursive[ArcSin[g_], _] := ReImExpr[ -I*Log[I*g+Sqrt[1-g^2]] ]
  483. ReImRecursive[ArcSinh[g_], _] := ReImExpr[ Log[g+Sqrt[1+g^2]] ]
  484. ReImRecursive[ArcTan[g_], _] := ReImExpr[ -I*Log[(1+I*g)/Sqrt[1+g^2]] ]
  485. ReImRecursive[ArcTanh[g_], _] := ReImExpr[ Log[(1+g)/Sqrt[1-g^2]] ]
  486.  
  487. ReImRecursive[Cos[g_], _] :=
  488.   Block[ {x, y}, {x, y} = ReImExpr[g]; {Cos[x]*Cosh[y],-Sin[x]*Sinh[y]} ]
  489. ReImRecursive[Cosh[g_], _] :=
  490.   Block[ {x, y}, {x, y} = ReImExpr[g]; {Cosh[x]*Cos[y],Sinh[x]*Sin[y]} ]
  491. ReImRecursive[Cot[g_], _] :=
  492.   Block[ {x, y}, {x, y} = ReImExpr[g];
  493.       {-Sin[2*x], Sinh[2*y]}/(Cos[2*x]-Cosh[2*y]) ]
  494. ReImRecursive[Coth[g_], _] :=
  495.   Block[ {x, y}, {x, y} = ReImExpr[g];
  496.       {-Sinh[2*x], Sin[2*y]}/(Cos[2*y]-Cosh[2*x]) ]
  497. ReImRecursive[Csc[g_], _] :=
  498.   Block[ {x, y}, {x, y} = ReImExpr[g];
  499.       {-2*Sin[x]*Cosh[y], 2*Cos[x]*Sinh[y]}/(Cos[2*x]-Cosh[2*y]) ]
  500. ReImRecursive[Csch[g_], _] :=
  501.   Block[ {x, y}, {x, y} = ReImExpr[g];
  502.       {-2*Cos[y]*Sinh[x], 2*Sin[y]*Cosh[x]}/(Cos[2*y]-Cosh[2*x]) ]
  503. ReImRecursive[Sec[g_], _] :=
  504.   Block[ {x, y}, {x, y} = ReImExpr[g];
  505.       {2*Cos[x]*Cosh[y], 2*Sin[x]*Sinh[y]}/(Cos[2*x]+Cosh[2*y]) ]
  506. ReImRecursive[Sech[g_], _] :=
  507.   Block[ {x, y}, {x, y} = ReImExpr[g];
  508.       {2*Cos[y]*Cosh[x], -2*Sin[y]*Sinh[x]}/(Cos[2*y]+Cosh[2*x]) ]
  509. ReImRecursive[Sin[g_], _] :=
  510.   Block[ {x, y}, {x, y} = ReImExpr[g]; {Sin[x]*Cosh[y],Cos[x]*Sinh[y]} ]
  511. ReImRecursive[Sinh[g_], _] :=
  512.   Block[ {x, y}, {x, y} = ReImExpr[g]; {Sinh[x]*Cos[y],Cosh[x]*Sin[y]} ]
  513. ReImRecursive[Tan[g_], _] :=
  514.   Block[ {x, y}, {x, y} = ReImExpr[g];
  515.       {Sin[2*x], Sinh[2*y]}/(Cos[2*x]+Cosh[2*y]) ]
  516. ReImRecursive[Tanh[g_], _] :=
  517.   Block[ {x, y}, {x, y} = ReImExpr[g];
  518.       {Sinh[2*x], Sin[2*y]}/(Cos[2*y]+Cosh[2*x]) ]
  519.  
  520. ReImRecursive[Literal[Conjugate[g_]], _] := {1, -1} ReImExpr[g]
  521.  
  522. ReImRecursive[Log[g_], _] := ConstantLog[g] /; ConstantQ[g]
  523. ReImRecursive[Log[g_], _] := {RealLog[Abs[g]],Arg[g]}
  524.  
  525. ReImRecursive[e_Plus, _] := ReImExpr /@ e (* thanks to listability of Plus *)
  526.  
  527. ReImRecursive[f_Times, _] :=
  528.     Block[{u,v,x,y,k},
  529.         {u,v}=ReImExpr[f[[1]]];
  530.         Do[    {x,y}=ReImExpr[f[[k]]];
  531.             {u,v}={u*x-v*y,u*y+v*x},
  532.            {k,2,Length[f]}];
  533.         {u,v} ]
  534.  
  535. ReImRecursive[g_ ^ n_Integer, _] :=
  536.   Block[{x,y,answer},
  537.     {x,y}=ReImExpr[g];
  538.     answer=Which[SameQ[y,0], {x^n,0},
  539.                   SameQ[x,0], If[EvenQ[n],
  540.                           {((-1)^(n/2))*(y^n),0},
  541.                   {0,((-1)^((n-1)/2))*(y^n)}],
  542.           True, Block[{sign,m,u,v},
  543.                   (* Note: N==0 should not occur. *)
  544.                   (* If N<0, use identity (X+I*Y)^N ==
  545.                  (X-I*Y)^(-N)/(X^2+Y^2)^(-N) *)
  546.               sign=Sign[n];
  547.               m=Abs[n];
  548.                   If[sign<0,y=-y];
  549.                           (* Apply Binomial Theorem for (X+I*Y)^N. *)
  550.                           u=Sum[((-1)^(i/2))*Binomial[m,i]*(x^(m-i))*(y^i),
  551.                                  {i,0,m,2}];
  552.                           v=Sum[((-1)^((i-1)/2))*Binomial[m,i]*(x^(m-i))*(y^i),
  553.                                  {i,1,m,2}];
  554.               If[sign<0,
  555.                  Block[{den},
  556.                    den=Abs[g]^(2*m);
  557.                    u=u/den;
  558.                    v=v/den]];
  559.               {u,v}]];
  560.     answer]
  561.  
  562. ReImRecursive[g_ ^ h_, _] :=
  563.     Block[{abs,arg,x,y,theta,power},
  564.         abs=Abs[g];
  565.         arg=Arg[g];
  566.         {x,y}=ReImExpr[h];
  567.         theta=y*RealLog[abs]+x*arg;
  568.         power=abs^x*E^(-y*arg);
  569.         {power*Cos[theta],power*Sin[theta]}
  570.     ]
  571.  
  572. ReImRecursive[f_?NumberQ, _] := {Re[f], Im[f]}
  573.  
  574. (* top level *)
  575.  
  576. ReImRecursive[e_, True] := {e, 0} /; !IsComplex[e]
  577. ReImRecursive[ _, True] := {ReImFail, ReImFail}
  578.  
  579. (* leaves *)
  580.  
  581. ReImRecursive[e_, toplevel_] := LeafReIm[e,toplevel]
  582.  
  583.  
  584. ReImInternalTimes[{x1_,y1_},{x2_,y2_}]:=
  585.   {x1*x2-y1*y2,x1*y2+x2*y1}
  586.  
  587. ReImInternalDivide[{x1_,y1_},{x2_,y2_}]:=
  588.   Block[{norm},
  589.     norm=x2^2+y2^2;
  590.     {(x1*x2+y1*y2)/norm,(x1*y2-x2*y1)/norm}]
  591.  
  592. (* ****************************************************************
  593. *
  594. *     Constant Functions
  595. *
  596. *     These functions implement Log and inverse trigonometric functions 
  597. *on constant args.  These are slightly clever about paying attention
  598. *to branch cuts.
  599. ***************************************************************** *)
  600.  
  601. ConstantSymbols={Catalan,Degree,E,EulerGamma,GoldenRatio,Pi}
  602.  
  603. ConstantAbs[g_]:=
  604.   Block[{x,y,answer},
  605.     {x,y}=ReImExpr[g];
  606.     answer=Which[SameQ[y,0], Switch[SignRe[N[x]],
  607.                                       1, x,
  608.                       0, 0,
  609.                       _, Foil[-1,x]],
  610.                   SameQ[x,0], Switch[SignRe[N[y]],
  611.                               1, y,
  612.                       0, 0,
  613.                       _, Foil[-1,y]],
  614.                   True, Sqrt[x^2+y^2]];
  615.     answer]
  616.  
  617. ConstantArg[g_]:=
  618.   Block[{x,y,answer},
  619.     {x,y}=ReImExpr[g];
  620.     answer=Which[SameQ[y,0], Switch[SignRe[N[x]],
  621.                                       1, 0,
  622.                       0, 0,
  623.                       _, Pi],
  624.                   SameQ[x,0], Switch[SignRe[N[y]],
  625.                               1, Pi/2,
  626.                       0, 0,
  627.                       _, -Pi/2],
  628.                   True, ReImInternalTimes[
  629.                   {0,-1},
  630.                   ConstantLog[ReImInternalDivide[{x,y},ConstantAbs[g]]]]];
  631.     answer]
  632.  
  633. ConstantSign[g_]:=
  634.   Block[{x,y,answer},
  635.     {x,y}=ReImExpr[g];
  636.     answer=Which[SameQ[y,0], SignRe[N[x]],
  637.                   SameQ[x,0], I*SignRe[N[y]],
  638.                   True, g/ConstantAbs[g]];
  639.     answer]
  640.  
  641. ConstantArcCot[g_]:=
  642.   Block[{x,y,u,v,k,answer},
  643.     {x,y}=ReImExpr[g];
  644.     {u,v}=If[SameQ[y,0],
  645.                {ArcCot[g],0},
  646.            Block[{arccot,log},
  647.              arccot=ArcCot[2*x/(-1+x^2+y^2)];
  648.          log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
  649.          (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
  650.          {3/4*Pi-1/2*arccot,-1/4*log}]];
  651.     (* Place on proper branch. *)
  652.     k=Round[Re[N[(ArcCot[g]-u)/(Pi/2)]]];
  653.     u=u+k*Pi/2;
  654.     answer={u,v};
  655.     answer]
  656.  
  657. ConstantArcTan[g_]:=
  658.   Block[{x,y,u,v,k,answer},
  659.     {x,y}=ReImExpr[g];
  660.     {u,v}=If[SameQ[y,0],
  661.                {ArcTan[g],0},
  662.            Block[{arctan,log},
  663.              arctan=ArcTan[2*x/(-1+x^2+y^2)];
  664.          log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
  665.          (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
  666.          {-1/2*arctan,1/4*log}]];
  667.     (* Place on proper branch. *)
  668.     k=Round[Re[N[(ArcTan[g]-u)/(Pi/2)]]];
  669.     u=u+k*Pi/2;
  670.     answer={u,v};
  671.     answer]
  672.  
  673. ConstantLog[g_]:=
  674.   Block[{x,y,answer},
  675.     {x,y}=ReImExpr[g];
  676.     answer=Which[SameQ[y,0], If[N[x]>0,
  677.                                   {Log[x],0},
  678.                   {Log[-x],Pi}],
  679.                   SameQ[x,0], If[N[y]>0,
  680.                           {Log[y],Pi/2},
  681.                   {Log[-y],-Pi/2}],
  682.                   True, Block[{norm,q,log,arctan},
  683.                       norm=x^2+y^2;
  684.                       q=y/x;
  685.                           log=Log[norm];
  686.               arctan=ArcTan[q];
  687.               {1/2*log,arctan}]];
  688.     answer]
  689.  
  690. (* ****************************************************************
  691. *
  692. *     Real Functions
  693. *
  694. ***************************************************************** *)
  695.  
  696. RealAbs[x_]:=
  697.   Which[ConstantQ[x], If[Re[N[x]]>=0,x,Foil[-1,x]],
  698.         SameQ[SyntacticSign[x],1], Abs[x],
  699.         True, Abs[Foil[-1,x]]]
  700.  
  701. RealLog[x_ ^ y_?NumberQ] := y Log[x] /; Im[y] == 0
  702. RealLog[x_] := Log[x]
  703.  
  704. RealSign[x_]:=
  705.   Which[ConstantQ[x], SignRe[N[x]],
  706.         SameQ[SyntacticSign[x],1], Sign[x],
  707.         True, -Sign[Foil[-1,x]]]
  708.  
  709. SyntacticSign[expr_]:=
  710.   (* If the first character of EXPR that will print is "-" then -1,
  711.      if SameQ[EXPR,0] then 0,
  712.      otherwise 1. *)
  713.   Block[{answer},
  714.     answer=Switch[Head[expr],
  715.                    Integer, Sign[expr],
  716.            Rational, Sign[expr],
  717.            Real, Sign[expr],
  718.                    Complex, If[SameQ[Re[expr],0],
  719.                        Sign[Im[expr]],
  720.                    Sign[Re[expr]]],
  721.            Plus, SyntacticSign[expr[[1]]],
  722.            Times, SyntacticSign[expr[[1]]],
  723.            _, 1];
  724.     answer]
  725.  
  726. Foil[expr1_,expr2_]:=
  727.   (* Distributive product. *)
  728.   Block[{i,j,answer},
  729.     answer=If[Not[SameQ[Head[expr1],Plus]],
  730.                If[Not[SameQ[Head[expr2],Plus]],
  731.               expr1*expr2,
  732.               Sum[expr1*expr2[[j]],
  733.               {j,Length[expr2]}]],
  734.            If[Not[SameQ[Head[expr2],Plus]],
  735.               Sum[expr1[[i]]*expr2,
  736.               {i,Length[expr2]}],
  737.               Sum[Sum[expr1[[i]]*expr2[[j]],
  738.                   {j,Length[expr2]}],
  739.               {i,Length[expr1]}]]];
  740.     answer]
  741.  
  742. (* ****************************************************************
  743. *
  744. *     Install our definitions of Abs, Arg, Re and Im
  745. *
  746. *     Info about the C Routines:
  747. *(1) Abs -- does numbers.
  748. *(2) Arg -- does infinities, real and imaginary exact numbers, and floats.
  749. *(3) Conjugate -- does numbers.
  750. *(4) Im -- does numbers.
  751. *(5) Re -- does numbers.
  752. *(6) Sign -- does infinities, numbers, constant symbols, integer Power's
  753. *recursively, Plus recursively if all args have same sign, and 
  754. *Times recursively.
  755. ***************************************************************** *)
  756.  
  757. ConstantQ[s_Symbol] := MemberQ[ConstantSymbols, s]
  758. ConstantQ[_?NumberQ] := True
  759. ConstantQ[a_ + b_] := ConstantQ[a] && ConstantQ[b]
  760. ConstantQ[a_ * b_] := ConstantQ[a] && ConstantQ[b]
  761. ConstantQ[a_ ^ b_] := ConstantQ[a] && ConstantQ[b]
  762. ConstantQ[_] = False
  763.  
  764. AllowedSix={Re, Im, Abs, Arg, Conjugate, Sign}
  765. AllowedFunctions=AllowedSix
  766. RecursiveFunctions={Sign} (* for top-level *)
  767.  
  768. AssumedComplex={ _ }  (* for top-level *)
  769.  
  770. IsComplex[e_] := Or @@ (MatchQ[e, #]& /@ AssumedComplex)
  771.  
  772. IsOK[f_List] := And @@ IsOK /@ f
  773. IsOK[f_] := MemberQ[AllowedFunctions, f]
  774.  
  775. IsRecursive[f_] := MemberQ[RecursiveFunctions, f]
  776.  
  777. PolarQ := IsOK[{Arg, Abs}]
  778. CartQ  := IsOK[{Re, Im}]
  779.  
  780. (* symbols ComplexExpand, TargetFunctions
  781.    have been created by internal code *)
  782.  
  783. Begin["`Private`"]
  784.  
  785. Unprotect[ComplexExpand]
  786.  
  787. Options[ComplexExpand] = {TargetFunctions -> AllowedSix}
  788.  
  789. ComplexExpand::exf = "Value of option TargetFunctions -> `1` does not
  790. contain any of the allowed functions `2`."
  791.  
  792. ComplexExpand[e_List, args___] := ComplexExpand[#, args]& /@ e
  793.  
  794. ComplexExpand[expr_, opts___Rule] := ComplexExpand[expr, {}, opts]
  795.  
  796. ComplexExpand[expr_, cvars_List, opts___Rule] :=
  797.     Block[{AssumedComplex, AllowedFunctions, RecursiveFunctions,
  798.            `x, `y, `funs},
  799.         AssumedComplex = cvars;
  800.         RecursiveFunctions = AllowedSix;
  801.         funs = TargetFunctions /. {opts} /. Options[ComplexExpand];
  802.         If[Head[funs] =!= List, funs = List[funs]];
  803.         AllowedFunctions = Select[funs, MemberQ[AllowedSix, #]&];
  804.         If[ Length[AllowedFunctions] < 1,
  805.             Message[ComplexExpand::exf, funs, AllowedSix];
  806.             Return[expr] ];
  807.         Update /@ AllowedSix;
  808.         {x, y} = ReImExpr[expr];
  809.         If[ x === ReImFail || y === ReImFail, expr, x + I y]
  810.     ]
  811.  
  812. ComplexExpand[expr_, var_, opts___Rule] := ComplexExpand[expr, {var}, opts]
  813.  
  814. ComplexExpand[args___] := $Failed /;
  815.     Message[General::argct, ComplexExpand, Length[{args}]]
  816.  
  817. Protect[ComplexExpand]
  818.  
  819. End[]
  820.  
  821. On[ General::shdw]
  822.  
  823. EndPackage[]
  824.  
  825. (* do not leave context on path *)
  826.  
  827. $ContextPath = Drop[ $ContextPath,
  828.     Position[$ContextPath, "System`ComplexExpand`"][[1]] ]
  829.  
  830. Unprotect[$Packages]
  831. $Packages = Drop[$Packages,
  832.     Position[$Packages, "System`ComplexExpand`"][[1]]]
  833. Protect[$Packages]
  834.  
  835. Null
  836.  
  837.