--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2009, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

-- SPAD files for the functional world should be compiled in the
-- following order:
--
--   op  kl  fspace ALGFUNC expr

)abbrev category ACF AlgebraicallyClosedField
++ Author: Manuel Bronstein
++ Date Created: 22 Mar 1988
++ Date Last Updated: 27 November 1991
++ Description:
++   Model for algebraically closed fields.
++ Keywords: algebraic, closure, field.

AlgebraicallyClosedField(): Category == Join(Field,RadicalCategory) with
    rootOf: Polynomial $ -> $
      ++ rootOf(p) returns y such that \spad{p(y) = 0}.
      ++ Error: if p has more than one variable y.
    rootOf: SparseUnivariatePolynomial $ -> $
      ++ rootOf(p) returns y such that \spad{p(y) = 0}.
    rootOf: (SparseUnivariatePolynomial $, Symbol) -> $
      ++ rootOf(p, y) returns y such that \spad{p(y) = 0}.
      ++ The object returned displays as \spad{'y}.
    rootsOf: Polynomial $ -> List $
      ++ rootsOf(p) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ Note: the returned symbols y1,...,yn are bound in the
      ++ interpreter to respective root values.
      ++ Error: if p has more than one variable y.
    rootsOf: SparseUnivariatePolynomial $ -> List $
      ++ rootsOf(p) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ Note: the returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
    rootsOf: (SparseUnivariatePolynomial $, Symbol) -> List $
      ++ rootsOf(p, y) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0};
      ++ The returned roots display as \spad{'y1},...,\spad{'yn}.
      ++ Note: the returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
    zeroOf: Polynomial $ -> $
      ++ zeroOf(p) returns y such that \spad{p(y) = 0}.
      ++ If possible, y is expressed in terms of radicals.
      ++ Otherwise it is an implicit algebraic quantity.
      ++ Error: if p has more than one variable y.
    zeroOf: SparseUnivariatePolynomial $ -> $
      ++ zeroOf(p) returns y such that \spad{p(y) = 0};
      ++ if possible, y is expressed in terms of radicals.
      ++ Otherwise it is an implicit algebraic quantity.
    zeroOf: (SparseUnivariatePolynomial $, Symbol) -> $
      ++ zeroOf(p, y) returns y such that \spad{p(y) = 0};
      ++ if possible, y is expressed in terms of radicals.
      ++ Otherwise it is an implicit algebraic quantity which
      ++ displays as \spad{'y}.
    zerosOf: Polynomial $ -> List $
      ++ zerosOf(p) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ The yi's are expressed in radicals if possible.
      ++ Otherwise they are implicit algebraic quantities.
      ++ The returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
      ++ Error: if p has more than one variable y.
    zerosOf: SparseUnivariatePolynomial $ -> List $
      ++ zerosOf(p) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ The yi's are expressed in radicals if possible, and otherwise
      ++ as implicit algebraic quantities.
      ++ The returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
    zerosOf: (SparseUnivariatePolynomial $, Symbol) -> List $
      ++ zerosOf(p, y) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ The yi's are expressed in radicals if possible, and otherwise
      ++ as implicit algebraic quantities
      ++ which display as \spad{'yi}.
      ++ The returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
 add
    SUP ==> SparseUnivariatePolynomial $

    assign  : (Symbol, $) -> $
    allroots: (SUP, Symbol, (SUP, Symbol) -> $) -> List $
    binomialRoots: (SUP, Symbol, (SUP, Symbol) -> $) -> List $

    zeroOf(p:SUP)            == assign(x := new(), zeroOf(p, x))
    rootOf(p:SUP)            == assign(x := new(), rootOf(p, x))
    zerosOf(p:SUP)           == zerosOf(p, new())
    rootsOf(p:SUP)           == rootsOf(p, new())
    rootsOf(p:SUP, y:Symbol) == allroots(p, y, rootOf)
    zerosOf(p:SUP, y:Symbol) == allroots(p, y, zeroOf)
    assign(x, f)             ==
      assignSymbol(x, f, $)$Foreign(Builtin)
      f

    zeroOf(p:Polynomial $) ==
      empty?(l := variables p) => error "zeroOf: constant polynomial"
      zeroOf(univariate p, first l)

    rootOf(p:Polynomial $) ==
      empty?(l := variables p) => error "rootOf: constant polynomial"
      rootOf(univariate p, first l)

    zerosOf(p:Polynomial $) ==
      empty?(l := variables p) => error "zerosOf: constant polynomial"
      zerosOf(univariate p, first l)

    rootsOf(p:Polynomial $) ==
      empty?(l := variables p) => error "rootsOf: constant polynomial"
      rootsOf(univariate p, first l)

    zeroOf(p:SUP, y:Symbol) ==
      zero?(d := degree p) => error "zeroOf: constant polynomial"
      zero? coefficient(p, 0) => 0
      a := leadingCoefficient p
      d = 2 =>
        b := coefficient(p, 1)
        (sqrt(b**2 - 4 * a * coefficient(p, 0)) - b) / (2 * a)
      (r := retractIfCan(reductum p)@Union($,"failed")) case "failed" =>
        rootOf(p, y)
      nthRoot(- (r::$ / a), d)

    binomialRoots(p, y, fn) ==
     -- p = a * x**n + b
      alpha := assign(x := new(y)$Symbol, fn(p, x))
      one?(n := degree p) =>  [ alpha ]
      cyclo := cyclotomic(n, monomial(1,1)$SUP)$NumberTheoreticPolynomialFunctions(SUP)
      beta := assign(x := new(y)$Symbol, fn(cyclo, x))
      [alpha*beta**i for i in 0..(n-1)::NonNegativeInteger]

    import PolynomialDecomposition(SUP,$)

    allroots(p, y, fn) ==
      zero? p => error "allroots: polynomial must be nonzero"
      zero? coefficient(p,0) =>
         concat(0, allroots(p quo monomial(1,1), y, fn))
      zero?(p1:=reductum p) => empty()
      zero? reductum p1 => binomialRoots(p, y, fn)
      decompList := decompose(p)
      # decompList > 1 =>
          h := last decompList
          g := leftFactor(p,h) :: SUP
          groots := allroots(g, y, fn)
          "append"/[allroots(h-r::SUP, y, fn) for r in groots]
      ans := nil()$List($)
      while not ground? p repeat
        alpha := assign(x := new(y)$Symbol, fn(p, x))
        q     := monomial(1, 1)$SUP - alpha::SUP
        if not zero?(p alpha) then
          p   := p quo q
          ans := concat(alpha, ans)
        else while zero?(p alpha) repeat
          p   := (p exquo q)::SUP
          ans := concat(alpha, ans)
      reverse! ans

)abbrev category ACFS AlgebraicallyClosedFunctionSpace
++ Author: Manuel Bronstein
++ Date Created: 31 October 1988
++ Date Last Updated: 7 October 1991
++ Description:
++   Model for algebraically closed function spaces.
++ Keywords: algebraic, closure, field.
AlgebraicallyClosedFunctionSpace(R: IntegralDomain):
 Category == Join(AlgebraicallyClosedField, FunctionSpace R) with
    rootOf : $ -> $
      ++ rootOf(p) returns y such that \spad{p(y) = 0}.
      ++ Error: if p has more than one variable y.
    rootsOf: $ -> List $
      ++ rootsOf(p, y) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0};
      ++ Note: the returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
      ++ Error: if p has more than one variable y.
    rootOf : ($, Symbol) -> $
      ++ rootOf(p,y) returns y such that \spad{p(y) = 0}.
      ++ The object returned displays as \spad{'y}.
    rootsOf: ($, Symbol) -> List $
      ++ rootsOf(p, y) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0};
      ++ The returned roots display as \spad{'y1},...,\spad{'yn}.
      ++ Note: the returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
    zeroOf : $ -> $
      ++ zeroOf(p) returns y such that \spad{p(y) = 0}.
      ++ The value y is expressed in terms of radicals if possible,and otherwise
      ++ as an implicit algebraic quantity.
      ++ Error: if p has more than one variable.
    zerosOf: $ -> List $
      ++ zerosOf(p) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ The yi's are expressed in radicals if possible.
      ++ The returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
      ++ Error: if p has more than one variable.
    zeroOf : ($, Symbol) -> $
      ++ zeroOf(p, y) returns y such that \spad{p(y) = 0}.
      ++ The value y is expressed in terms of radicals if possible,and otherwise
      ++ as an implicit algebraic quantity
      ++ which displays as \spad{'y}.
    zerosOf: ($, Symbol) -> List $
      ++ zerosOf(p, y) returns \spad{[y1,...,yn]} such that \spad{p(yi) = 0}.
      ++ The yi's are expressed in radicals if possible, and otherwise
      ++ as implicit algebraic quantities
      ++ which display as \spad{'yi}.
      ++ The returned symbols y1,...,yn are bound in the interpreter
      ++ to respective root values.
  add
    rootOf(p:$) ==
      empty?(l := variables p) => error "rootOf: constant expression"
      rootOf(p, first l)

    rootsOf(p:$) ==
      empty?(l := variables p) => error "rootsOf: constant expression"
      rootsOf(p, first l)

    zeroOf(p:$) ==
      empty?(l := variables p) => error "zeroOf: constant expression"
      zeroOf(p, first l)

    zerosOf(p:$) ==
      empty?(l := variables p) => error "zerosOf: constant expression"
      zerosOf(p, first l)

    zeroOf(p:$, x:Symbol) ==
      n := numer(f := univariate(p, kernel(x)$Kernel($)))
      positive? degree denom f => error "zeroOf: variable appears in denom"
      degree n = 0 => error "zeroOf: constant expression"
      zeroOf(n, x)

    rootOf(p:$, x:Symbol) ==
      n := numer(f := univariate(p, kernel(x)$Kernel($)))
      positive? degree denom f => error "roofOf: variable appears in denom"
      degree n = 0 => error "rootOf: constant expression"
      rootOf(n, x)

    zerosOf(p:$, x:Symbol) ==
      n := numer(f := univariate(p, kernel(x)$Kernel($)))
      positive? degree denom f => error "zerosOf: variable appears in denom"
      degree n = 0 => empty()
      zerosOf(n, x)

    rootsOf(p:$, x:Symbol) ==
      n := numer(f := univariate(p, kernel(x)$Kernel($)))
      positive? degree denom f => error "roofsOf: variable appears in denom"
      degree n = 0 => empty()
      rootsOf(n, x)

    rootsOf(p:SparseUnivariatePolynomial $, y:Symbol) ==
      (r := retractIfCan(p)@Union($,"failed")) case $ => rootsOf(r::$,y)
      rootsOf(p, y)$AlgebraicallyClosedField_&($)

    zerosOf(p:SparseUnivariatePolynomial $, y:Symbol) ==
      (r := retractIfCan(p)@Union($,"failed")) case $ => zerosOf(r::$,y)
      zerosOf(p, y)$AlgebraicallyClosedField_&($)

    zeroOf(p:SparseUnivariatePolynomial $, y:Symbol) ==
      (r := retractIfCan(p)@Union($,"failed")) case $ => zeroOf(r::$, y)
      zeroOf(p, y)$AlgebraicallyClosedField_&($)

)abbrev package AF AlgebraicFunction
++ Author: Manuel Bronstein
++ Date Created: 21 March 1988
++ Date Last Updated: 11 November 1993
++ Description:
++   This package provides algebraic functions over an integral domain.
++ Keywords: algebraic, function.

AlgebraicFunction(R, F): Exports == Implementation where
  R: IntegralDomain
  F: FunctionSpace R

  SE  ==> Symbol
  Z   ==> Integer
  Q   ==> Fraction Z
  OP  ==> BasicOperator
  K   ==> Kernel F
  P   ==> SparseMultivariatePolynomial(R, K)
  UP  ==> SparseUnivariatePolynomial F
  UPR ==> SparseUnivariatePolynomial R

  Exports ==> with
    rootOf  : (UP, SE) -> F
      ++ rootOf(p, y) returns y such that \spad{p(y) = 0}.
      ++ The object returned displays as \spad{'y}.
    operator: OP -> OP
      ++ operator(op) returns a copy of \spad{op} with the domain-dependent
      ++ properties appropriate for \spad{F}.
      ++ Error: if op is not an algebraic operator, that is,
      ++ an nth root or implicit algebraic operator.
    belong? : OP -> Boolean
      ++ belong?(op) is true if \spad{op} is an algebraic operator, that is,
      ++ an nth root or implicit algebraic operator.
    inrootof: (UP, F) -> F
      ++ inrootof(p, x) should be a non-exported function.
      -- un-export when the compiler accepts conditional local functions!
    droot : List F -> OutputForm
      ++ droot(l) should be a non-exported function.
      -- un-export when the compiler accepts conditional local functions!
    if R has RetractableTo Integer then
      **   : (F, Q) -> F
        ++ x ** q is \spad{x} raised to the rational power \spad{q}.
      minPoly: K  -> UP
        ++ minPoly(k) returns the defining polynomial of \spad{k}.
      definingPolynomial: F -> F
        ++ definingPolynomial(f) returns the defining polynomial of \spad{f}
        ++ as an element of \spad{F}.
        ++ Error: if f is not a kernel.
      iroot : (R, Z) -> F
        ++ iroot(p, n) should be a non-exported function.
        -- un-export when the compiler accepts conditional local functions!

  Implementation ==> add
    macro ALGOP == '%alg
    macro SPECIALDISP == '%specialDisp
    macro SPECIALDIFF == '%specialDiff
    ialg : List F -> F
    dvalg: (List F, SE) -> F
    dalg : List F -> OutputForm

    opalg  := operator('rootOf)$CommonOperators
    oproot := operator('nthRoot)$CommonOperators

    belong? op == has?(op, ALGOP)
    dalg l     == second(l)::OutputForm

    rootOf(p, x) ==
      k := kernel(x)$K
      (r := retractIfCan(p)@Union(F, "failed")) case "failed" =>
        inrootof(p, k::F)
      n := numer(f := univariate(r::F, k))
      positive? degree denom f => error "roofOf: variable appears in denom"
      inrootof(n, k::F)

    dvalg(l, x) ==
      p := numer univariate(first l, retract(second l)@K)
      alpha := kernel(opalg, l)
      - (map(differentiate(#1, x), p) alpha) / ((differentiate p) alpha)

    ialg l ==
      f := univariate(p := first l, retract(x := second l)@K)
      positive? degree denom f => error "roofOf: variable appears in denom"
      inrootof(numer f, x)

    operator op ==
      is?(op,'rootOf) => opalg
      is?(op,'nthRoot) => oproot
      error "Unknown operator"

    if R has AlgebraicallyClosedField then
      UP2R: UP -> Union(UPR, "failed")

      inrootof(q, x) ==
        monomial? q => 0

        (d := degree q) <= 0 => error "rootOf: constant polynomial"
        one? d=> - leadingCoefficient(reductum q) / leadingCoefficient q
        ((rx := retractIfCan(x)@Union(SE, "failed")) case SE) and
          ((r := UP2R q) case UPR) => rootOf(r::UPR, rx::SE)::F
        kernel(opalg, [q x, x])

      UP2R p ==
        ans:UPR := 0
        while p ~= 0 repeat
          (r := retractIfCan(leadingCoefficient p)@Union(R, "failed"))
            case "failed" => return "failed"
          ans := ans + monomial(r::R, degree p)
          p   := reductum p
        ans

    else
      inrootof(q, x) ==
        monomial? q => 0
        (d := degree q) <= 0 => error "rootOf: constant polynomial"
        one? d => - leadingCoefficient(reductum q) /leadingCoefficient q
        kernel(opalg, [q x, x])

    evaluate(opalg, ialg)$BasicOperatorFunctions1(F)
    setProperty(opalg, SPECIALDIFF,
                              dvalg@((List F, SE) -> F) pretend None)
    setProperty(opalg, SPECIALDISP,
                              dalg@(List F -> OutputForm) pretend None)

    if R has RetractableTo Integer then
      import PolynomialRoots(IndexedExponents K, K, R, P, F)

      dumvar := '%%var::F

      lzero   : List F -> F
      dvroot  : List F -> F
      inroot  : List F -> F
      hackroot: (F, Z) -> F
      inroot0 : (F, Z, Boolean, Boolean) -> F

      lzero l == 0

      droot l ==
        x := first(l)::OutputForm
        (n := retract(second l)@Z) = 2 => root x
        root(x, n::OutputForm)

      dvroot l ==
        n := retract(second l)@Z
        (first(l) ** ((1 - n) / n)) / (n::F)

      x ** q ==
        qr := divide(numer q, denom q)
        x ** qr.quotient * inroot([x, (denom q)::F]) ** qr.remainder

      hackroot(x, n) ==
        (n = 1) or (x = 1) => x
        (not one?(dx := denom x) and
           ((rx := retractIfCan(dx)@Union(Integer,"failed")) case Integer) and
           positive?(rx))
           => hackroot((numer x)::F, n)/hackroot(rx::Integer::F, n)
        (x = -1) and n = 4 =>
          ((-1::F) ** (1::Q / 2::Q) + 1) / ((2::F) ** (1::Q / 2::Q))
        kernel(oproot, [x, n::F])


      inroot l ==
        zero?(n := retract(second l)@Z) => error "root: exponent = 0"
        one?(x := first l) or one? n => x
        (r := retractIfCan(x)@Union(R,"failed")) case R => iroot(r::R,n)
        (u := isExpt(x)) case Record(var:K, exponent:Z) =>
          pr := u::Record(var:K, exponent:Z)
          is?(pr.var,oproot) and #argument(pr.var) = 2 =>
	    (first argument(pr.var)) **
		(pr.exponent /$Fraction(Z)
		     (n * retract(second argument(pr.var))@Z))
          inroot0(x, n, false, false)
        inroot0(x, n, false, false)

-- removes powers of positive integers from numer and denom
-- num? or den? is true if numer or denom already processed
      inroot0(x, n, num?, den?) ==
        rn:Union(Z, "failed") := (num? => "failed"; retractIfCan numer x)
        rd:Union(Z, "failed") := (den? => "failed"; retractIfCan denom x)
        (rn case Z) and (rd case Z) =>
          rec := qroot(rn::Z / rd::Z, n::NonNegativeInteger)
          rec.coef * hackroot(rec.radicand, rec.exponent)
        rn case Z =>
          rec := qroot(rn::Z::Fraction(Z), n::NonNegativeInteger)
          rec.coef * inroot0((rec.radicand**(n exquo rec.exponent)::Z)
                                / (denom(x)::F), n, true, den?)
        rd case Z =>
          rec := qroot(rd::Z::Fraction(Z), n::NonNegativeInteger)
          inroot0((numer(x)::F) /
                  (rec.radicand ** (n exquo rec.exponent)::Z),
                   n, num?, true) / rec.coef
        hackroot(x, n)

      if R has AlgebraicallyClosedField then iroot(r, n) == nthRoot(r, n)::F
      else
        iroot0: (R, Z) -> F

        if R has RadicalCategory then
          if R has imaginary:() -> R then iroot(r, n) == nthRoot(r, n)::F
          else
            iroot(r, n) ==
              odd? n or not before?(r,0) => nthRoot(r, n)::F
              iroot0(r, n)

        else iroot(r, n) == iroot0(r, n)

        iroot0(r, n) ==
          rec := rroot(r, n::NonNegativeInteger)
          rec.coef * hackroot(rec.radicand, rec.exponent)

      definingPolynomial x ==
        (r := retractIfCan(x)@Union(K, "failed")) case K =>
          is?(k := r::K, opalg) => first argument k
          is?(k, oproot) =>
            dumvar ** retract(second argument k)@Z - first argument k
          dumvar - x
        dumvar - x

      minPoly k ==
        is?(k, opalg)  =>
           numer univariate(first argument k,
                                           retract(second argument k)@K)
        is?(k, oproot) =>
           monomial(1,retract(second argument k)@Z :: NonNegativeInteger)
             - first(argument k)::UP
        monomial(1, 1) - k::F::UP

      evaluate(oproot, inroot)$BasicOperatorFunctions1(F)
      derivative(oproot, [dvroot, lzero])

    else   -- R is not retractable to Integer
      droot l ==
        x := first(l)::OutputForm
        (n := second l) = 2::F => root x
        root(x, n::OutputForm)

      minPoly k ==
        is?(k, opalg)  =>
           numer univariate(first argument k,
                                           retract(second argument k)@K)
        monomial(1, 1) - k::F::UP

    setProperty(oproot, SPECIALDISP,
                              droot@(List F -> OutputForm) pretend None)

