--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--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.

)abbrev domain IAN InnerAlgebraicNumber
++ Algebraic closure of the rational numbers
++ Author: Manuel Bronstein
++ Date Created: 22 March 1988
++ Date Last Updated: 4 October 1995 (JHD)
++ Description: Algebraic closure of the rational numbers.
++ Keywords: algebraic, number.
InnerAlgebraicNumber(): Exports == Implementation where
  Z   ==> Integer
  FE  ==> Expression Z
  K   ==> Kernel %
  P   ==> SparseMultivariatePolynomial(Z, K)
  SUP ==>  SparseUnivariatePolynomial

  Exports ==> Join(ExpressionSpace, AlgebraicallyClosedField,
                   RetractableTo Z, RetractableTo Fraction Z,
                   LinearlyExplicitRingOver Z, RealConstant,
                   LinearlyExplicitRingOver Fraction Z,
                   CharacteristicZero,
                   ConvertibleTo Complex Float, DifferentialRing,
                   CoercibleFrom P) with
    numer  : % -> P
      ++ numer(f) returns the numerator of f viewed as a
      ++ polynomial in the kernels over Z.
    denom  : % -> P
      ++ denom(f) returns the denominator of f viewed as a
      ++ polynomial in the kernels over Z.
    reduce : % -> %
      ++ reduce(f) simplifies all the unreduced algebraic numbers
      ++ present in f by applying their defining relations.
    trueEqual : (%,%) -> Boolean
      ++ trueEqual(x,y) tries to determine if the two numbers are equal
    norm : (SUP(%),Kernel %) -> SUP(%)
      ++ norm(p,k) computes the norm of the polynomial p
      ++ with respect to the extension generated by kernel k
    norm : (SUP(%),List Kernel %) -> SUP(%)
      ++ norm(p,l) computes the norm of the polynomial p
      ++ with respect to the extension generated by kernels l
    norm : (%,Kernel %) -> %
      ++ norm(f,k) computes the norm of the algebraic number f
      ++ with respect to the extension generated by kernel k
    norm : (%,List Kernel %) -> %
      ++ norm(f,l) computes the norm of the algebraic number f
      ++ with respect to the extension generated by kernels l
  Implementation ==> FE add
    macro ALGOP == '%alg

    Rep := FE

    -- private
    mainRatDenom(f:%):% ==
       ratDenom(f::Rep::FE)$AlgebraicManipulations(Integer, FE)::Rep::%
