(*
:Titulo / Title: RLSCF
*)
(*
Representación y localización de singularidades y ceros de funciones reales 
en 2D.
*)
(*
Representation and location of singularities and zeros of 2D real functions
*)
(*
:Autor / Author: Fernando Gomez Lanza
*)
(*
:email:  fgomez@xtec.cat
*)
(*
:URL:  http://www.xtec.cat/~fgomez
*)
(*
Version: 1.5.1
*)
(*
:Fecha / Date: March 2008
*)
(*
:Requisitos / Requirements: Mathematica 3.0 or later ,   InequalitySolve.m 
*)
(*
:Nota: El análisis y representación de funciones necesita Mathematica 
3.0 o superior. Este paquete proporciona buenos resultados para las funciones 
reales 2D más comunes, como, polinómicas o fracciones algebraicas (incluyendo 
Hermite, Laguerre y cualquier otro tipo de polinomio), raíces, 
trigonométricas e inversas, hiperbólicas e inversas, exponenciales y 
logarítmicas, utilizadas también dentro de Which[...] cuando sean definidas por intervalos. 
No usar expresiones Which[...]/Which[...]. Las funciones que 
presentan iteraciones como Tan[Tan[x]], y combinaciones trigonométricas, 1/(
Tan[x]+Sin[x]-x), dan el mismo resultado incorrecto que Plot, al aparecer 
ecuaciones trascendentes. En este caso debemos proceder manualmente. Podemos 
trabajar con otras funciones no definidas añadiendo sus propiedades a 
ListaCeros y ListaInfinitos. Ver el cuaderno de documentación para 
procedimiento manual y otras explicaciones. Este paquete no usa el método 
de la derivada segunda para localizar las discontinuidades.
*)
(*
:Note: The analysis and representation of functions needs Mathematica 3.0 or 
later. This package give us good results for common 2D real functions, such 
as, polynomial or algebric fraction (including Hermite, Laguerre, and all 
other polynomial type), roots, trigonometric and inverses, hiperbolic and 
inverses, exponential and logaritmic functions, using Which[..] also. 
Not use Which[..]/Which[...] expressions. 
Functions with iterations, such as, Tan[Tan[x]], and trigonometric 
combinations, 1/(Tan[x]+Sin[x]-x), give the same incorrect result with Plot, 
because appear trascendental equations. In this case we need to procede 
manually. We can work with other non defined functions adding the properties 
to ListaCeros and ListaInfinitos. See documentation notebook for manual 
procedure and other explanations. This package doesn't use the method of the 
one derived second to locate the discontinuities.
*)
BeginPackage["RLSCF`"];

ListaCeros::usage=
  "ListaCeros={...,{f,v,p},...}  \n\n Lista de funciones analizadas para \
encontrar ceros. f=head-function, v=valor que cumple f[v]=0, p=periodo de la \
función f (0 si no es periodica). Ejemplo {...,{Tan,0,Pi},...}. \n\n List of \
analysed functions to find zeros. f=head-function, v=value that verifies \
f[v]=0, p=function's period (0 if not period). Example {...,{Tan,0,Pi},...} \
."; 
  
ListaInfinitos::usage=
  "ListaInfinitos={...,{f,v,p},...}  \n\n Lista de funciones analizadas para \
encontrar infinitos. f=head-function, v=valor que cumple f[v]=0, p=periodo de \
la funcion f (0 si no es periodica). Ejemplo {...,{Tan,Pi/2,Pi},...}. \n\n \
List of analysed functions to find infinities. f=head-function, v=value that \
verifies f[v]=Infinity, p=function's period (0 if not period). Example \
{...,{Tan,Pi/2,Pi},...} .";

Ceros::usage=
  "Ceros[f[x],{x,xmin,xmax}] \n\n Calcula los ceros de la función f[x]  entre \
xmin y xmax. \n\n Calculate zeros of function f[x] from xmin to xmax.";

Infinitos::usage=
  "Infinitos[f[x],{x,xmin,xmax}] \n\n Calcula los infinitos que aparecen en \
la función  f[x] entre xmin y xmax. \n\n Calculate infinities of function \
f[x] from xmin to xmax.";

Argum::usage=
  "Argum[f,expr] \n\n Encuentra los argumentos de la función f[ ] contenidas \
en la expresión expr. \n\n Finds the arguments of function f[ ] in expr. \n\n \
Argum[Tan,3x+5Tan[4x/8]].";

ArguWhich::usage=
  "ArguWhich[expr] \n\n Encuentra los argumentos de la función Which en la \
expresión expr. \n\n Finds the arguments of function Which in expr.";
	
