
(* ****************************************************************
*
*     ComplexExpand, Roman Maeder, V2.0
*     Debugged/enhanced for V2.1 by Daniel Lichtblau
***************************************************************** *)

BeginPackage["System`ComplexExpand`"]

(*
This is not too dangerous since the System`ComplexExpand` context is removed at
the end
*)

Off[ General::shdw]

Unprotect[ {ComplexExpand, TargetFunctions, System`ComplexExpand`AbsExpr,
            System`ComplexExpand`ArgExpr, System`ComplexExpand`ConjugateExpr,
            System`ComplexExpand`ReImExpr, System`ComplexExpand`SignExpr}]


(********************************************************************** *)

(* ****************************************************************
*
*     Rigorous Definitions Of The Elementary Inverse Functions
*
*Definitions:
*
*(1) Log[z_]:=If[And[SameQ[Im[z],0],Re[z]<0],Log[-z]+Pi*I,Log[z]]
*(2) Power[u_,v_]:=E^(v*Log[u])
*(3) Sqrt[z_]:=Power[z,1/2]
*(4) ArcSin[z_]:=-I*Log[I*z+Sqrt[1-z^2]]
*(5) ArcCos[z_]:=Pi/2-ArcSin[z]
*(6) ArcTan[z_]:=-I*Log[(1+I*z)/Sqrt[1+z^2]]
*    ArcTan[u_,v_]:=-I*Log[(u+I*v)/Sqrt[u^2+v^2]]
*(7) ArcCsc[z_]:=ArcSin[1/z]
*(8) ArcSec[z_]:=ArcCos[1/z]
*(9) ArcCot[z_]:=If[SameQ[z,0],Pi/2,ArcTan[1/z]]
*(10) ArcSinh[z_]:=Log[z+Sqrt[1+z^2]]
*(11) ArcCosh[z_]:=Log[z+(z+1)*Sqrt[(z-1)/(z+1)]]
*(12) ArcTanh[z_]:=Log[(1+z)/Sqrt[1-z^2]]
*(13) ArcCsch[z_]:=ArcSinh[1/z]
*(14) ArcSech[z_]:=ArcCosh[1/z]
*(15) ArcCoth[z_]:=If[SameQ[z,0],Pi*I/2,ArcTanh[1/z]]
*
*Branch Cuts:
*     
*(1) Log
*     (-Infinity,0) attaches to quadrant II.
*     Im[Log[z]] == Pi     for z in (-Infinity,0)
*(2) Power  (Considering Power for constant v.)
*     (-Infinity,0) attaches to quadrant II, for u not an integer.
*     Power[u,v] == (-u)^v*(Cos[v]+I*Sin[v])     for v in (-Infinity,0)
*(3) Sqrt
*     (-Infinity,0) attaches to quadrant II.
*     Sign[Im[Sqrt[z]]] > 0     for z in (-Infinity,0)
*(4) ArcSin
*     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
*     Re[ArcSin[z]] == -Pi/2     for z in (-Infinity,-1)
*     Re[ArcSin[z]] == Pi/2      for z in (1,Infinity)
*(5) ArcCos
*     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
*     Re[ArcCos[z]] == Pi     for z in (-Infinity,-1)
*     Re[ArcCos[z]] == 0      for z in (1,Infinity)
*(6) ArcTan
*     (-I*Infinity,-I) attaches to quadrant IV and (I,I*Infinity) to quadrant II.
*     Re[ArcTan[z]] == -Pi/2     for z in (-I*Infinity,-I)
*     Re[ArcTan[z]] == Pi/2      for z in (I,I*Infinity)
*(7) ArcCsc
*     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
*     Re[ArcCsc[z]] == -Pi/2     for z in (-1,0).
*     Re[ArcCsc[z]] == Pi/2      for z in (0,1).
*(8) ArcSec
*     (-1,0) attaches to quadrant III and (0,1) to quadrant I.
*     Re[ArcSec[z]] == Pi     for z in (-1,0).
*     Re[ArcSec[z]] == 0      for z in (0,1).
*(9) ArcCot
*     (-I,0) attaches to quadrant III and (0,I) to quadrant I.
*     Re[ArcCot[z]] == Pi/2     for z in (-I,0).
*     Re[ArcCot[z]] == -Pi/2    for z in (0,I).
*(10) ArcSinh
*     (-I*Infinity,-I) attaches to quadrant I and (I,I*Infinity) to quadrant III.
*     Im[ArcSinh[z]] == -Pi/2   for z in (-I*Infinity,-I)
*     Im[ArcSinh[z]] == Pi/2    for z in (I,I*Infinity)
*(11) ArcCosh
*     (-Infinity,1) attaches to quadrants I and II.
*     Sign[Im[ArcCosh[z]]] > 0     for z in (-Infinity,1).
*(12) ArcTanh
*     (-Infinity,-1) attaches to quadrant II and (1,Infinity) to quadrant IV.
*     Im[ArcTanh[z]]] == Pi/2      for z in (-Infinity,-1).
*     Im[ArcTanh[z]]] == -Pi/2     for z in (1,Infinity).
*(13) ArcCsch
*     (-I,0) attaches to quadrant IV and (0,I) to quadrant I.
*     Sign[Re[ArcCsch[z]]] > 0     for z in (-I,0)
*     Sign[Re[ArcCsch[z]]] < 0     for z in (0,I)
*(14) ArcSech
*     (-Infinity,0) attaches to quadrant III and (1,Infinity) to quadrant IV.
*     Sign[Im[ArcSech[z]]] > 0     for z in (-Infinity,0)
*     Sign[Im[ArcSech[z]]] > 0     for z in (1,Infinity)
*(15) ArcCoth
*     (-1,0] attaches to quadrant III and (0,1) to quadrant I.
*     Im[ArcCoth[z]] == Pi/2        for z in (-1,0]
*     Im[ArcCoth[z]] == -Pi/2      for z in (0,1)
***************************************************************** *)

(* ****************************************************************
*
*     Useful Trigonometric Identities
*
*    Sinh[X] == -I*Sin[I*X] == -Sinh[-X] == I*Sin[-I*X]
*    Cosh[X] == Cos[I*X] == Cosh[-X] == Cos[-I*X]
*    Tanh[X] == -I*Tan[I*X] == -Tanh[-X] == I*Tan[-I*X]
*    ArcSinh[X] == I*ArcSin[-I*X] == Log[X+Sqrt[1+X^2]]
*    Sin[X] == -I*Sinh[I*X] == -Sin[-X] == I*Sinh[-I*X]
*    Cos[X] == Cosh[I*X] == Cos[-X] == Cosh[-I*X]	
*    Tan[X] == -I*Tanh[I*X] == -Tan[-X] == I*Tanh[-I*X]
*    ArcSin[X] == -I*ArcSinh[I*X] == -I*Log[I*X+Sqrt[1-X^2]]
*
*Note: Be very careful with ArcSinh and ArcSin.  The above equations
*are correct.  Cannot assume X is real and must pay attention to branch
*cuts.
***************************************************************** *)

(* ****************************************************************
*
*     Abs, Arg, Conjugate, Im, Re, Sign Identities
*
*For X!=0:
*
*     Abs[X] == Sqrt[X*Conjugate[X]]
*     Arg[X] == -I*Log[X/Sqrt[X*Conjugate[X]]]
*     Im[X] == (X-Conjugate[X])/2
*     Re[X] == (X+Conjugate[X])/2
*     Sign[X] == X/Sqrt[X*Conjugate[X]]
*
*All other relations follow from these relations.  Examples:
*
*     Conjugate[X]=Abs[X]^2/X
*     Re[X]=X/2+X/(2*Sign[X]^2)
*     Im[X]=X*(1-E^(-2*I*Arg[X]))/2
*
*Etc.
***************************************************************** *)

