# This is a program for real root isolation of multivariate real nonlinear systems
# Before use our program, one need to load the Maple package intpakX for interval computation
# We provide two commands:
# 1: [R,SR]=rootfinding_cotrans(sys,var,B0,err,k)
# 2: [R,SR]=rootfindingbox(sys,var,B0,err,k)
# The first one is based on coordinate transformation and then we need only to find real roots in subboxes of [-1,1]*[-1,1]
# The second one is that finding real roots in the given box B0
#######################
# Parameters:
# Input:
#       sys:  A multivariate real nonlinear systems (f1,f2,...,fn)
#       var:  list of variables
#       B0:   An initial box 
#       err:  A termination precision and err>0
#       k:    (optional, default=0) An integer for spliting the initial box into 2^k parts before computation
# Output:
#       R:    The isolating box set of the given systen
#       SR:   The suspected root box set of the given systen
######################
# Examples:
# libname:="D:\\intpakX", libname;
# read("D:\\rootfinding_analyticfun.txt");
#f1:=exp(x+y)+x^2-3;;f2:=cos(y)-x^2;
#A:=rootfinding_cotrans([f1,f2],[x,y],[[-10,10],[-10,10]],10^(-4));
#         [[[[0.750000000000000, 1.00000000000000], 
#           [-0.500000000000000, -0.250000000000000]], [
#           [-0.500000000000000, -0.375000000000000], 
#           [1.39130434782609, 1.50943396226415]]], []]
#    R:=A[1];SR:=A[2];



with(RootFinding):
with(ListTools):
with(VectorCalculus):
with(LinearAlgebra):
with(intpakX):
with(combinat):
Digits:=15:



rootfinding_cotrans:=proc(sys,var,B0,err)
local uB,n,Q,T,i,j,B,R,SR,t,L,temp,F,k;
if nargs<5 then
   k:=0;
else
   k:=args[5];
fi;
if nops(sys)<>2 or nops(var)<>2 then
  ERROR("The input system must be a one with two functions and two variables!");
fi;
R:=[];SR:=[];
L:=allbinary(n);
uB:=unitbox(B0);
for i from 1 to 4 do
    temp:=transformsys(sys,var,uB,L[i]);
    F:=temp[1];B:=temp[2];
    for j from 1 to nops(B) do
        if k=0 then 
           T:=rootfindingbox(F,var,B[j],err);
        else
           T:=rootfindingbox(F,var,B[j],err,k);
        fi;
        Q:=transformbox(T[1],L[i]);
        R:=[op(R),op(Q)];
        Q:=transformbox(T[2],L[i]);
        SR:=[op(SR),op(Q)];
    od;
od;
T:=SR;SR:=[];
for i from 1 to nops(T) do
    Q:=rootfindingbox(sys,var,T[i],err);
    R:=[op(R),op(Q[1])];
    SR:=[op(SR),op(Q[2])];
od;
return [R,SR];
end:


transformbox:=proc(B,L)
local n,Q,T,i,j,t;
n:=nops(L);Q:=[];
for i from 1 to nops(B) do;
    T:=B[i];
    for j from 1 to nops(L) do
        if L[j]=1 then
           t:=[1/T[j][2],1/T[j][1]];
           T[j]:=t;
        fi;
    od;
    Q:=[op(Q),T];
od;
return Q;
end:

unitbox:=proc(B)
local n,i,V,Q,T;
n:=nops(B);V:=[];
for i from 1 to n do
    T:=B[i];
    if T[1]<-1 and T[2]>1 then
       Q:=[[-1,1/T[1]],[-1,1],[1/T[2],1]];
    elif T[1]<-1 then
       if T[2]<=-1 then
          Q:=[[1/T[2],1/T[1]],[],[]];
       else
          Q:=[[-1,1/T[1]],[-1,T[2]],[]];
       fi;
    elif T[2]>1 then
       if T[1]>=1 then
          Q:=[[],[],[1/T[2],1/T[1]]];
       else
          Q:=[[],[T[1],1],[1/T[2],1]];
       fi;
    else
       Q:=[[],T,[]];
    fi;
    V:=[op(V),Q];