--        mv:= mainVariable denom f
--        mv case "failed" => f
--        algv:=mv::K
--	q:=univariate(f, algv, minPoly(algv))$PolynomialCategoryQuotientFunctions(IndexedExponents K,K,Integer,P,%)
--	q(algv::%)

    findDenominator(z:SUP %):Record(num:SUP %,den:%) ==
       zz:=z
       while not(zz=0) repeat
          dd:=(denom leadingCoefficient zz)::%
          not(dd=1) =>
             rec:=findDenominator(dd*z)
             return [rec.num,rec.den*dd]
          zz:=reductum zz
       [z,1]
    makeUnivariate(p:P,k:Kernel %):SUP % ==
      map(#1::%,univariate(p,k))$SparseUnivariatePolynomialFunctions2(P,%)
    -- public
    a,b:%
    differentiate(x:%):% == 0
    zero? a == zero? numer a
    one? a == one? numer a and one? denom a
    x:% / y:%        == mainRatDenom(x /$Rep y)
    x:% ** n:Integer ==
      negative? n => mainRatDenom (x **$Rep n)
      x **$Rep n
    trueEqual(a,b) ==
       -- if two algebraic numbers have the same norm (after deleting repeated
       -- roots, then they are certainly conjugates. Note that we start with a
       -- monic polynomial, so don't have to check for constant factors.
       -- this will be fooled by sqrt(2) and -sqrt(2), but the = in
       -- AlgebraicNumber knows what to do about this.
       ka:=reverse tower a
       kb:=reverse tower b
       empty? ka and empty? kb => retract(a)@Fraction Z = retract(b)@Fraction Z
       pa,pb:SparseUnivariatePolynomial %
       pa:=monomial(1,1)-monomial(a,0)
       pb:=monomial(1,1)-monomial(b,0)
       na:=map(retract,norm(pa,ka))$SparseUnivariatePolynomialFunctions2(%,Fraction Z)
       nb:=map(retract,norm(pb,kb))$SparseUnivariatePolynomialFunctions2(%,Fraction Z)
       (sa:=squareFreePart(na)) = (sb:=squareFreePart(nb)) => true
       g:=gcd(sa,sb)
       (dg:=degree g) = 0 => false
       -- of course, if these have a factor in common, then the
       -- answer is really ambiguous, so we ought to be using Duval-type
       -- technology
       dg = degree sa or dg = degree sb => true
       false
    norm(z:%,k:Kernel %): % ==
       p:=minPoly k
       n:=makeUnivariate(numer z,k)
       d:=makeUnivariate(denom z,k)
       resultant(n,p)/resultant(d,p)
    norm(z:%,l:List Kernel %): % ==
       for k in l repeat
           z:=norm(z,k)
       z
    norm(z:SUP %,k:Kernel %):SUP % ==
       p:=map(#1::SUP %,minPoly k)$SparseUnivariatePolynomialFunctions2(%,SUP %)
       f:=findDenominator z
       zz:=map(makeUnivariate(numer #1,k),f.num)$SparseUnivariatePolynomialFunctions2( %,SUP %)
       zz:=swap(zz)$CommuteUnivariatePolynomialCategory(%,SUP %,SUP SUP %)
       resultant(p,zz)/norm(f.den,k)
    norm(z:SUP %,l:List Kernel %): SUP % ==
       for k in l repeat
           z:=norm(z,k)
       z
    belong? op           == belong?(op)$ExpressionSpace_&(%) or has?(op, ALGOP)

    convert(x:%):Float ==
      retract map(#1::Float, x pretend FE)$ExpressionFunctions2(Z,Float)

    convert(x:%):DoubleFloat ==
      retract map(#1::DoubleFloat,
                  x pretend FE)$ExpressionFunctions2(Z, DoubleFloat)

    convert(x:%):Complex(Float) ==
      retract map(#1::Complex(Float),
                  x pretend FE)$ExpressionFunctions2(Z, Complex Float)

)abbrev domain AN AlgebraicNumber
++ Algebraic closure of the rational numbers
++ Author: James Davenport
++ Date Created: 9 October 1995
++ Date Last Updated: 10 October 1995 (JHD)
++ Description: Algebraic closure of the rational numbers, with mathematical =
++ Keywords: algebraic, number.
AlgebraicNumber(): Exports == Implementation where
  Z   ==> Integer
  P   ==> SparseMultivariatePolynomial(Z, Kernel %)
  SUP ==>  SparseUnivariatePolynomial

  Exports ==> Join(ExpressionSpace, AlgebraicallyClosedField,
                   RetractableTo Z, RetractableTo Fraction Z,
                   LinearlyExplicitRingOver Z, RealConstant,
                   LinearlyExplicitRingOver Fraction Z,
                   CharacteristicZero,
                   ConvertibleTo Complex Float, DifferentialRing,
                   CoercibleFrom P) with
    numer  : % -> P
      ++ numer(f) returns the numerator of f viewed as a
      ++ polynomial in the kernels over Z.
    denom  : % -> P
      ++ denom(f) returns the denominator of f viewed as a
      ++ polynomial in the kernels over Z.
    reduce : % -> %
      ++ reduce(f) simplifies all the unreduced algebraic numbers
      ++ present in f by applying their defining relations.
    norm : (SUP(%),Kernel %) -> SUP(%)
      ++ norm(p,k) computes the norm of the polynomial p
      ++ with respect to the extension generated by kernel k
    norm : (SUP(%),List Kernel %) -> SUP(%)
      ++ norm(p,l) computes the norm of the polynomial p
      ++ with respect to the extension generated by kernels l
    norm : (%,Kernel %) -> %
      ++ norm(f,k) computes the norm of the algebraic number f
      ++ with respect to the extension generated by kernel k
    norm : (%,List Kernel %) -> %
      ++ norm(f,l) computes the norm of the algebraic number f
      ++ with respect to the extension generated by kernels l
  Implementation == InnerAlgebraicNumber add
    zero? a == trueEqual(rep a, rep 0)
    one? a == trueEqual(rep a, rep 1)
    a=b == trueEqual(rep a - rep b,rep 0)

