axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / SOLVETRA.spad

)abbrev package SOLVETRA TransSolvePackage
++ Author: W. Wiwianka, Martin Rubey
++ Date Created: Summer 1991
++ Date Last Changed: 9/91
++ Description:
++ This package tries to find solutions of equations of type Expression(R).
++ This means expressions involving transcendental, exponential, logarithmic
++ and nthRoot functions.
++ After trying to transform different kernels to one kernel by applying
++ several rules, it calls zerosOf for the SparseUnivariatePolynomial in
++ the remaining kernel.
++ For example the expression  \spad{sin(x)*cos(x)-2}  will be transformed to
++ \spad{-2 tan(x/2)**4 -2 tan(x/2)**3 -4 tan(x/2)**2 +2 tan(x/2) -2}
++ by  using the function  normalize  and then to
++ \spad{-2 tan(x)**2 + tan(x) -2}
++ with help of subsTan. This function tries to express the given function
++ in terms of \spad{tan(x/2)}  to express in terms of \spad{tan(x)} .
++ Other examples are the expressions    \spad{sqrt(x+1)+sqrt(x+7)+1}   or
++ \spad{sqrt(sin(x))+1} .

TransSolvePackage(R) : SIG == CODE where
  R : Join(OrderedSet, EuclideanDomain, RetractableTo Integer, 
            LinearlyExplicitRingOver Integer, CharacteristicZero)

  I   ==> Integer
  NNI ==> NonNegativeInteger
  RE  ==> Expression R
  EQ  ==> Equation
  S   ==> Symbol
  V   ==> Variable
  L   ==> List
  K   ==> Kernel RE
  SUP ==> SparseUnivariatePolynomial
  C   ==> Complex 
  F   ==> Float
  INT ==> Interval
  SMP ==> SparseMultivariatePolynomial


  SIG ==> with

    solve : RE -> L EQ RE
      ++ solve(expr) finds the solutions of the equation expr = 0
      ++ where expr is a function of type Expression(R)
      ++ with respect to the unique symbol x appearing in eq.
 
    solve : EQ RE -> L EQ RE
      ++ solve(eq) finds the solutions of the equation  eq
      ++ where  eq  is  an equation of functions of type Expression(R)
      ++ with respect to the unique symbol x appearing in eq.

    solve : (EQ RE, S) -> L EQ RE
      ++ solve(eq,x) finds the solutions of the equation  eq
      ++ where  eq  is  an equation of functions of type Expression(R)
      ++ with respect to the symbol x.

    solve : (RE, S) ->  L EQ RE
      ++ solve(expr,x) finds the solutions of the equation expr = 0
      ++ with respect to the symbol x where expr is a function
      ++ of type Expression(R).
      ++
      ++X solve(1/2*v*v*cos(theta+phi)*cos(theta+phi)+g*l*cos(phi)=g*l,phi)
      ++X definingPolynomial %phi0
      ++X definingPolynomial %phi1

    solve : (L EQ RE, L S) -> L L EQ RE
      ++ solve(leqs, lvar) returns a list of solutions to the list of 
      ++ equations leqs with respect to the list of symbols lvar.

  CODE ==> add

     import ACF
     import HomogeneousAggregate(R)
     import AlgebraicManipulations(R, RE)
     import TranscendentalManipulations(R, RE)
     import TrigonometricManipulations(R, RE)
     import ElementaryFunctionStructurePackage(R, RE)
     import SparseUnivariatePolynomial(R)
     import LinearSystemMatrixPackage(RE,Vector RE,Vector RE,Matrix RE)
     import TransSolvePackageService(R)
     import MultivariateFactorize(K, IndexedExponents K, R, SMP(R, K))

        ---- Local Function Declarations ----

     solveInner : (RE, S) -> L EQ RE
     tryToTrans : ( RE , S)     ->  RE

     eliminateKernRoot: (RE , K) -> RE
     eliminateRoot: (RE , S) -> RE

     combineLog : ( RE , S ) -> RE
     testLog : ( RE , S ) -> Boolean
     splitExpr : ( RE ) -> L RE
     buildnexpr : ( RE , S ) -> L RE
     logsumtolog : RE -> RE
     logexpp : ( RE , RE ) -> RE

     testRootk : ( RE, S) -> Boolean
     testkernel : ( RE , S ) -> Boolean
     funcinv : ( RE , RE ) -> Union(RE,"failed")
     testTrig : ( RE , S ) -> Boolean
     testHTrig : ( RE , S ) -> Boolean
     tableXkernels : ( RE , S ) -> L RE
     subsTan : ( RE , S )  -> RE
 
     -- exported functions

     solve(oside: RE) : L EQ RE ==
       zero? oside => error "equation is always satisfied"
       lv := variables oside
       empty? lv => error "inconsistent equation"
       #lv>1 => error "too many variables"
       solve(oside,lv.first)

     solve(equ:EQ RE) : L EQ RE ==
       solve(lhs(equ)-rhs(equ))

     solve(equ:EQ RE, x:S) : L EQ RE ==
       oneside:=lhs(equ)-rhs(equ)
       solve(oneside,x)

     testZero?(lside:RE,sol:EQ RE):Boolean ==
       if R has QuotientFieldCategory(Integer) then
         retractIfCan(rhs sol)@Union(Integer,"failed") case "failed" => true
       else
         retractIfCan(rhs sol)@Union(Fraction Integer,"failed") _
            case "failed" => true
       zero? eval(lside,sol) => true
       false

     solve(lside: RE, x:S) : L EQ RE ==
       [sol for sol in solveInner(lside,x) | testZero?(lside,sol)]

     solveInner(lside: RE, x:S) : L EQ RE ==
       lside:=eliminateRoot(lside,x)
       ausgabe1:=tableXkernels(lside,x)

       X:=new()@Symbol
       Y:=new()@Symbol::RE
       (#ausgabe1) = 1 =>
          bigX:= (first ausgabe1)::RE
          eq1:=eval(lside,bigX=(X::RE))
              -- Type  :  Expression R
          f:=univariate(eq1,first kernels (X::RE))
              -- Type  :  Fraction SparseUnivariatePolynomial Expression R
          lfatt:= factors factorPolynomial numer f
          lr:L RE := "append" /[zerosOf(fatt.factor,x) for fatt in lfatt]
              -- Type  :  List Expression R
          r1:=[]::L RE
          for i in 1..#lr repeat
             finv := funcinv(bigX,lr(i))
             if finv case RE then r1:=cons(finv::RE,r1)
          bigX_back:=funcinv(bigX,bigX)::RE
          if not testkernel(bigX_back,x) then
            if bigX = bigX_back then return []::L EQ RE
            return
              "append"/[solve(bigX_back-ri, x) for ri in r1]
          newlist:=[]::L EQ RE

          for i in 1..#r1 repeat
             elR :=  eliminateRoot((numer(bigX_back - r1(i))::RE ),x)
             f:=univariate(elR, kernel(x))
              -- Type  :  Fraction SparseUnivariatePolynomial Expression R
             lfatt:= factors factorPolynomial numer f
             secondsol:="append" /[zerosOf(ff.factor,x) for ff in lfatt]
             for j in 1..#secondsol repeat
                newlist:=cons((x::RE)=rootSimp( secondsol(j) ),newlist)
          newlist
       newlside:=tryToTrans(lside,x) ::RE
       listofkernels:=tableXkernels(newlside,x)
       (#listofkernels) = 1 => solve(newlside,x)
       lfacts := factors factor(numer lside)
       #lfacts > 1 =>
          sols : L EQ RE := []
          for frec in lfacts repeat
              sols := append(solve(frec.factor :: RE, x), sols)
          sols
       return []::L EQ RE

    -- local functions

     --  This function was suggested by Manuel Bronstein as a simpler 
     --  alternative to normalize.
     simplifyingLog(f:RE):RE ==
       (u:=isExpt(f,"exp"::Symbol)) case _
              Record(var:Kernel RE,exponent:Integer) =>       
         rec := u::Record(var:Kernel RE,exponent:Integer)
         rec.exponent * first argument(rec.var)
       log f


     testkernel(var1:RE,y:S) : Boolean ==
       var1:=eliminateRoot(var1,y)
       listvar1:=tableXkernels(var1,y)
       if (#listvar1 = 1) and ((listvar1(1) = (y::RE))@Boolean ) then
            true
       else if #listvar1 = 0 then true
            else false

     solveRetract(lexpr:L RE, lvar:L S):Union(L L EQ RE, "failed") ==
        nlexpr : L Fraction Polynomial R := []
        for expr in lexpr repeat
           rf:Union(Fraction Polynomial R, "failed") := retractIfCan(expr)$RE
           rf case "failed" => return "failed"
           nlexpr := cons(rf, nlexpr)
        radicalSolve(nlexpr, lvar)$RadicalSolvePackage(R)

     tryToTrans(lside: RE, x:S) : RE ==
       if testTrig(lside,x) or testHTrig(lside,x) then
          convLside:=( simplify(lside) )::RE
          resultLside:=convLside
          listConvLside:=tableXkernels(convLside,x)
          if (#listConvLside) > 1  then
            NormConvLside:=normalize(convLside,x)
            NormConvLside:=( NormConvLside ) :: RE
            resultLside:=subsTan(NormConvLside , x)

       else if testLog(lside,x) then
              numlside:=numer(lside)::RE
              resultLside:=combineLog(numlside,x)
            else
              NormConvLside:=normalize(lside,x)
              NormConvLside:=( NormConvLside ) :: RE
              resultLside:=NormConvLside
              listConvLside:=tableXkernels(NormConvLside,x)
              if  (#listConvLside) > 1  then
                cnormConvLside:=complexNormalize(lside,x)
                cnormConvLside:=cnormConvLside::RE
                resultLside:=cnormConvLside
                listcnorm:=tableXkernels(cnormConvLside,x)
                if (#listcnorm) > 1 then
                  if testLog(cnormConvLside,x) then
                    numlside:=numer(cnormConvLside)::RE
                    resultLside:=combineLog(numlside,x)
       resultLside


     subsTan(exprvar:RE,y:S) : RE ==
       Z:=new()@Symbol
       listofkern:=tableXkernels(exprvar,y)
       varkern:=(first listofkern)::RE
       Y:=(numer first argument first (kernels(varkern)))::RE
       test : Boolean := varkern=tan(((Y::RE)/(2::RE))::RE)
       if not( (#listofkern=1) and test) then
         return exprvar
       fZ:=eval(exprvar,varkern=(Z::RE))
       fN:=(numer fZ)::RE
       f:=univariate(fN, first kernels(Z::RE))
       secondfun:=(-2*(Y::RE)/((Y::RE)**2-1) )::RE
       g:=univariate(secondfun,first kernels(y::RE))
       H:=(new()@Symbol)::RE
       newH:=univariate(H,first kernels(Z::RE))
       result:=decomposeFunc(f,g,newH)
       if not ( result = f ) then
         result1:=result( H::RE )
         resultnew:=eval(result1,H=(( tan((Y::RE))::RE ) ))
       else return exprvar


     eliminateKernRoot(var: RE, varkern: K) : RE ==
       X:=new()@Symbol
       var1:=eval(var, (varkern::RE)=(X::RE) )
       var2:=numer univariate(var1, first kernels(X::RE))
       var3:= monomial(1, ( retract( second argument varkern)@I )::NNI)@SUP RE_
              - monomial(first argument varkern, 0::NNI)@SUP RE
       resultvar:=resultant(var2, var3)

     eliminateRoot(var:RE, y:S) : RE ==
       var1:=var
       while testRootk(var1,y) repeat
         varlistk1:=tableXkernels(var1,y)
         for i in varlistk1 repeat
            if is?(i, "nthRoot"::S) then
              var1:=eliminateKernRoot(var1,first kernels(i::RE))
       var1


     logsumtolog(var:RE) : RE ==
       (listofexpr:=isPlus(var)) case "failed" => var
       listofexpr:= listofexpr ::L RE
       listforgcd:=[]::L R
       for i in listofexpr repeat
          exprcoeff:=leadingCoefficient(numer(i))
          listforgcd:=cons(exprcoeff, listforgcd)
       gcdcoeff:=gcd(listforgcd)::RE
       newexpr:RE :=0
       for i in listofexpr repeat
          exprlist:=splitExpr(i::RE)
          newexpr:=newexpr + logexpp(exprlist.2, exprlist.1/gcdcoeff)
       kernelofvar:=kernels(newexpr)
       var2:=1::RE
       for i in kernelofvar repeat
          var2:=var2*(first argument i)
       gcdcoeff * log(var2)


     testLog(expr:RE,Z:S) : Boolean ==
       testList:=[log]::L S
       kernelofexpr:=tableXkernels(expr,Z)
       if #kernelofexpr = 0 then
         return false
       for i in kernelofexpr repeat
          if not member?(name(first kernels(i)),testList) or _
             not testkernel( (first argument first kernels(i)) ,Z) then
            return false
       true

     splitExpr(expr:RE) : L RE ==
       lcoeff:=leadingCoefficient((numer expr))
       exprwcoeff:=expr
       listexpr:=isTimes(exprwcoeff)
       if listexpr case "failed" then
         [1::RE , expr]
       else
         listexpr:=remove_!(lcoeff::RE , listexpr)
         cons(lcoeff::RE , listexpr)

     buildnexpr(expr:RE, Z:S) : L RE ==
       nlist:=splitExpr(expr)
       n2list:=remove_!(nlist.1, nlist)
       anscoeff:RE:=1
       ansmant:RE:=0
       for i in n2list repeat
          if freeOf?(i::RE,Z) then
            anscoeff:=(i::RE)*anscoeff
          else
            ansmant:=(i::RE)
       [anscoeff, ansmant * nlist.1 ]

     logexpp(expr1:RE, expr2:RE) : RE ==
       log( (first argument first kernels(expr1))**expr2 )

     combineLog(expr:RE,Y:S) : RE ==
       exprtable:Table(RE,RE):=table()
       (isPlus(expr)) case "failed" => expr
       ans:RE:=0
       while expr ^= 0 repeat
         loopexpr:RE:=leadingMonomial(numer(expr))::RE
         if testLog(loopexpr,Y) and (#tableXkernels(loopexpr,Y)=1) then
           exprr:=buildnexpr(loopexpr,Y)
           if search(exprr.1,exprtable) case "failed" then
             exprtable.(exprr.1):=0
           exprtable.(exprr.1):= exprtable.(exprr.1) + exprr.2
         else
           ans:=ans+loopexpr
         expr:=(reductum(numer expr))::RE
       ansexpr:RE:=0
       for i in keys(exprtable) repeat
          ansexpr:=ansexpr + logsumtolog(exprtable.i) * (i::RE)
       ansexpr:=ansexpr + ans


     testRootk(varlistk:RE,y:S) : Boolean ==
       testList:=[nthRoot]::L S
       kernelofeqnvar:=tableXkernels(varlistk,y)
       if #kernelofeqnvar = 0 then
         return false
       for i in kernelofeqnvar repeat
          if member?(name(first kernels(i)),testList) then
            return true
       false

     tableXkernels(evar:RE,Z:S) : L RE ==
       kOfvar:=kernels(evar)
       listkOfvar:=[]::L RE
       for i in kOfvar repeat
          if not freeOf?(i::RE,Z) then
              listkOfvar:=cons(i::RE,listkOfvar)
       listkOfvar

     testTrig(eqnvar:RE,Z:S) : Boolean ==
       testList:=[sin , cos , tan , cot , sec , csc]::L S
       kernelofeqnvar:=tableXkernels(eqnvar,Z)
       if #kernelofeqnvar = 0 then
         return false
       for i in kernelofeqnvar repeat
          if not member?(name(first kernels(i)),testList) or _
             not testkernel( (first argument first kernels(i)) ,Z) then
            return false
       true


     testHTrig(eqnvar:RE,Z:S) : Boolean ==
       testList:=[sinh , cosh , tanh , coth , sech , csch]::L S
       kernelofeqnvar:=tableXkernels(eqnvar,Z)
       if #kernelofeqnvar = 0 then
         return false
       for i in kernelofeqnvar repeat
          if not member?(name(first kernels(i)),testList) or _
             not testkernel( (first argument first kernels(i)) ,Z) then
            return false
       true

     -- Auxiliary local function for use in funcinv.
     makeInterval(l:R):C INT F ==
       if R has complex and R has ConvertibleTo(C F) then
         map(interval$INT(F),convert(l)$R)$ComplexFunctions2(F,INT F)
       else
         error "This should never happen"

     funcinv(k:RE,l:RE) : Union(RE,"failed") ==
       is?(k, "sin"::Symbol)   => asin(l)
       is?(k, "cos"::Symbol)   => acos(l)
       is?(k, "tan"::Symbol)   => atan(l)
       is?(k, "cot"::Symbol)   => acot(l)
       is?(k, "sec"::Symbol)   =>
           l = 0 => "failed"
           asec(l)
       is?(k, "csc"::Symbol)   =>
           l = 0 => "failed"
           acsc(l)
       is?(k, "sinh"::Symbol)  => asinh(l)
       is?(k, "cosh"::Symbol)  => acosh(l)
       is?(k, "tanh"::Symbol)  => atanh(l)
       is?(k, "coth"::Symbol)  => acoth(l)
       is?(k, "sech"::Symbol)  => asech(l)
       is?(k, "csch"::Symbol)  => acsch(l)
       is?(k, "atan"::Symbol)  => tan(l)
       is?(k, "acot"::Symbol)  =>
           l = 0 => "failed"
           cot(l)
       is?(k, "asin"::Symbol)  => sin(l)
       is?(k, "acos"::Symbol)  => cos(l)
       is?(k, "asec"::Symbol)  => sec(l)
       is?(k, "acsc"::Symbol)  =>
           l = 0 => "failed"
           csc(l)
       is?(k, "asinh"::Symbol) => sinh(l)
       is?(k, "acosh"::Symbol) => cosh(l)
       is?(k, "atanh"::Symbol) => tanh(l)
       is?(k, "acoth"::Symbol) =>
           l = 0 => "failed"
           coth(l)
       is?(k, "asech"::Symbol) => sech(l)
       is?(k, "acsch"::Symbol) =>
           l = 0 => "failed"
           csch(l)
       is?(k, "exp"::Symbol)   =>
           l = 0 => "failed"
           simplifyingLog l
       is?(k, "log"::Symbol)   =>
         if R has complex and R has ConvertibleTo(C F) then
           -- We will check to see if the imaginary part lies in [-Pi,Pi)
           ze : Expression C INT F
           ze := map(makeInterval,l)$ExpressionFunctions2(R,C INT F)
           z : Union(C INT F,"failed") := retractIfCan ze
           z case "failed" => exp l
           im := imag z
           fpi : Float := pi()
           (-fpi < inf(im)) and (sup(im) <= fpi) => exp l
           "failed"
         else -- R not Complex or something which doesn't map to Complex Floats
           exp l 
       is?(k, "%power"::Symbol)   => 
            (t:=normalize(l)) = 0 => "failed"
            log t
       l

     import SystemSolvePackage(RE)

     ker2Poly(k:Kernel RE, lvar:L S):Polynomial RE ==
        member?(nm:=name k, lvar) => nm :: Polynomial RE
        k :: RE :: Polynomial RE

     smp2Poly(pol:SMP(R,Kernel RE), lvar:L S):Polynomial RE ==
        map(x +-> ker2Poly(x, lvar),
            y +-> y::RE::Polynomial RE, pol)$PolynomialCategoryLifting(
              IndexedExponents Kernel RE, Kernel RE, R, SMP(R, Kernel RE), 
                      Polynomial RE)

     makeFracPoly(expr:RE, lvar:L S):Fraction Polynomial RE ==
        smp2Poly(numer expr, lvar) / smp2Poly(denom expr, lvar)

     makeREpol(pol:Polynomial RE):RE ==
        lvar := variables pol
        lval : List RE := [v::RE for v in lvar]
        ground eval(pol,lvar,lval)

     makeRE(frac:Fraction Polynomial RE):RE ==
        makeREpol(numer frac)/makeREpol(denom frac)

     solve1Pol(pol:Polynomial RE, var: S, sol:L EQ RE):L L EQ RE ==
        repol := eval(makeREpol pol, sol)
        vsols := solve(repol, var)
        [cons(vsol, sol) for vsol in vsols]

     solve1Sys(plist:L Polynomial RE, lvar:L S):L L EQ RE ==
        rplist := reverse plist
        rlvar := reverse lvar
        sols : L L EQ RE := list(empty())
        for p in rplist for v in rlvar repeat
           sols := "append"/[solve1Pol(p,v,sol) for sol in sols]
        sols


     solveList(lexpr:L RE, lvar:L S):L L EQ RE ==
        ans1 := solveRetract(lexpr, lvar)
        not(ans1 case "failed") => ans1 :: L L EQ RE
        lfrac:L Fraction Polynomial RE :=
           [makeFracPoly(expr, lvar) for expr in lexpr]
        trianglist := triangularSystems(lfrac, lvar)
        l: L L L EQ RE := [solve1Sys(plist, lvar) for plist in trianglist]
        reduce(append, l, [])
        
     solve(leqs:L EQ RE, lvar:L S):L L EQ RE ==
        lexpr:L RE := [lhs(eq)-rhs(eq) for eq in leqs]
        solveList(lexpr, lvar)