od;
return V;
end:

transformsys:=proc(sys,var,B,L)
local n,F,T,Q,t,i,j,k,a;
n:=nops(sys);F:=sys;Q:=[];k:=1;
if L[1]=0 then
   if B[1][2]=[] then
      k:=0;
   else
      Q:=[[B[1][2]]];
   fi;
fi;
if L[1]=1 then
   F:=subs(var[1]=1/var[1],F);
   if B[1][1]=[] and B[1][3]=[] then
      k:=0
   else
      if B[1][1]<>[] then
         Q:=[[B[1][1]]];
      fi;
      if B[1][3]<>[] then
         Q:=[op(Q),[B[1][3]]];
      fi; 
   fi;
fi;
i:=2;
while i<=n and k=1 do
    if L[i]=0 then
      if B[i][2]=[] then
         k:=0;
      else
         T:=[];
         for j from 1 to nops(Q) do
             a:=[op(Q[j]),B[i][2]];
             T:=[op(T),a];
         od;
         Q:=T;
      fi;
    fi;
    if L[i]=1 then
       F:=subs(var[i]=1/var[i],F);
       if B[i][1]=[] and B[i][3]=[] then
          k:=0
       else
          T:=[];
          if B[i][1]<>[] then
             for j from 1 to nops(Q) do
                 a:=[op(Q[j]),B[i][1]];
                 T:=[op(T),a];
              od;
          fi;
          if B[i][3]<>[] then
             for j from 1 to nops(Q) do
                 a:=[op(Q[j]),B[i][3]];
                 T:=[op(T),a];
             od;
          fi;
          Q:=T;
       fi; 
     fi;
     i:=i+1
od;
F:=simplify(F);
F:=numer(F);
if k=0 then
   Q:=[];
fi;
Q:=evalf(Q);
return [F,Q];
end:



transformsys_old:=proc(sys,var,N,L)
local n,F,B,T,Q,t,i,j;
n:=nops(sys);F:=sys;B:=[[seq([-1.,1.],i=1..n)]];
for i from 1 to n do
    if L[i]=1 then
       F:=subs(var[i]=1/var[i],F);
       Q:=[];
       for j from 1 to nops(B) do
           T:=B[j];
           t:=evalf([-1,-1/N]);T[i]:=t;
           Q:=[op(Q),T];
           t:=evalf([1/N,1]);T[i]:=t;
           Q:=[op(Q),T];
       od;
       B:=Q;
     fi;
od;
F:=simplify(F);
return [F,B];
end:

allbinary:=proc(n)
local i,j,T,Q;
T:=[[0],[1]];
for i from 2 to n do
    Q:=[];
    for j from 1 to nops(T) do
        Q:=[op(Q),[op(T[j]),0]];
        Q:=[op(Q),[op(T[j]),1]];
    od;
    T:=Q;
od;
return T;
end:

    
rootfindingbox:=proc(sys,var,B0,err)
local BS,T,R,SR,n,a1,a2,a3,w,B,F,m,sys_horner,F_horner;
if nargs<5 then
   BS:=evalf([B0]);
else
   BS:=splitB(evalf(B0),args[5]);
fi;
R:=[];SR:=[];
n:=nops(BS);m:=nops(BS[1]);
sys_horner:=convert(sys,horner,var);
while n>0 do
    B:=BS[1];
    BS:=subsop(1=NULL,BS);
    a1:=Exclusion(sys_horner,var,B);
    if a1=0 then
       F:=preconditioner(sys,var,B);
       #F_horner:=convert(F,horner,var);
       #a2:=IsOMsys(F_horner,var,B);
       a2:=IsOMsys(F,var,B); 
       if a2=1 then
          a3:=Existence(F,var,B);
          if a3=1 then
             R:=[op(R),B];
          elif a3=1/2 then
             SR:=[op(SR),B];
          fi;
       else
           w:=widthB(B,m);
           if w>err then
              T:=splitB(B,2);
              BS:=[op(BS),op(T)];
            else
              SR:=[op(SR),B];
            fi;
        fi;
    fi;
    n:=nops(BS);