(* ****************************************************************
*
*     Leaf Functions
*
***************************************************************** *)

LeafAbs[f_?NumberQ, True] := Abs[f] /; IsOK[Abs,f]
LeafAbs[f_?ConstantQ, toplevel_] := ConstantAbs[f, toplevel]

LeafAbs[f_, True] := ReImFail /; IsOK[Abs,f]
(* LeafAbs[f_, True] :=Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im},f] *)

LeafAbs[f_?RealValQ, False] := Sqrt[f^2]

LeafAbs[f_, False] := Abs[f] /; IsOK[Abs,f]

(* cartesian *)
LeafAbs[f_, False] := Module[{x,y,norm},
	{x,y}=ReImExpr[f,False];
	norm=GetNorm[x,y];
	Which[norm===1,1,
		True, Sqrt[x^2 + y^2]]
	] /; IsOK[{Re, Im},f]
	
LeafAbs[f_, True]:= Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im},f]

(* polar: n/a *)

(* other cases *)
LeafAbs[f_, False]:=
	Which[  IsOK[Arg, f], f/E^(I*ArgExpr[f, False]),
		IsOK[Conjugate, f], Sqrt[f*ConjugateExpr[f, False]],
		IsOK[Im, f], Sqrt[f^2-2*I*f*ReImExpr[f,False][[2]]],
		IsOK[Re, f], Sqrt[2*f*ReImExpr[f, False][[1]]-f^2],
		IsOK[Sign, f], f/SignExpr[f, False],
		True, ReImFail]

LeafAbs[f_, True]:=
	Which[	IsOK[Arg, f], f/E^(I*Arg[f]),
		IsOK[Conjugate, f], Sqrt[f*Conjugate[f]],
		IsOK[Im, f], Sqrt[2*f*Re[f]-f^2],
		IsOK[Re, f], Sqrt[f^2-2*I*f*Im[f]],
		IsOK[Sign, f], f/Sign[f],
		True, ReImFail]


LeafArg[f_?NumberQ, _] := Arg[f]

LeafArg[f_?ConstantQ, toplevel_] := ConstantArg[f, toplevel]

LeafArg[f_, True]:= ReImFail /; IsOK[Arg,f]

LeafArg[f_, False]:= Arg[f] /; IsOK[Arg,f] 

(* cartesian *)
LeafArg[f_, True]:= ArcTan[Re[f],Im[f]] /; IsOK[{Re, Im},f]
LeafArg[f_, False]:= 
  Module[{x,y,q},
	{x,y}=ReImExpr[f, False];
	(*q=y/x;
	If [Head[q]===Tan && N[-Pi < N[q[[1]]] <= Pi], q[[1]],]*) 	
	MyArcTan[x, y, False]
	   ]  /; IsOK[{Re, Im},f] 

(* polar: n/a *)

(* other cases *)
LeafArg[f_, False] := 
	Which[  IsOK[Im, f], With[{im=ReImExpr[f, False][[2]]},
			MyArcTan[f-I*im, im, False]],
		IsOK[Re, f], With[{re=ReImExpr[f, False][[1]]},
			MyArcTan[re, -I*(f-re), False]],
		IsOK[Sign, f], -I*Log[SignExpr[f, False]],
		IsOK[Abs, f], -I*Log[f/AbsExpr[f, False]],
		IsOK[Conjugate, f], -I*Log[f/Sqrt[f*ConjugateExpr[f, False]]],
		True, ReImFail]

LeafArg[f_, True] :=
	Which[  IsOK[Im, f], ArcTan[f-I*Im[f],Im[f]],
		IsOK[Re, f], ArcTan[Re[f],-I*(f-Re[f])],
		IsOK[Sign, f], -I*Log[Sign[f]],
		IsOK[Abs, f], -I*Log[f/Abs[f]],
		IsOK[Conjugate, f], -I*Log[f/Sqrt[f*Conjugate[f]]],
		True, ReImFail]


LeafConjugate[f_?RealValQ, False] := f

LeafConjugate[f_?NumberQ, _] := Conjugate[f]
LeafConjugate[g_?ConstantQ, True] := Re[g] - I Im[g]

LeafConjugate[f_, True]:= ReImFail /; IsOK[Conjugate,f]
LeafConjugate[f_, False]:= Conjugate[f] /; IsOK[Conjugate,f]

(* cartesian *)
LeafConjugate[f_, True]:= Re[f]-I*Im[f] /; IsOK[{Re, Im}]
LeafConjugate[f_, False]:= Module[{x,y},
	{x,y} = ReImExpr[f, False];
	x-I*y ]  /; IsOK[{Re, Im},f] 

(* polar *)
LeafConjugate[f_, True]:= Abs[f] Exp[-I Arg[f]] /; IsOK[{Arg, Abs}]
LeafConjugate[f_, False]:= AbsExpr[f, False] * 
	Exp[-I ArgExpr[f, False]] /; IsOK[{Arg, Abs},f]

(* other cases *)
LeafConjugate[f_, False]:=
	Which[  IsOK[Arg, f], f*E^(-2*I*ArgExpr[f, False]),
		IsOK[Im, f], f-2*I*ReImExpr[f, False][[2]],
		IsOK[Re, f], 2*ReImExpr[f, False][[1]]-f,
		IsOK[Abs, f], AbsExpr[f, False]^2/f,
		IsOK[Sign, f], f*SignExpr[f, False]^(-2),
		True, ReImFail]

LeafConjugate[f_, True]:=
	Which[  IsOK[Arg, f], f*E^(-2*I*Arg[f]),
		IsOK[Im, f], f-2*I*Im[f],
		IsOK[Re, f], 2*Re[f]-f,
		IsOK[Abs, f], Abs[f]^2/f,
		IsOK[Sign, f], f*Sign[f]^(-2),
		True, ReImFail]


LeafReIm[f_?NumberQ, _] := {Re[f], Im[f]}
LeafReIm[f_?ConstSymbQ, _] := {Re[f], Im[f]}

LeafReIm[f_, True]:= {ReImFail, ReImFail} /; IsOK[{Re, Im},f]
LeafReIm[f_, True]:= {ReImFail, -I*(f-Re[f])} /; IsOK[Re,f]
LeafReIm[f_, True]:= {f-I*Im[f], ReImFail} /; IsOK[Im,f]

(* cartesian *)

LeafReIm[f_?RealValQ, False]:= {f, 0}
LeafReIm[f_?((I * RealValQ) || (RealValQ * I)), False]:= {0, f}

LeafReIm[f_, False]:= {Re[f], Im[f]} /; IsOK[{Re, Im},f]

(* mixed mode, evaluate only once *)
LeafReIm[f_, False]:= With[{re = Re[f]}, {re, -I*(f-re)}] /;
IsOK[Re,f]
LeafReIm[f_, False]:= With[{im = Im[f]}, {f-I*im, im}] /;
IsOK[Im,f]

(* polar *)
LeafReIm[f_, False]:=
	With[{arg = ArgExpr[f, False], abs = AbsExpr[f, False]},
		{abs Cos[arg], abs Sin[arg]}
	] /; IsOK[{Arg, Abs},f]
LeafReIm[f_, True]:=
	With[{arg = Arg[f], abs = Abs[f]},
		{abs Cos[arg], abs Sin[arg]}
	] /; IsOK[{Arg, Abs},f]

