(* Title: Real Only *)

(* Author: George Beck *)

(* :Summary:
Two ideas are implemented here. Odd roots of negative numbers are defined to be negative, and unavoidable complex numbers are signaled by the symbol Nonreal, in a way similar to the built-in object Indeterminate.
*)

(* :Context: Miscellaneous`RealOnly` *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 3.0 *)

(* :Copyright: Copyright 1995-1999, Wolfram Research, Inc. *)

(* :History: Revised May 16, 1995 to fix the bug: if Solve or Roots yields
powers of negative numbers, the redefinition of Power kicks in before
there is a chance to see if any of the solutions involve complex numbers.
Thus, Solve[x^3==-8, x] gave {{x->-2},{x->-2},{x->-2}} instead of
the intended {{x->-2},{x->Nonreal},{x->Nonreal}}. Therefore the code for
Solve and Root was added at the bottom.

Revised November 6, 1997 to fix bug 26470.
Adam Strzebonski added two lines in the code
"Unprotect[Power]; (* added by adams Nov. 5, 1997 *)"
saying "This should fix the autoloading vs. setting DownValues[Power] problem."

*)

(* :Keywords:
algebra, root, radical
*)

(* :Source:
*)

(* :Warning:
1. Loading this package redefines Power, Solve, Roots, and $Post; complex results are modified.

2. Starting with cubics, there are polynomial equations whose solutions, when expressed in terms of radicals, essentially involve complex numbers, even though one or more of the solutions may ultimately be real. (Historically, this seems to have been the main reason for the adoption of complex numbers.)
An example is the cubic x^3 + 6 x^2 + 3 - 12 == 0, which has three real roots. For such cases this package will incorrectly claim that some solutions are not real, while in fact they are real. Perhaps the wisest course is to check with NSolve.
*)

(* :Limitation: Complex results are not modified for special functions. *)

(* :Discussion:
Two ideas are implemented here. Odd roots of negative numbers are defined to be negative, and calculations with unavoidable complex numbers are condensed to the symbol Nonreal. This is done by redefining the built-in functions Power and $Post. The effect of the redefinition of Power is held off until after the nature of the results of Solve or Roots have been taken into account.
*)

BeginPackage["Miscellaneous`RealOnly`"]

Nonreal::usage = "Nonreal is a symbol that replaces a calculation result that involves an unavoidable complex number.";

Begin["`Private`"]

protected = Unprotect[Power];
Power[b_?Negative, Rational[m_, n_?OddQ]] := (-(-b)^(1/n))^m
Protect[Evaluate[protected]];

(* Please read the comment in Solve before changing this modification of Power. *)

NonrealRule::usage = "NonrealRule is a rule that, when applied to an expression, drops small imaginary parts, replaces large imaginary parts by the object Nonreal, and forces Nonreal up through most elementary calculations.";

NonrealRule =
{
	Complex[x_, _?(Chop @ # == 0&)] :> x,
	Complex[_, _?(Chop @ # =!= 0&)] -> Nonreal,
	(
		Plus | Times | Minus | Subtract |
		Sqrt | Power | PowerMod |
		Abs | Exp | Log |
		Sin | Cos | Tan |
		Cot | Sec | Csc |
		ArcSin | ArcCos | ArcTan |
		ArcCot | ArcSec | ArcCsc |
		Floor | Ceiling | Round |
		Mod | Quotient | Prime |
		Min | Max | LCM | GCD |
		Random | Rationalize	
	)
		[___, Nonreal, ___] :> Nonreal
};

Nonreal::warning = "Nonreal number encountered.";

NonrealAux[result_] := Module[
	{preliminary},
	preliminary = (result //. NonrealRule);
	If[
		Not[FreeQ[preliminary, Nonreal,
			{0, Infinity}, Heads -> True]],
		Message[Nonreal::warning]
	];
	preliminary
]

protected = Unprotect[Solve, Roots]

Solve[args___] :=
     Block[{$InsideSolve = True,
            powerVals = DownValues[Power],
            protected,
            result
           },
           protected = Unprotect[Power];
           (* Note that the presence of the symbol b in our rule for Power
              is used as the flag to identify the rule. This will break if
              that symbol name is ever changed.
           *)
           DownValues[Power] = DeleteCases[DownValues[Power],
                                           _?(!FreeQ[#, b]&)];
           result = Solve[args] /.
                    (x_ -> _?(!FreeQ[Chop[N[#]], _Complex]&)) :> x -> Nonreal;
	   Unprotect[Power]; (* added by adams Nov. 5, 1997 *)
           DownValues[Power] = powerVals;
           Protect[Evaluate[protected]];
           result
    ] /; $InsideSolve =!= True

(* In Roots, the only differences from Solve are
      the replacement of Solve by Roots and of "->" by "==". *)

Roots[args___] :=
     Block[{$InsideRoots = True,
            powerVals = DownValues[Power],
            protected,
            result
           },
           protected = Unprotect[Power];
           DownValues[Power] = DeleteCases[DownValues[Power],
                                           _?(!FreeQ[#, b]&)];
           result = Roots[args] /.
                 (x_ == _?(!FreeQ[Chop[N[#]], _Complex]&)) :> x == Nonreal;
	   Unprotect[Power]; (* added by adams Nov. 5, 1997 *)
           DownValues[Power] = powerVals;
           Protect[Evaluate[protected]];
           result
    ] /; $InsideRoots =!= True

Protect[Evaluate[protected]];

If[
	FreeQ[$Post, NonrealAux],
	$Post = Composition[NonrealAux, If[ValueQ[$Post], $Post, Identity]]
];

End[]
EndPackage[]