od;
return [R,SR];
end:



middleB:=proc(B,n)
local p,i;
p:=[];
for i from 1 to n do
    p:=[op(p),evalf((B[i][1]+B[i][2])/2)];
od;
return p;
end:

widthB:=proc(B,n)
local p,i;
p:=[];
for i from 1 to n do
    p:=[op(p),B[i][2]-B[i][1]];
od;
p:=sort(p);
return p[n];
end:

mysubs:=proc(f,var,p)
local A,a,i;
A:=[];
for i from 1 to nops(p) do
    A:=[op(A),var[i]=p[i]];
od;
a:=subs(op(A),f);
return evalf(a);
end:


intervalfun:=proc(f,var,B)
local F,A;
if type(f,numeric) then
   A:=evalf([f,f]);
else
   F:=inapply(f,op(var));
   A:=F(op(B));
fi;
return A;
end:

intervalfun_mean:=proc(f,var,B)
local n,p,fp,A,a,i;
n:=nops(B);p:=middleB(B,n);
fp:=mysubs(f,var,p);
a:=intervalfun(diff(f,var[1]),var,B) &* (B[1] &- p[1]);
A:=fp &+ a;
for i from 2 to n do
    a:=intervalfun(diff(f,var[i]),var,B) &* (B[i] &- p[i]);
    A:=A &+ a;
od;
return A;
end:


Exclusion:=proc(F,var,B)
local ind,i,A;
ind:=0;i:=1;
while i<=nops(F) and ind=0 do
      A:=intervalfun(F[i],var,B);
      if A[1]>0 or A[2]<0 then
         ind:=1;
      fi;
      i:=i+1;
od;
return ind;
end:

generateOMmatrix:=proc(n)
local i,U,k;
k:=0;
while k=0 do
     U:=RandomMatrix(n,generator=-1. .. 1. );
     for i from 1 to n do
           U[[i],[i]]:=sign(U[i][i])*n;
     od;
     k:=IsOMmatrix(U,n);
od;
return U;
end:

Isnumericmatrix:=proc(M,n)
local i,j,t;
t:=1;i:=1;j:=1;
while i<=n and t=1 do
     while j<=n and t=1 do
           if type(M[i][j],numeric) then
           else
               t:=0;
           fi;
           j:=j+1;
      od;
      i:=i+1;
od;
return t;
end:

preconditioner:=proc(sys0,var,B)
local p,Jp,A,sys,U,n,a;
n:=nops(B);
p:=middleB(B,n);
Jp:=Jacobian(sys0 ,var = p);
if Rank(Jp)=n then
      if n=2 then
         U:=Matrix(2,2,[[2,1],[1,-2]]);
      elif n=3 then
         U:=Matrix(3,3,[[3,1,1],[1,-3,1],[1,1,3]]);
      elif n=4 then
         U:=Matrix(4,4,[[4,1,1,-1],[1,-4,1,1],[1,1,4,1],[1,1,1,4]]);
      elif n=5 then
         U:=Matrix(5,5,[[5,1,1,-1/2,1],[1,-5,1,-1,-1],[-1,1,5,1,-1],[1,1,1,5,-1],[1,1,1,1,5]]);
      else
         U:=generateOMmatrix(n);
      fi;
      A:=U.MatrixInverse(Jp);
      a:=Isnumericmatrix(A,n);
      if a=0 then
         A:=Matrix(n, shape = identity);
      fi;
else
   A:=Matrix(n, shape = identity);
fi;
sys:=A.<sys0>;
sys:=convert(sys,list);
return sys;
end:

IsOMmatrix:=proc(M,n)
local ind,i,j,W,T;
ind:=1;i:=1;
while i<=n and ind=1 do
      j:=1;
      while j<=n and ind=1 do
            if M[i][j]=0 then
               ind:=0;
            fi;
            j:=j+1;
      od;
      i:=i+1;
od;
i:=n;
while i>=2 and ind=1 do
      T:=choose(n,i);
      j:=1;
      while j<=nops(T) and ind=1 do
            W:=Determinant(M[[1..i],T[j]]);
            if W=0 then
               ind:=0;
            fi;
            j:=j+1;
      od;
      i:=i-1;
od;
return ind;
end:


IsOMsys:=proc(F,var,B)
local ind,i,j,S,T,B0,f,svar,W,M,J,n;
ind:=1;i:=1;J:=Jacobian(F,var);n:=nops(F);
M:=Matrix(n,n);
while i<=n and ind=1 do
      j:=1;
      while j<=n and ind=1 do
            M[[i],[j]]:=intervalfun(J[i][j],var,B);
            if M[i][j][1]<=0 and  M[i][j][2]>=0 then
               ind:=0;
            fi;
            j:=j+1;
      od;
      i:=i+1;
od;
i:=n;
while i>=2 and ind=1 do
      S:=Matrix(i,i,symbol=s);
      f:=Determinant(S);
      svar:=convert(S,list);
      T:=choose(n,i);
      j:=1;
      while j<=nops(T) and ind=1 do
            B0:=convert(M[[1..i],T[j]],list);
            W:=intervalfun(f,svar,B0);
            if W[1]<=0 and  W[2]>=0 then
               ind:=0;
            fi;
            j:=j+1;
      od;
      i:=i-1;
od;
return ind;
end:

splitB:=proc(B,z)
local n,T,m,Q,temp,i,j,a,k,w;
n:=nops(B);T:=[B];
for i from 1 to n do
    m:=nops(T);w:=(B[i][2]-B[i][1])/z;Q:=[];
    for j from 1 to m do
        for k from 0 to z-1 do;
            a:=[B[i][1]+k*w,B[i][1]+(k+1)*w];
            temp:=subsop(i=a,T[j]);
            Q:=[op(Q),temp];
        od;
    od;
    T:=Q;
od;
return T;
end:

vertexB:=proc(B)
local p,n,P,i,j;
n:=nops(B);p:=[[B[1][1]],[B[1][2]]];
for i from 2 to n do
    P:=[];
    for j from 1 to nops(p) do
        P:=[op(P),[op(p[j]),B[i][1]]];
        P:=[op(P),[op(p[j]),B[i][2]]];
    od;
    p:=P;
od;
return p;
end:


mysubs2:=proc(f,var,p)
local n,q,ind,i;
n:=nops(p);
q:=mysubs(f,var,p[1]);
if q=0 then
   ind:=0;
else
   ind:=sign(q);
fi;
i:=2;
while i<=n and ind<>0 do
      q:=mysubs(f,var,p[i]);
       if q=0 then
          ind:=0;
       elif sign(q)<>ind then
          ind:=0;
       fi;
       i:=i+1;
od;
return ind;
end:

enlargeB:=proc(B0,B)
local n,w1,w2,A,i;
n:=nops(B);A:=[];
for i from 1 to n do
    w1:=(B0[i][1]-B[i][1])/2;
    w2:=(B[i][2]-B0[i][2])/2; 
    A:=[op(A),[B0[i][1]-w1,B0[i][2]+w2]];
od;
return A;
end:

Existence:=proc(F,var,B)
local n,i,T,G,B0,var0,vB,s1,s2,ind,A,w,setB,bigB,k;
n:=nops(B);
if n=2 then
   A:=Existence_2d(F,var,B);
   ind:=A[1];