(* other cases *)
LeafReIm[f_, False]:=
	Which[ 	IsOK[Conjugate, f], With[{conj = ConjugateExpr[f, False]},
		      {(f+conj)/2, (f-conj)/(2*I)}],
		IsOK[Arg, f], With[{arg = ArgExpr[f, False]},
		      {f/2+f*E^(-2*I*arg)/2,
		       f/(2*I)-f*E^(-2*I*arg)/(2*I)}],
		IsOK[Abs, f], With[{abs = AbsExpr[f, False]},
		      {f/2+abs^2/(2*f), f/(2*I)-abs^2/(2*I*f)}],
		IsOK[Sign, f], With[{sign = SignExpr[f, False]},
		      {f/2+f*sign^(-2)/2,
		       f/(2*I)-f*sign^(-2)/(2*I)}],
		True, ReImFail]

LeafReIm[f_, True]:=
	Which[  IsOK[Conjugate, f], With[{con = Conjugate[f]},
		      {(f+con)/2, (f-con)/(2*I)}],
		IsOK[Arg, f], With[{arg = Arg[f]},
		      {f/2+f*E^(-2*I*Arg[f])/2,
		       f/(2*I)-f*E^(-2*I*Arg[f])/(2*I)}],
		IsOK[Abs, f], With[{abs = Abs[f]},
		      {f/2+abs^2/(2*f), f/(2*I)-abs^2/(2*I*f)}],
		IsOK[Sign, f], With[{sign = Sign[f]},
		      {f/2+f*sign^(-2)/2,
		       f/(2*I)-f*sign^(-2)/(2*I)}],
		True, ReImFail]


LeafSign[f_?NumberQ, True] := Sign[f]
(*LeafSign[f_?ConstSymbQ, True] := Sign[f]*)

LeafSign[f_?ConstantQ, toplevel_] := ConstantSign[f, toplevel]

LeafSign[f_, True]:= ReImFail /; IsOK[Sign,f]
LeafSign[f_, False]:= Sign[f] /; IsOK[Sign,f]

(* cartesian *)
LeafSign[f_, True]:= (Re[f] + I Im[f])/Sqrt[Re[f]^2+Im[f]^2] /; IsOK[{Re, Im}]
LeafSign[f_, False]:= Module[{x,y,norm,normsqrt}, 
			{x,y}=ReImExpr[f, False];
			norm=GetNorm[x,y];
			If[norm===1,normsqrt=1,normsqrt=Sqrt[norm]];
			(x + I y)/normsqrt
		] /; IsOK[{Re, Im},f]

(* other cases *)
LeafSign[f_, False]:=
	Which[  IsOK[Arg, f], Exp[I ArgExpr[f, False]],
		IsOK[Abs, f], f/AbsExpr[f, False],
		IsOK[Conjugate, f], 
			f/Sqrt[f*ConjugateExpr[f, False]],
		IsOK[Re, f], f/Sqrt[2*f*ReImExpr[f, False][[1]]-f^2],
		IsOK[Im, f], f/Sqrt[f^2-2*I*f*ReImExpr[f, False][[2]]],
		True, ReImFail] 

LeafSign[f_, True]:=
	Which[  IsOK[Arg, f], Exp[I Arg[f]],
		IsOK[Abs, f], f/Abs[f],
		IsOK[Conjugate, f], f/Sqrt[f*Conjugate[f]],
		IsOK[Re, f], f/Sqrt[2*f*Re[f]-f^2],
		IsOK[Im, f], f/Sqrt[f^2-2*I*f*Im[f]],
		True, ReImFail]


(* End of leaf code ********************************************** *)

(* ****************************************************************
*
*     Arg
*
***************************************************************** *)

ArgExpr[x_?NonNegValQ, False] := 0
ArgExpr[Abs[_], _] := 0
ArgExpr[Arg[x_?RealValQ], False] := 0
ArgExpr[Arg[Arg[_]], False] := 0
ArgExpr[Sign[x_], False] := ArgExpr[x, False]

ArgExpr[Re[x_], False] := With[{re=ReImExpr[x, False][[1]]},
	Which[re===Re[x], ArgRecursive[re, False],
	True, ArgExpr[re, False] ] ]

ArgExpr[Im[x_], False] := With[{im=ReImExpr[x, False][[2]]},
	Which[im===Im[x], ArgRecursive[im, False],
	True, ArgExpr[im, False] ] ]

ArgExpr[Arg[x_], False] := With[{arg=ArgExpr[x, False]},
	Which[arg===Arg[x], ArgRecursive[arg, False],
	True, ArgExpr[arg, False] ] ]

ArgExpr[Conjugate[x_], False] := With[{conj=ConjugateExpr[x, False]},
	Which[conj===Conjugate[x], ArgRecursive[conj, False],
	True, ArgExpr[conj, False] ] ]

ArgExpr[f_?ConstantQ, False] := ConstantArg[f, False]

ArgExpr[e_, True] := (ArgRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
	IsRecursive[Arg] || ConstantQ[e]
ArgExpr[e_, False] := ArgRecursive[e, False] /;
	IsRecursive[Arg] || ConstantQ[e]

ArgExpr[e_, True] := LeafArg[e, True] /. ArcTanTemp -> ArcTan
ArgExpr[e_, False] := LeafArg[e, False]

ArgRecursive[e_ ^ r_Rational, toplevel_] :=
	r ArgExpr[e,toplevel] /; r < 1 (* principal root *)

ArgRecursive[e_ ^ (_Rational|_Integer), False] := 0 /;
	 ArgExpr[e, toplevel] == 0
ArgRecursive[e_ ^ (_Rational|_Integer), True] := 0 /;
	 Arg[e] == 0

ArgRecursive[e_?RealValQ ^ p_Integer, False] := 0 /; EvenQ[p]

(* internal code makes -1 < Im[c] <= 1 *)
ArgRecursive[ Exp[c_Complex Pi], _ ] := Im[c] Pi
ArgRecursive[ Exp[e_], _ ] := 0 /; Im[e] === 0 (* Exp[e] positive for real e *)

ArgRecursive[e_ ^ p_, toplevel_:False] := 
	Module[{abs,arg,x,y,theta},
		abs=AbsExpr[e,toplevel];
		arg=ArgExpr[e,toplevel];
		{x,y}=ReImExpr[p,toplevel];
		theta = y * RealLog[abs] + x * arg;
		MyArcTan[Cos[theta], Sin[theta], toplevel]
	]

(* special exponentials *)

ArgRecursive[e_, toplevel_] := LeafArg[e, toplevel]

(* ****************************************************************
*
*     Abs
*
***************************************************************** *)

AbsExpr[t:Abs[_], True] := t

AbsExpr[Abs[x_], False] := AbsExpr[x,False]

AbsExpr[Re[x_], False] := With[{re=ReImExpr[x, False][[1]]},
	Which[re===Re[x], AbsRecursive[re, False],
	True, AbsExpr[re, False] ] ]

AbsExpr[Im[x_], False] := With[{im=ReImExpr[x, False][[2]]},
	Which[im===Im[x], AbsRecursive[im, False],
	True, AbsExpr[im, False] ] ]

AbsExpr[Sign[e_], False] := Which[
	ZeroQ[e], 0,
	PosValQ[e] || NegValQ[e], 1,
	True, With[{sign=SignExpr[e, False]},
		If[sign===Sign[e],AbsRecursive[sign, False],
			AbsExpr[sign, False]] ] ]	

AbsExpr[x_?NonNegValQ, False] := x

AbsExpr[Conjugate[e_], toplevel_] := AbsExpr[e, toplevel]
AbsExpr[Literal[Arg[e_]], False] := Sqrt[ArgExpr[e, False]^2]

AbsExpr[e_, True] := (AbsRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
	IsRecursive[Abs] || ConstantQ[e]
AbsExpr[e_, False] := AbsRecursive[e, False] /;
	IsRecursive[Abs] || ConstantQ[e]

AbsExpr[e_, True] := LeafAbs[e, True] /. ArcTanTemp -> ArcTan
AbsExpr[e_, False] := LeafAbs[e, False]

AbsRecursive[e_Plus, toplevel_] := AbsPlus[e, toplevel]
AbsRecursive[e_Times, toplevel_] :=  AbsExpr[#, False]& /@ e

AbsRecursive[g_ ^ n_Integer, toplevel_] := 
	Module[{x, y},
		{x,y} = ReImExpr[g, toplevel];
		Which[	y === 0, If[EvenQ[n], g^n, AbsExpr[g, toplevel]^n],
			x === 0, If[Mod[n,4] === 0, g^n,
			    AbsExpr[g, toplevel]^n],
			True, AbsExpr[g, toplevel]^n]
           ]

AbsRecursive[t:g_ ^ h_, toplevel_] := 
	Module[{x, y}, 
		{x, y} = ReImExpr[h, toplevel];
		AbsExpr[g, toplevel]^x E^(-y*ArgExpr[g, toplevel])
           ]

AbsRecursive[e_, toplevel_] := LeafAbs[e, toplevel]

AbsPlus[f_Plus, toplevel_] := AbsExpr[#, False]& /@ f /;
	SameSignsQ[f, toplevel] =!= ReImFail

AbsPlus[e_, toplevel_] := LeafAbs[e, toplevel]

(* ****************************************************************
*
*     ConjugateExpr
*
***************************************************************** *)

(* commute with Conjugate -> real valued for real arg *)

ConjugateFunctions={
  ArcCot,  ArcCsch, ArcSinh, ArcTan, 
  Cos, Cosh, Cot, Coth, Csc, Csch,
  Sec, Sech, Sin, Sinh, Tan, Tanh,
  Plus, Times, Sign(* Sign --Handle this a bit differently in ReImExpr *)
}

(* are real valued always: *)

RealFunctions = { Abs, Arg, Im, Re }

ConjugateExpr[z_?(!ConstantQ), toplevel_] := Module[{x,y}, 
	{x,y} =  ReImExpr[z, toplevel]; 
	Which[ (x===ReImFail) || (y===ReImFail), ReImFail,
	True, x - I y] ] /; PolarQ[z] || CartQ[z] || ! IsOK[Conjugate, z]

ConjugateExpr[Conjugate[f_], _] := f

ConjugateExpr[t:f_Symbol[_], True] := t /; MemberQ[RealFunctions,f]
ConjugateExpr[t:f_Symbol[_], False] := ReImExpr[t, False][[1]] /;
	MemberQ[RealFunctions,f]

ConjugateExpr[(-1) ^ y_?RealValQ,_] := (-1)^(-y)
ConjugateExpr[x_?NegValQ ^ y_?RealValQ,_] := (-1)^(-y) * (-x)^y
ConjugateExpr[x_?PosValQ ^ Complex[0,z_], _] := x^(-Complex[0,z])
ConjugateExpr[x_?PosValQ ^ (y_?RealValQ * Complex[0,z_]), _] :=
	x^(-y*Complex[0,z])

ConjugateExpr[e_, True] :=
	(ConjugateRecursive[e, True] /.ArcTanTemp -> ArcTan) /;
		IsRecursive[Conjugate] || ConstantQ[e]
ConjugateExpr[e_, False] := ConjugateRecursive[e, False] /;
	IsRecursive[Conjugate] || ConstantQ[e]
	
ConjugateExpr[e_, True] := LeafConjugate[e, True] /. ArcTanTemp -> ArcTan
ConjugateExpr[e_, False] := LeafConjugate[e, False]

ConjugateRecursive[t:f_Symbol[__], True] := Conjugate /@ t /; 
		MemberQ[ConjugateFunctions,f] && ConstantQ[t]

ConjugateRecursive[t:f_Symbol[__],toplevel_] := 
	ConjugateExpr[#, False]& /@ t /;
		MemberQ[ConjugateFunctions,f] && f =!= Sign && !ConstantQ[t]

(* real powers *)
ConjugateRecursive[e_?PosValQ ^ (r_Rational|r_Integer), _] := e^r

ConjugateRecursive[e_?RealValQ, toplevel_] := e

ConjugateRecursive[e_, toplevel_] := LeafConjugate[e, toplevel]


(* ****************************************************************
*
*     SignExpr
*
***************************************************************** *)

SignExpr[t:Sign[_], True] := t
SignExpr[Sign[x_], False] := SignExpr[x,False]

SignExpr[t:f_Symbol[_], False] := With[{x=ReImExpr[t,
False][[1]]},
	Which[ZeroQ[x], 0,
	PosValQ[x], 1,
	NegValQ[x], -1,
	x===t, SignRecursive[x, False],
	True, SignExpr[x, False]
	 ] ] /; MemberQ[RealFunctions,f]

SignExpr[e_, True] := (SignRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
	IsRecursive[Sign] || ConstantQ[e]

SignExpr[e_, False] := SignRecursive[e, False] /;
	IsRecursive[Sign] || ConstantQ[e]

SignExpr[e_, True] := LeafSign[e, True] /. ArcTanTemp -> ArcTan
SignExpr[e_, False] := LeafSign[e, False]

SignRecursive[e_Times, False] := SignExpr[#, False]& /@ e
SignRecursive[e_Times, True] := Sign /@ e

SignRecursive[e_Plus, toplevel_] :=
	Module[{answer = SameSignsQ[e, toplevel]},
		If[ answer =!= ReImFail, answer, LeafSign[e, toplevel] ]
	] 

SignRecursive[g_ ^ n_Integer, True]:=
	Module[{x,y},
		{x,y}=ReImExpr[g, False];
		Which[	y === 0, If[EvenQ[n], Sign[g]^2, Sign[g]],
			x === 0, Switch[Mod[n,4],
		                      0, Sign[g]^2,
				      1, I*Sign[g],
				      2, -(Sign[g]^2),
				      _, -I*Sign[g]],
			True, Sign[g]^n]
	]
	
SignRecursive[g_ ^ n_Integer, False]:=
	Module[{x,y},
		{x,y}=ReImExpr[g,False];
		Which[	y === 0, If [EvenQ[n],
			SignExpr[g, False]^2, SignExpr[g, False]],
			x === 0, Switch[Mod[n,4],
		                      0, SignExpr[g, False]^2,
				      1, I*SignExpr[g, False],
				      2, -(SignExpr[g, False]^2),
				      _, -I*SignExpr[g, False]], 
			True, SignExpr[g, False] ^ n]
	]

(* number roots *)

SignRecursive[g_ ^ r_?RealValQ, toplevel_] :=
	Which[toplevel, Sign[g]^r,
	True, SignExpr[g, toplevel] ^ r]
	
(* exponentials *)

SignRecursive[ Exp[e_], _ ] := Exp[ I Im[e] ]

SignRecursive[e_, toplevel_] := LeafSign[e, toplevel]

SameSignsQ[f_, toplevel_]:=
  (* If all args have same sign (skipping 0), then return sign. *)
	Module[{answer = 0, sign, i = 1, n = Length[f]},
		While[i<=n,
			(*sign=SignExpr[f[[i]], toplevel];*)
			sign=Sign[f[[i]]];
			Which[	SameQ[answer,0], answer=sign,
				SameQ[answer, ReImFail], Break[],
				SameQ[sign,0], _,
				SameQ[answer,sign], _,
				True, answer=ReImFail; Break[] ];
			i=i+1];
		answer]

SignRe[n_]:=
  (* Exact sign of real part of number. *)
  Module[{sign},
    sign=Sign[Re[n]];
    sign=Which[sign<0, -1,
                sign==0, 0,
		True, 1];
    sign]

SignIm[n_]:=
  (* Exact sign of imaginary part of number. *)
  Module[{sign},
    sign=Sign[Im[n]];
    sign=Which[sign<0, -1,
                sign==0, 0,
		True, 1];
    sign]

(* ****************************************************************
*
*     Re and Im 
*
***************************************************************** *)

ReImExpr[t:f_Symbol[_], True] := {t, 0} /; MemberQ[RealFunctions,f]

ReImExpr[Literal[Conjugate[g_]], True] := {1, -1} {Re[g],Im[g]} /.
	ArcTanTemp -> ArcTan

ReImExpr[Literal[Abs[g_]], False] := {AbsExpr[g, False],0}
ReImExpr[Literal[Arg[g_]], False] := {ArgExpr[g, False],0}

ReImExpr[Re[g_], False] := {(ReImExpr[g,False])[[1]],0}
ReImExpr[Im[g_], False] := {(ReImExpr[g,False])[[2]],0}

ReImExpr[Literal[Sign[g_?RealValQ]], False] := {SignExpr[g, False],0}(*try*)
ReImExpr[Literal[Sign[g_]], False] := 
	ReImRecursive[SignExpr[g, False], False] /; IsOK[Sign, g]

ReImExpr[Literal[Sign[g_]], False] := Module[{x, y},
	{x, y} = ReImExpr[g, False];
	{x, y} / Sqrt[x^2 + y^2]		
	 ] /; ! IsOK[Sign, g] 

ReImExpr[Literal[Conjugate[g_]], False] := 
	{1, -1} ReImExpr[g, False];

ReImExpr[e_, True] := (ReImRecursive[e, True] /. ArcTanTemp -> ArcTan) /;
	IsRecursive[Re] || IsRecursive[Im] || ConstantQ[e]
ReImExpr[e_, toplevel_:False] := ReImRecursive[e, toplevel] /;
	IsRecursive[Re] || IsRecursive[Im] || ConstantQ[e]
	
ReImExpr[e_, True] := LeafReIm[e, True] /. ArcTanTemp -> ArcTan
ReImExpr[e_, False] := LeafReIm[e, False]


ReImRecursive[t:f_?IsConjugate[g_?RealValQ], True] := {t, 0}
ReImRecursive[t:f_?IsConjugate [g_?(ConstantQ&&RealValQ)], False] := 
	{t, 0} (*/; MemberQ[ConjugateFunctions, f] && ConstantQ[g]*)

ReImRecursive[ArcCos[g_], tl_] :=
	ReImExpr[ Pi/2+I*Log[I*g+Sqrt[1-g^2]] ,tl]
ReImRecursive[ArcCosh[g_?(ZeroQ[# + 1]&)], tl_] := {0, Pi}
ReImRecursive[ArcCosh[g_?(NotZeroQ[# + 1]&)], tl_] := 
	ReImExpr[Log[g+(g+1)*Sqrt[(g-1)/(g+1)]], tl]
ReImRecursive[ArcCot[g_?ZeroQ], tl_] := {Pi/2, 0}
ReImRecursive[ArcCot[g_?NotZeroQ], tl_] :=
	ReImExpr[ -I*Log[(1+I/g)/Sqrt[1+g^(-2)]], tl]
ReImRecursive[ArcCoth[g_?ZeroQ], tl_] := {0, Pi/2}
ReImRecursive[ArcCoth[g_?NotZeroQ], tl_] :=
	ReImExpr[ Log[(1+1/g)/Sqrt[1-g^(-2)]], tl]
ReImRecursive[ArcCsc[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
ReImRecursive[ArcCsc[g_?NotZeroQ], tl_] :=
	ReImExpr[ -I*Log[I/g+Sqrt[1-g^(-2)]], tl]
ReImRecursive[ArcCsch[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
ReImRecursive[ArcCsch[g_?NotZeroQ], tl_] :=
	ReImExpr[ Log[1/g+Sqrt[1+g^(-2)]], tl]
ReImRecursive[ArcSec[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
ReImRecursive[ArcSec[g_?NotZeroQ], tl_] :=
	ReImExpr[ ArcCos[1/g], tl]
ReImRecursive[ArcSech[g_?ZeroQ], tl_] := {Indeterminate, Indeterminate}
ReImRecursive[ArcSech[g_?NotZeroQ], tl_] :=
	ReImExpr[Log[1/g+(1/g+1)*Sqrt[(1-g)/(1+g)]], tl]
ReImRecursive[ArcSin[g_], tl_] :=
	ReImExpr[ -I*Log[I*g+Sqrt[1-g^2]], tl]
ReImRecursive[ArcSinh[g_], tl_] := ReImExpr[ Log[g+Sqrt[1+g^2]], tl]
ReImRecursive[ArcTan[g_?(ZeroQ[#^2 + 1]&)], tl_] := 
	{Indeterminate, Indeterminate}
ReImRecursive[ArcTan[g_?(NotZeroQ[#^2 + 1]&)], tl_] :=
	ReImExpr[ -I*Log[(1+I*g)/Sqrt[1+g^2]], tl]
ReImRecursive[ArcTan[g_,h_], tl_] := Module[{norm}, norm=GetNorm[g,h];
	ReImExpr[-I*Log[(g+I*h)/(Which[norm===1,1,
		 True, Sqrt[norm]])],tl]] /;  
		!TrueQ[RealValQ[g]] || !TrueQ[RealValQ[h]]	

ReImRecursive[ArcTanh[g_?(ZeroQ[#^2 - 1]&)], tl_] :=
	{Indeterminate, Indeterminate}
ReImRecursive[ArcTanh[g_?(NotZeroQ[#^2 - 1]&)], tl_] :=
	ReImExpr[ Log[(1+g)/Sqrt[1-g^2]], tl]

ReImRecursive[Cos[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Cos[x]*Cosh[y],-Sin[x]*Sinh[y]} ]
ReImRecursive[Cosh[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Cosh[x]*Cos[y],Sinh[x]*Sin[y]} ]
ReImRecursive[Cot[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{-Sin[2*x], Sinh[2*y]}/(Cos[2*x]-Cosh[2*y]) ]
ReImRecursive[Coth[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{-Sinh[2*x], Sin[2*y]}/(Cos[2*y]-Cosh[2*x]) ]
ReImRecursive[Csc[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{-2*Sin[x]*Cosh[y], 2*Cos[x]*Sinh[y]}/(Cos[2*x]-Cosh[2*y]) ]
ReImRecursive[Csch[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{-2*Cos[y]*Sinh[x], 2*Sin[y]*Cosh[x]}/(Cos[2*y]-Cosh[2*x]) ]
ReImRecursive[Sec[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{2*Cos[x]*Cosh[y], 2*Sin[x]*Sinh[y]}/(Cos[2*x]+Cosh[2*y]) ]
ReImRecursive[Sech[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{2*Cos[y]*Cosh[x], -2*Sin[y]*Sinh[x]}/(Cos[2*y]+Cosh[2*x]) ]
ReImRecursive[Sin[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Sin[x]*Cosh[y],Cos[x]*Sinh[y]} ]
ReImRecursive[Sinh[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl]; {Sinh[x]*Cos[y],Cosh[x]*Sin[y]} ]
ReImRecursive[Tan[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{Sin[2*x], Sinh[2*y]}/(Cos[2*x]+Cosh[2*y]) ]
ReImRecursive[Tanh[g_], tl_] :=
  Module[ {x, y}, {x, y} = ReImExpr[g, tl];
  	{Sinh[2*x], Sin[2*y]}/(Cos[2*y]+Cosh[2*x]) ]


ReImRecursive[Log[g_?ConstantQ], toplevel_] := ConstantLog[g, toplevel]
ReImRecursive[Log[g_], toplevel_] := 
	{RealLog[AbsExpr[g,toplevel]], ArgExpr[g, toplevel]}

ReImRecursive[e_Plus, toplevel_] :=  ReImExpr[#, False]& /@ e

ReImRecursive[f_Times, toplevel_] :=
	Module[{u,v,x,y,k},
		{u,v}=ReImExpr[f[[1]], toplevel];
		Do[	({x,y}=ReImExpr[f[[k]], toplevel];
			{u,v}={u*x-v*y,u*y+v*x}),
		   {k,2,Length[f]}];
		{u,v} ]

ReImRecursive[g_ ^ n_Integer, toplevel_] :=
  Module[{x,y,answer},
    {x,y}=ReImExpr[g, toplevel];
    answer=Which[SameQ[y,0], {x^n,0},
                  SameQ[x,0], If[EvenQ[n],
		                  {((-1)^(n/2))*(y^n),0},
				  {0,((-1)^((n-1)/2))*(y^n)}],
		  True, Module[{sign,m,u,v},
		          (* Note: N==0 should not occur. *)
		          (* If N<0, use identity (X+I*Y)^N ==
			     (X-I*Y)^(-N)/(X^2+Y^2)^(-N) *)
			  sign=Sign[n];
			  m=Abs[n];
		          If[sign<0,y=-y];
                          (* Apply Binomial Theorem for (X+I*Y)^N. *)
                          u=Sum[((-1)^(i/2))*Binomial[m,i]*(x^(m-i))*(y^i),
                                 {i,0,m,2}];
                          v=Sum[((-1)^((i-1)/2))*Binomial[m,i]*(x^(m-i))*(y^i),
                                 {i,1,m,2}];
			  If[sign<0,
			     Module[{den},
			       den=AbsExpr[g, toplevel]^(2*m);
			       u=u/den;
			       v=v/den]];
			  {u,v}]];
    answer]

ReImRecursive[g_ ^ h_, toplevel_] :=
	Module[{abs},
		abs=AbsExpr[g, toplevel];
		If[ZeroQ[abs] && PosValQ[h], {0, 0}, Module[{arg,x,y,theta,power},
			arg=ArgExpr[g, toplevel];
			{x,y}=ReImExpr[h, toplevel];
			theta=y*RealLog[abs]+x*arg;
			power=abs^x*E^(-y*arg);
			{power*Cos[theta],power*Sin[theta]}  ]
		]
	]

ReImRecursive[f_?NumberQ, _] := {Re[f], Im[f]}

(* top level *)

ReImRecursive[e_?RealValQ, True] := {e, 0}

ReImRecursive[e_, True] := {ReImFail, ReImFail}

(* leaves *)

ReImRecursive[e_, toplevel_] := LeafReIm[e,toplevel]

ReImInternalTimes[{x1_,y1_},{x2_,y2_}]:=
  {x1*x2-y1*y2,x1*y2+x2*y1}

ReImInternalDivide[{x1_,y1_},{x2_,y2_}]:=
  Module[{norm},
    norm=x2^2+y2^2;
    {(x1*x2+y1*y2)/norm,(x2*y1-x1*y2)/norm}]

(* ****************************************************************
*
*     Constant Functions
*
*     These functions implement Log and inverse trigonometric functions 
*on constant args.  These are slightly clever about paying attention
*to branch cuts.
***************************************************************** *)

ConstantSymbols={Catalan,Degree,E,EulerGamma,GoldenRatio,Pi}

ConstantAbs[g_, toplevel_:False]:=
  Module[{x,y,answer,norm},
    {x,y}=ReImExpr[g, toplevel]; 
    answer=Which[SameQ[y,0], Switch[SignRe[N[x]],
                                      1, x,
				      0, 0,
				      _, Foil[-1,x]],
                  SameQ[x,0], Switch[SignRe[N[y]],
		                      1, y,
				      0, 0,
				      _, Foil[-1,y]],
                  True, If[(norm=GetNorm[x,y])===1,1,Sqrt[norm]]];
    answer]

ConstantSign[g_, toplevel_]:= Module[{x,y},
    {x,y}=ReImExpr[g, toplevel];
    Which[SameQ[y,0], SignRe[N[x]],
                  SameQ[x,0], I*SignRe[N[y]],
                  True, If[(norm=GetNorm[x,y])===1, g, g/Sqrt[norm]] ] ]

ConstantArg[g_, toplevel_:False]:= (*ConstantLog[g, toplevel][[2]] ]*)
	Module[{x,y},
	{x,y}=ReImExpr[g, toplevel];
	Which [SameQ[y,0], Switch[SignRe[N[x]],
                                      1, 0,
				      0, 0,
				      _, Pi],
         SameQ[x,0], Switch[SignRe[N[y]],
		                      1, Pi/2,
				      0, 0,
				      _, -Pi/2],
         True, MyArcTan[x, y, toplevel] ] ]
   
ConstantLog[g_, toplevel_:False]:= Module[{x,y}, 
	{x,y} = ReImExpr[g, toplevel];
	Which[SameQ[y,0], If[N[x]>0,
                                  {Log[x],0},
				  {Log[-x],Pi}],
                  SameQ[x,0], If[N[y]>0,
		                  {Log[y],Pi/2},
				  {Log[-y],-Pi/2}],
                  True, Module[{norm,q,log,arctan,r},
			norm = GetNorm[x,y];
                        log=Which[ norm===1, 0,
			  	True, RealLog[norm] ];
			arctan = MyArcTan[x, y, toplevel];
			{1/2*log,arctan}] ] ]
   
ConstantArcCot[g_]:=
  Module[{x,y,u,v,k,answer},
    {x,y}=ReImExpr[g];
    {u,v}=If[SameQ[y,0],
               {ArcCot[g],0},
	       Module[{arccot,log},
	         arccot=ArcCot[2*x/(-1+x^2+y^2)];
		 log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
		 (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
		 {3/4*Pi-1/2*arccot,-1/4*log}]];
    (* Place on proper branch. *)
    k=Round[Re[N[(ArcCot[g]-u)/(Pi/2)]]];
    u=u+k*Pi/2;
    answer={u,v};
    answer]

ConstantArcTan[g_]:=
  Module[{x,y,u,v,k,answer},
    {x,y}=ReImExpr[g];
    {u,v}=If[SameQ[y,0],
               {ArcTan[g],0},
	       Module[{arctan,log},
	         arctan=ArcTan[2*x/(-1+x^2+y^2)];
		 log=Log[(x^2+(y+1)^2)/(x^2+(y-1)^2)];
		 (* log=Log[(1+x^2+2*y+y^2)/(1+x^2-2*y+y^2)]; *)
		 {-1/2*arctan,1/4*log}]];
    (* Place on proper branch. *)
    k=Round[Re[N[(ArcTan[g]-u)/(Pi/2)]]];
    u=u+k*Pi/2;
    answer={u,v};
    answer]

(* *****************************

Utility Functions

*******************************)
    
GetNorm[(Cos|-Cos)[(theta_|-theta_)], (Sin|-Sin)[(theta_|-theta_)]]:= 1
GetNorm[(Sin|-Sin)[(theta_|-theta_)], (Cos|-Cos)[(theta_|-theta_)]]:= 1
GetNorm[r_ (Cos|-Cos)[(theta_|-theta_)], r_ (Sin|-Sin)[(theta_|-theta_)]]:= r^2
GetNorm[r_ (Sin|-Sin)[(theta_|-theta_)], r_ (Cos|-Cos)[(theta_|-theta_)]]:= r^2
GetNorm[re_,im_]:= re^2 + im^2

MyArcTan[x_, y_, toplevel_] :=
	Which[
	NonNegValQ[x] && ZeroQ[y], 0,
	ConstantQ[x] && ConstantQ[y], MyArcTan2[x, y, toplevel],
	True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]
	]
	
MyArcTan2[x_, y_, toplevel_] :=
	Which[
		PosValQ[x], Which[y===x, Pi/4, y===-x, -Pi/4,
			y/x === Sqrt[3], Pi/3, x/y === Sqrt[3], Pi/6,
			y/x === -Sqrt[3], -Pi/3, x/y === -Sqrt[3], -Pi/6,
			True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]],
 		NegValQ[x], Which[ZeroQ[y], Pi, y===x, -3 Pi/4, y===-x, 3 Pi/4,
			y/x === Sqrt[3], -2 Pi/3, x/y === Sqrt[3], -5 Pi/6,
			y/x === -Sqrt[3], 2 Pi/3, x/y === -Sqrt[3], 5 Pi/6,
			True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]]],
		PosValQ[y] && ZeroQ[x], Pi/2,
		NegValQ[y] && ZeroQ[x], -Pi/2,
 		True, If [toplevel, ArcTan[x, y], ArcTanTemp[x, y]] ]
 
CESimplify[w_] := w
	
(* ****************************************************************
*
*     Real Functions
*
***************************************************************** *)

(*RealAbs is currently not used by ComplexExpand*)
RealAbs[x_]:=
  Which[ConstantQ[x], If[Re[N[x]]>=0,x,Foil[-1,x]],
        SameQ[SyntacticSign[x],1], Abs[x],
        True, Abs[Foil[-1,x]]]

RealLog[x_?NonNegValQ ^ y_?NumberQ] := y RealLog[x] /;
	Im[y] == 0 (*&& AbsExpr[x, toplevel]==x*)

RealLog[x_] := 
  If [ Head[x]===Power && x[[1]]===E,  x[[2]], Log[x]]
	
RealSign[x_]:=
  Which[ConstantQ[x], SignRe[N[x]],
        SameQ[SyntacticSign[x],1], Sign[x],
        True, -Sign[Foil[-1,x]]]

SyntacticSign[expr_]:=
  (* If the first character of EXPR that will print is "-" then -1,
     if SameQ[EXPR,0] then 0,
     otherwise 1. *)
  Module[{answer},
    answer=Switch[Head[expr],
                   Integer, Sign[expr],
		   Rational, Sign[expr],
		   Real, Sign[expr],
                   Complex, If[SameQ[Re[expr],0],
		               Sign[Im[expr]],
			       Sign[Re[expr]]],
		   Plus, SyntacticSign[expr[[1]]],
		   Times, SyntacticSign[expr[[1]]],
		   _, 1];
    answer]

Foil[expr1_,expr2_]:=
  (* Distributive product. *)
  Module[{i,j,answer},
    answer=If[Not[SameQ[Head[expr1],Plus]],
               If[Not[SameQ[Head[expr2],Plus]],
	          expr1*expr2,
	          Sum[expr1*expr2[[j]],
		      {j,Length[expr2]}]],
	       If[Not[SameQ[Head[expr2],Plus]],
	          Sum[expr1[[i]]*expr2,
		      {i,Length[expr2]}],
	          Sum[Sum[expr1[[i]]*expr2[[j]],
		          {j,Length[expr2]}],
		      {i,Length[expr1]}]]];
    answer]

(* ****************************************************************
*
*     Install our definitions of Abs, Arg, Re and Im
*
*     Info about the C Routines:
*(1) Abs -- does numbers.
*(2) Arg -- does infinities, real and imaginary exact numbers, and floats.
*(3) Conjugate -- does numbers.
*(4) Im -- does numbers.
*(5) Re -- does numbers.
*(6) Sign -- does infinities, numbers, constant symbols, integer Power's
*recursively, Plus recursively if all args have same sign, and 
*Times recursively.
***************************************************************** *)

ConstSymbQ[s_Symbol] := MemberQ[ConstantSymbols, s]

ConstantQ[s_Symbol] := MemberQ[ConstantSymbols, s]
ConstantQ[s_Symbol[x_]] := MemberQ[AllowedSix,s] && ConstantQ[x]
ConstantQ[_?NumberQ] := True
ConstantQ[_?ZeroQ] := True
ConstantQ[a_ + b_] := ConstantQ[a] && ConstantQ[b]
ConstantQ[a_ * b_] := ConstantQ[a] && ConstantQ[b]
ConstantQ[a_ ^ b_] := ConstantQ[a] && ConstantQ[b]
ConstantQ[_] := False

RealValQ[s_Symbol] := MemberQ[ConstantSymbols, s] || 
	(AtomQ[s] && !IsComplex[s])

RealValQ[s_Symbol[x_]] := MemberQ[RealFunctions,s] ||
	(MemberQ[ConjugateFunctions,s] && RealValQ[x])

RealValQ[Log[_?PosValQ]]:= True
RealValQ[ArcTan[a_,b_]]:= RealValQ[a] && RealValQ[b]
RealValQ[ArcTanTemp[_,_]]:= True
RealValQ[(_Real | _Rational | _Integer)] := True
RealValQ[x_ - I*Im[x_]] := True
RealValQ[I*Im[x_] - x_] := True
RealValQ[I*z_ + Im[z_]] := True
RealValQ[I*z_ - I*Re[z_]] := True
RealValQ[I*Re[z_] - I*z_] := True
RealValQ[a_ - Im[a_]] := True
RealValQ[-a_ + Im[a_]] := True
RealValQ[a_ + Literal[Conjugate[a_]]] := True
RealValQ[a_ + b_] := (RealValQ[a] && RealValQ[b]) (*|| a === -Im[b] ||
	-a === Im[b] || b===Conjugate[a]*)
RealValQ[a_ * Literal[Conjugate[a_]]] := True
RealValQ[a_ * b_] := (RealValQ[a] && RealValQ[b]) (*|| b===Conjugate[a]*) ||
	(ImagValQ[a] && ImagValQ[b])
RealValQ[a_?PosValQ ^ b_] := RealValQ[b]
RealValQ[a_^b_?IntegerQ] := RealValQ[a] || (ImagValQ[a] && EvenQ[b])
RealValQ[_] := False

ImagValQ[x_ - Re[x_]] := True
ImagValQ[Re[x_] - x_] := True
ImagValQ[a_ - Literal[Conjugate[a_]] ] := True
ImagValQ[-a_ + Literal[Conjugate[a_]] ] := True
ImagValQ[a_ + Literal[Conjugate[-a_]] ] := True
ImagValQ[a_ - Re[a]] := True
ImagValQ[-a_ + Re[a]] := True
ImagValQ[a_ + Re[-a]] := True
ImagValQ[a_ + b_ ] := (ImagValQ[a] && ImagValQ[b]) (*|| b === -Conjugate[a] ||
	-b === Conjugate[a] || b === -Re[a] || -b === Re[a] ||
	a === -Re[b] || -a === Re[b]*)
ImagValQ[a_ * b_] := (ImagValQ[a] && RealValQ[b]) ||
	(ImagValQ[b]&&RealValQ[a])
ImagValQ[_] := False

NonNegValQ[Literal[Conjugate[a_]]] := NonNegValQ[a]
NonNegValQ[Literal[Arg[_?NonNegValQ]]] := True
NonNegValQ[Literal[Re[_?ImagValQ]]] := True
NonNegValQ[Literal[Abs[_?ZeroQ]]] := True
NonNegValQ[Literal[Im[_?RealValQ]]] := True
NonNegValQ[Literal[Sign[a_]]] := NonNegValQ[a]
NonNegValQ[ArcTanTemp[_,_?ZeroQ]]:= True
NonNegValQ[a_ + b_] := NonNegValQ[a] && NonNegValQ[b]
NonNegValQ[a_ * Literal[Conjugate[a_]]] := True
NonNegValQ[a_ * b_] := (NonNegValQ[a] && NonNegValQ[b]) ||
	(NonPosValQ[a] && NonPosValQ[b]) (*|| b===Conjugate[a]*)
NonNegValQ[a_ ^ b_] := NonNegValQ[a] && RealValQ[b]
NonNegValQ[a_?RealValQ ^ b_?EvenQ] := True
NonNegValQ[Literal[Abs[_]]] := True
NonNegValQ[Literal[Im[_?RealValQ]]] := True

NonNegValQ[a_Real | a_Rational | a_Integer] := NonNegative[a]
NonNegValQ[a_?ConstSymbQ] := True

NonNegValQ[_] := False

NonPosValQ[a_Real | a_Rational | a_Integer] := Not[Positive[a]]
NonPosValQ[Literal[Arg[_?NonNegValQ]]] := True
NonPosValQ[Literal[Conjugate[a_]]] := NonPosValQ[a]
NonPosValQ[Literal[Sign[a_]]] := NonPosValQ[a]
NonPosValQ[Literal[Re[_?ImagValQ]]] := True
NonPosValQ[Literal[Abs[_?ZeroQ]]] := True
NonPosValQ[Literal[Im[_?RealValQ]]] := True
NonPosValQ[a_ + b_] := NonPosValQ[a] && NonPosValQ[b]
NonPosValQ[a_ * Literal[Conjugate[-a_]]] := True
NonPosValQ[-a_ * Literal[Conjugate[a_]]] := True
NonPosValQ[a_ * (-Literal[Conjugate[a_]])] := True
NonPosValQ[a_ * b_] := (NonPosValQ[a] && NonNegValQ[b]) ||
	(NonNegValQ[a] && NonPosValQ[b])
NonPosValQ[_] := False

ZeroQ[Literal[Arg[_?NonNegValQ]]] := True
ZeroQ[Literal[Im[_?RealValQ]]] := True
ZeroQ[Literal[Re[_?ImagValQ]]] := True
ZeroQ[Literal[(Abs | Conjugate | Sign)[_?ZeroQ]]] := True
ZeroQ[ArcTanTemp[_?ZeroQ,_?ZeroQ]]:= True
ZeroQ[x_] := NonNegValQ[x] && NonPosValQ[x]
ZeroQ[_] := False
NotZeroQ[x_] := ! ZeroQ[x]

Epsilon = 10^(-30)

PosValQ[a_?RealValQ] := NonNegValQ[N[a]-Epsilon]
NegValQ[a_?RealValQ] := NonPosValQ[N[a]+Epsilon]
PosValQ[_] := False
NegValQ[_] := False

AllowedSix={Re, Im, Abs, Arg, Conjugate, Sign}
AllowedFunctions=AllowedSix
RecursiveFunctions={Sign} (* for top-level *)

AssumedComplex={ _ }  (* for top-level *)

IsComplex[e_] := Or @@ (MatchQ[e, #]& /@ AssumedComplex)

PolarQ[x_] := IsOK[{Arg, Abs},x]
CartQ[x_]  := IsOK[{Re, Im},x]

IsOK[f_List,x_] := And @@ (IsOK[#, x]& /@ f )
IsOK[f_,x_] := MemberQ[AllowedFunctions, f]

IsAnyOK[f_List,x_] := Or @@ (IsOK[#, x]& /@ f )
IsRecursive[f_] := MemberQ[RecursiveFunctions, f]
IsConjugate[f_] := MemberQ[ConjugateFunctions, f]

(* symbols ComplexExpand, TargetFunctions
   have been created by internal code *)

Begin["`Private`"]

Unprotect[ComplexExpand]

Options[ComplexExpand] = {TargetFunctions -> AllowedSix}

ComplexExpand::exf = "Value of option TargetFunctions -> `1` does not
contain any of the allowed functions `2`."

ComplexExpand::rlst="Warning: ComplexExpand does not support `1`. It only supports the rule TargetFunctions -> <function or list of functions>."

ComplexExpand[e_List, args___] := ComplexExpand[#, args]& /@ e

ComplexExpand[expr_, opts___Rule] := ComplexExpand[expr, {}, opts]

ComplexExpand[expr_, cvars_List, opts___Rule] :=
	Block[{AssumedComplex, AllowedFunctions, RecursiveFunctions,
	       badrulelist, x, `y, `funs},
		AssumedComplex = cvars;
		RecursiveFunctions = AllowedSix;
		badrulelist=Select[{opts},
			Function[!MatchQ[#,TargetFunctions->_]]];
		Do[Message[ComplexExpand::rlst,badrulelist[[i]]],
			{i, Length[badrulelist]}];
		funs = TargetFunctions /. {opts} /. Options[ComplexExpand];
		If[Head[funs] =!= List, funs = List[funs]];
		AllowedFunctions = Select[funs, MemberQ[AllowedSix, #]&];
		If[ Length[AllowedFunctions] < 1,
			Message[ComplexExpand::exf, funs, AllowedSix];
			Return[expr] ];
		{x, y} = ReImExpr[expr, False];
		If[ x === ReImFail || y === ReImFail, expr,
		 Block[{BadFunctions, Badpatterns, pos, return = x + I y},
		    BadFunctions = Select[AllowedSix, 
			! MemberQ[AllowedFunctions, #]& ];
		    BadPatterns = Apply[Alternatives,
			Table[BadFunctions[[j]] ,
			{j, Length[BadFunctions]}]];
		    pos =  Position[return, BadPatterns[_]];
		RecursiveFunctions={Sign};(* for top-level *)
		AllowedFunctions=AllowedSix;
		return = MapAt[ComplexExpand[#, cvars, opts]&,
		    	return, pos] /. {ArcTanTemp->ArcTan};
		     CESimplify[return]
		     	] 
		]
	]


ComplexExpand[expr_, var_, opts___Rule] := ComplexExpand[expr, {var}, opts]

ComplexExpand[args___] := $Failed /;
	Message[General::argct, ComplexExpand, Length[{args}]]
	

Protect[ComplexExpand]

End[]

On[ General::shdw]

EndPackage[]

(* do not leave context on path *)

$ContextPath = DeleteCases[$Packages, "System`ComplexExpand`"]

Unprotect[$Packages]
$Packages = DeleteCases[$Packages, "System`ComplexExpand`"]
Protect[$Packages]

Null