LSFun::usage=
  "LSFun[f[x],{x,xmin,xmax}] \n\n Proporciona la lista de posibles \
singularidades de la función f[x] en el intervalo xmin a xmax. \n\n Provide \
the list of probable singularities of function f[x] from xmin to xmax. \n\n \
LSFun[Tan[x],{x,-5,5}]";

DibujarIntervalos::usage=
  "DibujarIntervalos[f[x],{x,xmin,xmax},opciones] \n\n Representa la función  \
f[x]  entre xmin y xmax. No usar con raíces. \n\n Plots f[x] from \
xmin to xmax. Not use with roots.";

SoloR::usage=
  "SoloR[f[x]] \n\n Transforma la función f[x] para que considere solo los \
valores reales, no considera las raíces imaginarias de x^(1/3), por ejemplo. \
\n\n Transform the function f[x] to consider real values only, example: \
complex roots of x^(1/3) are not considered.";

PlotFR::usage=
  "PlotFR[f[x],{x,xmin,xmax}] \n\n Representa la función \
f[x] desde xmin a xmax, evitando las discontinuidades. \n\n Plots the function f[x] from xmin to xmax, avoiding discontinuities.";
 
vlkmax::usage=
  "vlkmax=10 \n\n Cuando aparece un número infinito de \
discontinuidades PlotFR solo muestra vlkmax. \n\n When appear an infinite \
number of discontinuities PlotFR shows vlkmax only.";

vlkmax=10;
		
ListaCeros={{Null,0,0},{Log,1,0},{Tan,0,Pi},{Cot,Pi/2,Pi},{Sin,0,Pi},{Cos,
      Pi/2,Pi},{Sinh,0,0},{Tanh,0,0}};

ListaInfinitos={{Log,0,0},{Tan,Pi/2,Pi},{Cot,0,Pi},{Csc,0,Pi},{Sec,Pi/2,Pi},{
      Csch,0,0},{Coth,0,0}};