else
   i:=n;T:=[];k:=0;
   #### find two systems in faces ####
   while  i>=1 and nops(T)<2 do
         G:=subs(var[i]=B[i][1],F);
         B0:=subsop(i=NULL,B);
         var0:=subsop(i=NULL,var);
         A:=Existence(G[1..n-1],var0,B0);
         if A=1 then
            T:=[op(T),[G,var0,B0]];
         elif A=1/2 then
            k:=1;
         fi;
         G:=subs(var[i]=B[i][2],F);
         A:=Existence(G[1..n-1],var0,B0);
         if A=1 then
            T:=[op(T),[G,var0,B0]];
         elif A=1/2 then
            k:=1;
         fi;          
         i:=i-1;
    od;
    ##### compute teh evaluation of fn in two faces ######
    if nops(T)=2 then
       G:=T[1][1];var0:=T[1][2];B0:=T[1][3];
       w:=widthB(B0,n-1);setB:=[B0];
       vB:=vertexB(B0);
       s1:=mysubs2(G[n],var0,vB);
       while s1=0 and w>10^(-6) do
             i:=1;setB:=splitB(B0,2);A:=0;
             while i<=nops(setB) do
                   A:=Existence(G[1..n-1],var0,setB[i]);
                   if A=1 then
                      B0:=setB[i];
                      i:=nops(setB)+1;
                   elif A=1/2 then
                      bigB:=enlargeB(setB[i],B0);
                      A:=Existence(G[1..n-1],var0,bigB);
                      if A=1 then
                         B0:=bigB;
                         i:=nops(setB)+1;
                      fi;
                   fi;
                   i:=i+1;
             od;  
             if A=1 then
                vB:=vertexB(B0);
                s1:=mysubs2(G[n],var0,vB);
                if s1=0 then
                   w:=w/2;
                fi;
             else
                w:=10^(-7);
             fi;         
       od;
       G:=T[2][1];var0:=T[2][2];B0:=T[2][3];
       w:=widthB(B0,n-1);setB:=[B0];
       vB:=vertexB(B0);
       s2:=mysubs2(G[n],var0,vB);
       while s2=0 and w>10^(-6) do
             i:=1;setB:=splitB(B0,2);A:=0;
             while i<=nops(setB) do
                   A:=Existence(G[1..n-1],var0,setB[i]);
                   if A=1 then
                      B0:=setB[i];
                      i:=nops(setB)+1;
                   elif A=1/2 then
                      bigB:=enlargeB(setB[i],B0);
                      A:=Existence(G[1..n-1],var0,bigB);
                      if A=1 then
                         B0:=bigB;
                         i:=nops(setB)+1;
                      fi;
                   fi;
                   i:=i+1;
             od;  
             if A=1 then
                vB:=vertexB(B0);
                s2:=mysubs2(G[n],var0,vB);
                if s2=0 then
                   w:=w/2;
                fi;
             else
                w:=10^(-7);
             fi;         
       od;
       if s1*s2<0 then
          ind:=1;
       elif s1*s2>0 then
          ind:=0;
       else
          ind:=1/2;
       fi;
    elif k=0 then
         ind:=0;
    else
       ind:=1/2;       
    fi;
fi;
return ind;
end:













Findindex:=proc(A,a)
local n,i,t;
n:=nops(A);i:=1;t:=0;
while i<=n and t=0 do
  if A[i]=a then
     t:=i;
  fi;
  i:=i+1;
od;
return t;
end:


whichcase:=proc(f,g)
local sf,sg,t,ind;
sf:=sort(f);
sg:=sort(g);
if sf[1]>0 or sf[4]<0 or sg[1]>0 or sg[4]<0 then
   ind:=0;

elif (sf[1]=0 and sf[2]>0) or (sf[4]=0 and sf[3]<0) then
   t:=Findindex(f,0);   
   if g[t]=0 then
       ind:=1;
   else
       ind:=0;
   fi;