Unprotect[Log];
Log[Global`x_/Global`y_]=Log[Global`x]-Log[Global`y];
Log[Global`b_,Global`x_/Global`y_]=
  Log[Global`b,Global`x]-Log[Global`b,Global`y];
Protect[Log];
 
Begin["`privado`"];

AF[f_,var_]:=
  If[AtomQ[f]&&(NumberQ[f]==False)&&(Head[f]=!=String)&&(f=!=var),f,
    Function[var,f]];

Argum[f_,expr_]:=
  If[f===Null,{expr},
    Union[Map[Apply[Part,Flatten[{expr,#,1}]]&,Position[expr,f[_]]]]];

ArguWhich[f_]:=
  Flatten[Map[Partition[Apply[List,Apply[Part,#]],2]&,
      Map[Drop[Prepend[#,f],-1]&,Position[f,Which]]],1];

DescomponerWhich[f_/;FreeQ[f,Which]==False]:=
  ReplaceAll[f,Which[__,___]->EliRepes[Map[Last[#]&,ArguWhich[f]]]];

ArguPower[f_]:=
  Flatten[Map[Apply[List,#]&,
      DeleteCases[
        Map[Apply[Part,#]&,Map[Drop[Prepend[#,f],-1]&,Position[f,Power]]],
        Power[_,_Integer]]],1];

Needs["Algebra`InequalitySolve`"];

Kintervalo[expr_,k_,l_List:{-10,10}]:=
  Module[{lista,lmax,lmin},
    lista=Apply[List,
        Flatten[Map[InequalitySolve[l[[1]]<=#<=l[[2]],k]&,Flatten[{expr}]]]];
    lista=N[lista];
	 lmax=Floor[
        Max[Select[
            Map[Apply[Part,#]&,
              Union[Map[Append[#,2]&,
                  Map[Prepend[#,lista]&,Position[lista,k==_]]],
                Map[Append[#,2]&,Map[Prepend[#,lista]&,Position[lista,k<=_]]],
                Map[Append[#,5]&,
                  Map[Prepend[#,lista]&,Position[lista,_Inequality]]],
                Map[Append[#,3]&,
                  Map[Prepend[#,lista]&,Position[lista,_<=k<=_]]]]],
            Im[#]==0&]]];
    lmin=Ceiling[
        Min[Select[
            Map[Apply[Part,#]&,
              Union[Map[Append[#,2]&,
                  Map[Prepend[#,lista]&,Position[lista,k==_]]],
                Map[Append[#,2]&,Map[Prepend[#,lista]&,Position[lista,k>=_]]],
                Map[Append[#,1]&,
                  Map[Prepend[#,lista]&,Position[lista,_Inequality]]],
                Map[Append[#,1]&,
                  Map[Prepend[#,lista]&,Position[lista,_<=k<=_]]]]],
            Im[#]==0&]]];
    If[lmin<=lmax,{{lmin,lmax}},{{lmax-vlkmax,lmax},{lmin,lmin+vlkmax}}]];

EliRepes[l_List]:=Union[l,l];

LSGen[ff_,l_List:{Global`x,-10,10},g_List:{Null,0,0}]:=
  Module[{p,s,f,x,k,lk},x=l[[1]];f=AF[ff,l[[1]]];p=Argum[g[[1]],f[x]];
    p=ReplaceAll[
        p,{Tan[a_]->Tan[a+k*Pi],Sin[a_]->Sin[a+2Pi*k],Cos[a_]->Cos[a+2Pi*k],
          Sec[a_]->Sec[a+2Pi*k],Csc[a_]->Csc[a+2Pi*k],Cot[a_]->Cot[a+Pi*k]}];
    s=Flatten[Map[Solve[#==g[[2]]+k*g[[3]],x]&,p]];s=Map[(x/.#)&,s];
    s=Chop[N[s],10^(-$MachinePrecision+1)];
    lk=If[FreeQ[s,k],{{1,1}},Kintervalo[s,k,{l[[2]],l[[3]]}]];
    EliRepes[Select[
        N[Flatten[Union[Map[Table[s,{k,#[[1]],#[[2]],1}]&,lk]]]],(
            l[[2]]<=#<=l[[3]])&]]];

Ceros[f_,l_List:{Global`x,-10,10}]:=
  EliRepes[Flatten[Map[LSGen[f,l,#]&,ListaCeros]]];

Infinitos[f_,l_List:{Global`x,-10,10}]:=
  EliRepes[Flatten[Map[LSGen[f,l,#]&,ListaInfinitos]]];

LsLsum[ff_,l_List:{Global`x,-10,10}]:=
  Module[{lslsum,lsum,ls,f,x},x=l[[1]];f=AF[ff,x][x];lsum=Apart[f];
    lsum=Switch[lsum,_Plus,Apply[List,lsum],_,{lsum}];
			ls=If[FreeQ[f,Which],{},
        EliRepes[
          Select[Flatten[Apply[List,Map[First[#]&,ArguWhich[f]],10]],
            NumberQ]]];
	  lsum=Flatten[Map[If[FreeQ[#,Which],#,DescomponerWhich[#]]&,lsum]];
    lsum=Union[lsum,Flatten[Map[ArguPower[#]&,lsum],1]];
		lsum=Union[lsum,Flatten[Map[Apply[List,#]&,lsum]]];
	 lslsum=Map[LsLsum[#,l]&,Cases[lsum,x_+y_]];
	lsum=DeleteCases[lsum,x_+y_];
	{Union[ls,Flatten[Map[First[#]&,lslsum]]],
      Union[lsum,Flatten[Map[Part[#,2]&,lslsum]]]}];

LSFun[ff_,l_List:{Global`x,-10,10}]:=
  Module[{lslsum,lsum,ls,f,x},x=l[[1]];
	  lslsum=LsLsum[ff,l];
		ls=First[lslsum];
		lsum=Last[lslsum];
		EliRepes[
      Flatten[{ls,Map[Ceros[#,l]&,EliRepes[Map[Denominator[#]&,lsum]]],
          Map[Infinitos[#,l]&,EliRepes[Map[Numerator[#]&,lsum]]]}]]];

DibujarIntervalos[f_,l_List:{Global`x,-10,10},opciones___]:=
  Module[{linter,lsingu,ff,x},x=l[[1]];ff=AF[f,x];lsingu=Evaluate[LSFun[f,l]];
    linter=Map[Prepend[#,x]&,
        Partition[Sort[Union[N[Drop[l,1]],N[lsingu]]],2,1]];
    Show[Map[Plot[ff[x],#,DisplayFunction->Identity,opciones]&,linter],
      DisplayFunction->$DisplayFunction]];

SoloR[f_]:=
  f//.a_^n_:>
      ((Sign[a])^Numerator[Rationalize[n]])*Abs[a]^
            n/;((Head[a]=!=Abs)&&(Mod[Denominator[Rationalize[n]],2]=!=0)&&(
              IntegerQ[n]==False));

PlotFR[f_,l_List:{Global`x,-10,10},opciones___]:=
  Module[{linter,lsingu,ff,x},x=l[[1]];ff=AF[f,x];lsingu=Evaluate[LSFun[f,l]];
    linter=Map[Prepend[#,x]&,
        Partition[Sort[Union[N[Drop[l,1]],N[lsingu]]],2,1]];
    Show[Map[Plot[Evaluate[SoloR[ff[x]]],#,DisplayFunction->Identity,
            opciones]&,linter],DisplayFunction->$DisplayFunction]];

End[];

EndPackage[];