elif (sg[1]=0 and sg[2]>0) or (sg[4]=0 and sg[3]<0) then
   t:=Findindex(g,0);   
   if f[t]=0 then
       ind:=1;
   else
       ind:=0;
   fi;
else 
   ind:=2;
fi;
return [ind,t];
end:

bisection:=proc(f,g,var,L0,F0,G0)
local F,G,s,L,w,fp,gp,p;
F:=F0;G:=G0;s:=0;L:=L0;
if G[1]*G[2]>0 then
   s:=G[1];
else
   w:=abs(L[2]-L[1]);
   while w>10^(-10) and s=0 do
       p:=(L[1]+L[2])/2;
       fp:=evalf(subs(var=p,f));
       gp:=evalf(subs(var=p,g));
       w:=abs(L[2]-L[1]);
  
       if F[1]*fp<=0 and G[1]*gp<=0 then
          L[2]:=p;
       elif F[2]*fp<=0 and G[2]*gp<=0 then
          L[1]:=p;
       else
          s:=gp;
          
       fi;
   od;
fi;
return s;
end:




Existence_2d:=proc(sys,var,B)
local f,g,uf,ug,p1,p2,p3,p4,f1,f2,f3,f4,g1,g2,g3,g4,s,c,t,ind,P,G,F,L,S;
f:=sys[1];g:=sys[2];
p1:=[B[1][1],B[2][2]];
p2:=[B[1][2],B[2][2]];
p3:=[B[1][2],B[2][1]];
p4:=[B[1][1],B[2][1]];
P:=[p1,p2,p3,p4];
f1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],f));
f2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],f));
f3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],f));
f4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],f));
g1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],g));
g2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],g));
g3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],g));
g4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],g));
F:=[f1,f2,f3,f4];G:=[g1,g2,g3,g4];
c:=whichcase([f1,f2,f3,f4],[g1,g2,g3,g4]);
################################################
if c[1]=0 then
   ind:=0;
elif c[1]=1 then    
   ind:=1/2;
   t:=c[2]; 
else
   S:=[];
   if f1=0 then
      f1:=1;
   fi;
   if f2=0 then
      f2:=1;
   fi;
   if f3=0 then
      f3:=1;
   fi;
   if f4=0 then
      f4:=1;
   fi;
   if f1*f2<0 then
      L:=B[1];
      uf:=subs(var[2]=B[2][2],f);
      ug:=subs(var[2]=B[2][2],g);
      s:=bisection(uf,ug,var[1],L,F[1..2],G[1..2]);
      S:=[op(S),s];
      if s=0 then
         t:=12;
      fi;
   fi;
   if f2*f3<0 then
      L:=[B[2][2],B[2][1]];
      uf:=subs(var[1]=B[1][2],f);
      ug:=subs(var[1]=B[1][2],g);
      s:=bisection(uf,ug,var[2],L,F[2..3],G[2..3]);
      S:=[op(S),s];
      if s=0 then
         t:=23;
      fi;
   fi;
   if f3*f4<0 then
      L:=[B[1][2],B[1][1]];
      uf:=subs(var[2]=B[2][1],f);
      ug:=subs(var[2]=B[2][1],g);
      s:=bisection(uf,ug,var[1],L,F[3..4],G[3..4]);
      S:=[op(S),s];
      if s=0 then
         t:=34;
      fi;
   fi;
   if f4*f1<0 then
      L:=B[2];
      uf:=subs(var[1]=B[1][1],f);
      ug:=subs(var[1]=B[1][1],g);
      s:=bisection(uf,ug,var[2],L,[F[4],F[1]],[G[4],G[1]]);
      S:=[op(S),s];
      if s=0 then
         t:=41;
      fi;
   fi;
   if nops(S)>1 then
      if S[1]*S[2]<0 then
         ind:=1;
      elif S[1]*S[2]>0 then 
         ind:=0;
      else
         ind:=1/2;
      fi;
   else
         ind:=1/2;
         t:=0;
   fi;
fi;
return [ind,t];
end: