---- dalaw.spad

-- Dynamic algebraic laws

-- ATENCION TERESA : flag y context en dynev.spad

)abbrev category DALAWC    DynamicAlgebraicLawCategory
  ++ DynamicAlgebraicLawCategory(K:DFIELDC, Context:CTXTC)

)abbrev domain   DALAW     DynamicAlgebraicLaw
  ++ DynamicAlgebraicLaw(K:DFIELDC, Context:CTXTC)

--% DynamicAlgebraicLawCategory
 
  ++ A dynamic algebraic law la gives a computation rule
  ++ for its kernel a of the form: p(a)=0, 
  ++ where p is a univariate polynomial over the field K
  ++ called the DEFINING POLYNOMIAL of a (or of la)
  ++ It also gives a symbol that will represent a in outputs
  ++ p is monic, reduced, with >0 degree
  ++ and all its roots have the same multiplicity
  ++ If K has characteristic 0 then p is squarefree
  ++ If K is finite then p is squarefree and its prime factors 
  ++ have constant degree

  ++ Dynamic algebraic laws may be of 2 different types:
  ++   * linear if p has degree 1
  ++   * nonlinear otherwise
  ++ For a nonlinear law the polynomial p may be known to be
  ++ prime resp squarefree, 
  ++ then it is associated with the flag "prime" resp "squarefree",
  ++ else it is associated with the flag "noFlag"
  ++ If K has characteristic 0 or is finite, 
  ++ then the flag cannot be "noFlag"
  ++ The flag of a linear law is "prime"

  ++ In addition la gives some information :
  ++   * the multiplicity of a as a root of p,
  ++   * the degree of a over the underlying field,
  ++     if it is known (otherwise it is 0)
  ++   * some context (but the domain of contexts
  ++     may be empty, see DACLOS)

  ++ Two dynamic algebraic laws are roughly equal if and only if either:
  ++   * both are linear with roughly equal constants,
  ++   * or both are nonlinear with roughly equal polynomials
 
DynamicAlgebraicLawCategory(K,Context): Category == Exports where
 
  K: DynamicFieldCategory 
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  OUTFORM ==> OutputForm
  Poly ==> DynamicUnivariatePolynomial(K)
  Err ==> Union(trueLaw:$,failed:String)
  Split ==> Record(nzlaw:Err,zlaw:Err)
 
  Exports == DynamicSetCategory with
 
    -- operations
 
    -- no reduction, no splitting
 
      linear?: $ -> Boolean
        ++ linear?(la) returns true if and only if la is linear
      polynomial: $ -> Poly
        ++ polynomial(la) returns the defining polynomial of la
      constantIfLinear: $ -> K
        ++ constant(la) returns r such that 
        ++ the defining polynomial of la is x-r
        ++ Note: la must be linear
      polynomialIfNonlinear: $ -> Poly
        ++ polynomial(la) returns the defining polynomial of la
        ++ Note: la must be nonlinear
      refinePolynomial: ($,Poly) -> $
        ++ refinePolynomial(la,p) returns a law similar to la 
        ++ except for the polynomial that is set to p
        ++ Note: p should be a factor of the polynomial of la,
        ++ p must be monic and reduced
        ++ error if p has degree 0
      flag: $ -> Flag
        ++ flag(la) returns the flag of la
      refineFlag: ($,Flag) -> $
        ++ refineFlag(la,fl) returns a law similar to la 
        ++ except for the flag that is set to fl if fl>flag(la),
        ++ otherwise it returns la
      symbol: $ -> Symbol
        ++ symbol(la) returns the symbol of la
      multiplicity: $ -> NNI
        ++ multiplicity(la) returns the multiplicity part of la
      setMultiplicity: ($,NNI) -> $
        ++ setMultiplicity(la,nn) returns la except for
        ++ the multiplicity part that is set to nn
      degree: $ -> NNI
        ++ degree(la) returns the degree part of la
      setDegree: ($,NNI) -> $
        ++ setDegree(la,nn) returns la except for
        ++ the degree part that is set to nn
      context: $ -> Context
        ++ context(la) returns the context part of la
      setContext: ($,Context) -> $
        ++ setContext(la,ctxt) returns la except for
        ++ the context part that is set to ctxt
      build: (Poly,Symbol,Flag,NNI,NNI,Context) -> Err
        ++ build(p,fl,sy,mu,de,ct) returns the law associated with 
        ++ the polynomial p, the symbol sy, and the context ct
        ++ Note: p must be monic and reduced,
        ++ fl must say whatever is known about p,
        ++ mu must be the common multiplicity of every root of p in p, or 0,
        ++ de must be the common degree of every root of p over 
        ++ the current case, or 0
        ++ Note: returns "failed" if p has degree 0

    -- splitting
 
      split: ($,Poly) -> Split
        ++ split(la,q) returns (nzl,zl) where
        ++ nzl and zl are dynamic laws,
        ++ but one of them may be "failed"
        ++ this result is called a SPLITTING of la wrt q
        ++ it is such that "q(a)^=0" in nzl and "q(a)=0" in zl 
        ++ (where a is the kernel of la)
  
--% DynamicAlgebraicLaw
 
  ++ For a description of dynamic algebraic laws,
  ++ see comments for DynamicAlgebraicLawCategory
 
DynamicAlgebraicLaw(K,Context): Exports == Implementation where
 
  K: DynamicFieldCategory 
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  OUTFORM ==> OutputForm
  Poly ==> DynamicUnivariatePolynomial(K)
  Err ==> Union(trueLaw:$,failed:String)
  Split ==> Record(nzlaw:Err,zlaw:Err)
  Rul ==> Union(linear:K,nonlinear:Poly)
 
  Exports == DynamicAlgebraicLawCategory(K,Context)
 
  Implementation == add
 
    -- representation

      Rep:= Record(rul:Rul,sym:Symbol,fla:Flag,mul:NNI,deg:NNI,ctxt:Context)
 
      -- Let la be a law with defining polynomial p:Poly
      -- If la is linear, say p=x-a, 
      --   then the rule for la is the coefficient a:K,
      --   else the rule for la is the polynomial p
 
    -- declarations
      la,la1,la2: $
      p,q: Poly
      sy: Symbol
      fl: Flag
      nn,mu,de: NNI
      ct: Context

    -- define
 
    -- no reduction, no splitting
 
      roughEqual?(la1,la2) ==
        ru1:Rul:= la1.rul
        ru2:Rul:= la2.rul
        ru1 case linear =>
          (ru2 case linear) and (roughEqual?(ru1.linear,ru2.linear)$K)
        ru2 case linear =>
          false
        (roughEqual?(ru1.nonlinear,ru2.nonlinear)$Poly)
 
      linear?(la) == 
        la.rul case linear
 
      polynomial(la) ==
        la.rul case nonlinear => la.rul.nonlinear
        monomial(1,1)$Poly - monomial(la.rul.linear,0)$Poly
 
      constantIfLinear(la) == 
        la.rul.linear
 
      polynomialIfNonlinear(la) == 
        la.rul.nonlinear
 
      refinePolynomial(la,p) ==
        ans:$:= copy la
        n0:NNI:= roughDegree(p)$Poly
        n0 = 0 => 
          error "POLYNOMIAL MUST BE NON-CONSTANT IN refinePolynomial$DALAW"
        if n0 = 1 then
          r0:K:= -$K roughLeadingCoefficient(roughReductum(p)$Poly)$Poly
          ru0:Rul:= [r0]::Rul
        else
          ru0:Rul:= [p]::Rul
        ans.rul:= ru0
        ans
 
      flag(la) ==
        la.fla

      refineFlag(la,fl) ==
        ans:$:= copy la
        if ans.fla < fl then
          ans.fla:= fl
        ans

      symbol(la) ==
        la.sym

      multiplicity(la) == 
        la.mul

      setMultiplicity(la,nn) ==
        ans:$:= copy la
        ans.mul:= nn
        ans

      degree(la) == 
        la.deg

      setDegree(la,nn) ==
        ans:$:= copy la
        ans.deg:= nn
        ans

      context(la) == 
        la.ctxt

      setContext(la,ct) ==
        ans:$:= copy la
        ans.ctxt:= ct
        ans

      build(p,sy,fl,mu,de,ct) ==
        n0:NNI:= roughDegree(p)$Poly
        n0 = 0 => ["failed"]::Err
        if n0 = 1 then
          r0:K:= -$K roughLeadingCoefficient(roughReductum(p)$Poly)$Poly
          ru0:Rul:= [r0]::Rul
        else
          ru0:Rul:= [p]::Rul
        [[ru0,sy,fl,mu,de,ct]::Rep]::Err
 
    -- reduction, no splitting
 
      reduce(la) ==
        ans:$:= copy la
        ru0:Rul:= ans.rul
        if ru0 case linear then 
          ru0:= [reduce(ru0.linear)$K]::Rul
        else 
          ru0:= [reduce(ru0.nonlinear)$Poly]::Rul
        ans.rul:= ru0
        ans
 
      coerce(la):OUTFORM ==
        sy0:Symbol:= la.sym
        le0:List OUTFORM:=
          if la.rul case linear then
            [sy0::OUTFORM,"="::OUTFORM,(la.rul.linear)::OUTFORM]
          else 
            [output(la.rul.nonlinear,sy0)$Poly,"= 0"::OUTFORM]
        blankSeparate le0
 
    -- splitting
 
      split(la,q) ==
        ru0:Rul:= la.rul
        ru0 case linear =>
          if eval(q,ru0.linear)$Poly =$K 0
            then [["failed"]::Err,[la]::Err]::Split
          else [[la]::Err,["failed"]::Err]::Split
        p0:Poly:= ru0.nonlinear
        fl0:Flag:= la.fla
        fl0 = prime()$Flag =>
          if monicRemainder(q,p0)$Poly = 0 then
            [["failed"]::Err,[la]::Err]
          else [[la]::Err,["failed"]::Err]
        dp0:= 
          if fl0 = squarefree()$Flag then splitSquarefree(p0,q)$Poly
          else split(p0,q)$Poly
        d0:Poly:= dp0.gcdiv
        p0:= dp0.cofact1
        [build(p0,la.sym,fl0,la.mul,la.deg,la.ctxt),
          build(d0,la.sym,fl0,la.mul,la.deg,la.ctxt)]::Split
 

---- daiclos.spad

-- Dynamic algebraic internal closure

)abbrev category DAICLFC  DynamicAlgebraicInternalClosedFieldCategory
  ++ DynamicAlgebraicInternalClosedFieldCategory(Context:CTXTC) 

)abbrev category DAICLC  DynamicAlgebraicInternalClosureCategory
  ++ DynamicAlgebraicInternalClosureCategory(K0:DFIELDC,Context:CTXTC)

)abbrev domain   DAICL1  DynamicAlgebraicInternalClosure1
)abbrev domain   DAICL2  DynamicAlgebraicInternalClosure2
)abbrev domain   DAICL3  DynamicAlgebraicInternalClosure3
  ++ DAICL1, DAICL2 and DAICL3 are just parts of DAICL
  ++ DAICL1 is MUTABLE
)abbrev domain   DAICL   DynamicAlgebraicInternalClosure
  ++ DynamicAlgebraicInternalClosure(K0:DFIELDC,Context:CTXTC)


--% DynamicAlgebraicInternalClosedFieldCategory
 
DynamicAlgebraicInternalClosedFieldCategory(Context): Category == Exports where

  Context: ContextCategory

  NNI ==> NonNegativeInteger
  SPoly ==> SparseUnivariatePolynomial($)

  Exports == DynamicFieldCategory with
 
    -- operations
 
    -- splitting

      rootOf: (SPoly, Symbol) -> $
        ++ rootOf(p,a) adds a new level with law "p(a)=0"
        ++ and returns the correponding algebraic kernel
        ++ p is called the DEFINING POLYNOMIAL of this kernel
        ++ Error if p is constant
      rootOf: (SPoly, Symbol, Flag, NNI, NNI, Context) -> $
        ++ rootOf(p,sym,fl,mu,de,ct) adds a new level with law "p(a)=0"
        ++ and returns the corresponding algebraic kernel 
        ++ Note: This is a technical version of rootOf, sometimes useful
        ++   * fl says whatever is known about p,
        ++   * mu is the multiplicity of a over the current extension
        ++     if it is known, otherwise it must be 0,
        ++   * de is the degree of a over the current extension
        ++     if it is known, otherwise it must be 0,
        ++   * ct is any element in Context
        ++ If some of these informations are not known, use default values:
        ++   * noFlag for flag
        ++   * 0 for multiplicity and degree
        ++   * empty for context

--% DynamicAlgebraicInternalClosureCategory

DynamicAlgebraicInternalClosureCategory(K0,Context): Category == Exports where

  K0: DynamicFieldCategory
  Context: ContextCategory

  PI ==> PositiveInteger
  NNI ==> NonNegativeInteger
  SPoly ==> SparseUnivariatePolynomial($)
  Law ==> DynamicAlgebraicLaw($,Context)
  Tower ==> DynamicTower(Law)

  Exports == 
   Join(DynamicAlgebraicInternalClosedFieldCategory(Context), Algebra(K0)) with
 
    -- operations
 
    -- no reduction, no splitting

      level: () -> NNI
        ++ level() returns the level of the current case
      root: PI -> $
        ++ root(i) returns the ALGEBRAIC KERNEL of level i in $
      level: $ -> NNI
        ++ level(x) returns the level of x
      law: PI -> Law
        ++ law(l) returns the law of level l in the current case
      degree: NNI -> NNI
        ++ degree(nn) returns the degree of the algebraic
        ++ kernel of level nn over the basic case
        ++ if this degree is known, otherwise it returns 0
      degree: () -> NNI
        ++ degree() returns the degree of the maximal extension
        ++ in the current case if this degree is known,
        ++ (the degree is always known if $ is finite)
        ++ if the degree is not known then the answer is 0
      if K0 has Finite then 
        size: NNI -> NNI
          ++ size(nn) returns the size of the extension of level nn
          ++ in the current case
        size: () -> NNI
          ++ size() returns the size of the maximal extension 
          ++ in the current case
      roughRetractable?: $ -> Boolean
      roughRetract: $ -> K0
 
    -- splitting

      if $ has Finite then
        degree: ($,NNI) -> NNI
          ++ degree(x,nn) returns the smallest d such that x**(nn**d)=x
          ++ where q = size0**nn
      definingPolynomial: PI -> SPoly
        ++ definingPolynomial(i) returns the DEFINING POLYNOMIAL 
        ++ for the i-th kernel
      representingPolynomial: $ -> SPoly
        ++ representingPolynomial(x) returns a REPRESENTING POLYNOMIAL of x
        ++ If x has level l>0 it is a polynomial q such that x=q(a) 
        ++ where a is the algebraic kernel of level l
        ++ and the coefficients of q are algebraic numbers of level < l
        ++ If p denotes the defining polynomial of level l and
        ++ if x is reduced then q is uniquely defined mod p
        ++ and deg(q) < deg(p)
        ++ If the level of x is 0 then q is the constant polynomial x

    -- assertion

      if K0 has Finite then Finite
        ++ Note: of course the "mathematical" algebraic closure of
        ++ a finite field is NOT finite,
        ++ but here the "current extension" is always finite!

    -- "internal" operations

    -- no reduction, no splitting

      current : () -> Tower
      next : () -> List Tower
      setCurrent : Tower -> Void
      setNext : List Tower -> Void
      refreshCurrent : () -> Void
      refreshNext : () -> Void
      refresh : () -> Void
      addNext: Tower -> Void

--% DynamicAlgebraicInternalClosure1
 
  ++ Part of DynamicAlgebraicInternalClosure that is directly related
  ++ the mutable variables

DynamicAlgebraicInternalClosure1(K0,Context): Exports == Implementation where
 
  K0: DynamicFieldCategory
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  I ==> Integer
  Poly ==> DynamicUnivariatePolynomial($)
  SPoly ==> SparseUnivariatePolynomial($)
  Law ==> DynamicAlgebraicLaw($,Context)
  Tower ==> DynamicTower(Law)
   
  Exports == DynamicAlgebraicInternalClosureCategory(K0,Context) 
 
  Implementation == add
 
    -- representation

      CK:= Record(lev:PI,poly:Poly)
      Rep:= Union(rat:K0,alg:CK)
 
    -- about reduction
      -- let x:$ and i:NNI:=level(x) (i=0 <=> x case rat)
      -- if i=0 then x must be reduced as an element of K
      -- if i>0 then its representing polynomial q must be reduced in Poly
      -- and the level of x must be as small as possible
      -- which means that roughDegree(q) > 0
 
    -- mutable

      mutcas: Tower := empty()$Tower
        -- the value of mutcas is the current case for $
        -- its laws are reduced up to level(mutcas)
      mutnext: List Tower := []
        -- the value of mutnext is the list of the next cases to 
        -- be considered by the control package
        -- the towers in mutnext have level 0, and their laws
        -- are not reduced 
 
    -- declarations

      tow: Tower
      ltow: List Tower
 
    -- define
 
    -- no reduction, no splitting
 
      current() == 
        mutcas
 
      next() == 
        mutnext
 
      setCurrent(tow) ==
        mutcas:= tow
        void
     
      setNext(ltow) ==
        mutnext:= ltow
        void
     
      refreshCurrent() ==
        mutcas:= empty()$Tower
        void

      refreshNext() ==
        mutnext:= empty()$(List(Tower))
        void

      refresh() ==
        mutcas:= empty()$Tower
        mutnext:= empty()$(List(Tower))
        void

      addNext(tow) ==
        mutnext:= [tow,:mutnext]
        void

--% DynamicAlgebraicInternalClosure2
 
  ++ Part of DynamicAlgebraicInternalClosure that does not
  ++ make use of the mutable variables

DynamicAlgebraicInternalClosure2(K0,Context): Exports == Implementation where
 
  K0: DynamicFieldCategory
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  I ==> Integer
  Poly ==> DynamicUnivariatePolynomial($)
  Law ==> DynamicAlgebraicLaw($,Context)
  Tower ==> DynamicTower(Law)
   
  Exports == DynamicAlgebraicInternalClosureCategory(K0,Context) 

  Implementation == DynamicAlgebraicInternalClosure1(K0,Context) add
 
    -- representation

      CK:= Record(lev:PI,poly:Poly)
      Rep:= Union(rat:K0,alg:CK)
 
    -- about reduction

      -- see DAICLOS1

    -- local operations
 
    -- no reduction, no splitting

      roughZero?: $ -> Boolean

    -- declarations

      x,x1,x2: $
      r: K0
      np: PI
      nn: NNI
      n: I
 
    -- define
 
    -- no reduction, no splitting
 
      roughCoerce(x) ==
        x case rat => roughCoerce(x.rat)$K0
        la0:Law:= law(current,x.alg.lev)$Tower
        sy0:Symbol:= symbol(la0)$Law
        output(x.alg.poly,sy0)$Poly

      coerce(r):$ == [r]::Rep
 
      coerce(n):$ == [n::K0]::Rep
 
      roughRetractable?(x) == x case rat
 
      roughRetract(x) ==
        x case alg =>
          error "ARGUMENT OF roughRetract IS NOT roughRetractable"
        x.rat
 
      characteristic() == characteristic()$K0
 
      0 == [0$K0]::Rep
 
      1 == [1$K0]::Rep
 
      x1 + x2 ==
        x1 case rat =>
          if x2 case rat then [x1.rat +$K0 x2.rat]::Rep
          else [[x2.alg.lev,
            makePolynomial(x1)$Poly +$Poly x2.alg.poly]::CK]::Rep
        x2 case rat => [[x1.alg.lev,
          makePolynomial(x2)$Poly +$Poly x1.alg.poly]::CK]::Rep
        x1.alg.lev < x2.alg.lev =>
          [[x2.alg.lev, makePolynomial(x1)$Poly +$Poly x2.alg.poly]::CK]::Rep
        x2.alg.lev < x1.alg.lev =>
          [[x1.alg.lev, makePolynomial(x2)$Poly +$Poly x1.alg.poly]::CK]::Rep
        [[x1.alg.lev, x1.alg.poly +$Poly x2.alg.poly]::CK]::Rep
 
      x1 - x2 == x1 + (- x2)
 
      - x ==
        x case rat => [-$K0 (x.rat)]::Rep
        [[x.alg.lev, -$Poly x.alg.poly]::CK]::Rep
 
      n * x ==
        x case rat => [n *$K0 (x.rat)]::Rep
        [[x.alg.lev, n *$Poly x.alg.poly]::CK]::Rep
 
      r * x ==
        x case rat => [r *$K0 (x.rat)]::Rep
        [[x.alg.lev, (([r]::Rep)::$) *$Poly x.alg.poly]::CK]::Rep
 
      x1 * x2 ==
        x1 case rat =>
          if x2 case rat then [(x1.rat) *$K0 (x2.rat)]::Rep
          else (x1.rat) * x2
        x2 case rat => (x2.rat) * x1
        (x1.alg.lev < x2.alg.lev) =>
          [[x2.alg.lev, x1 *$Poly x2.alg.poly]::CK]::Rep
        (x2.alg.lev < x1.alg.lev) =>
          [[x1.alg.lev, x2 *$Poly x1.alg.poly]::CK]::Rep
        [[x1.alg.lev, x1.alg.poly *$Poly x2.alg.poly]::CK]::Rep
 
      power(x,nn) == 
      -- to improve when K0 has Finite
        nn=0 => 1
        nn=1 => x
        x case rat => [(x.rat) **$K0 nn]::Rep
        [[x.alg.lev,x.alg.poly **$Poly nn]::CK]::Rep
 
      roughZero?(x) ==
        x case rat => roughEqual?(x.rat,0)$K0
        roughEqual?(x.alg.poly,0)$Poly
 
      roughEqual?(x1,x2) == roughZero?(x1-x2)
 
      level(x) ==
        x case rat => 0
        (x.alg.lev)::NNI
 
      law(i) == law(current,i)$Tower

      root(np) == 
        [[np,monomial(1$$,1$NNI)$Poly]::CK]::Rep

      roughExquo(x1,x2) ==
        x1 case rat =>
          if x2 case rat then [roughExquo(x1.rat,x2.rat)$K0]::Rep
          else error "ERROR IN roughExquo$DACL (x1 case rat, x2 case alg)"
        (x2 case rat) or (x1.alg.lev > x2.alg.lev) =>
          p0:Poly:= roughExquo(x1.alg.poly,x2)$Poly
          [[x1.alg.lev, p0]::CK]::Rep
        x1.alg.lev < x2.alg.lev =>
          error "ERROR IN roughExquo$DACL (level(x1)<level(x2))"
        p0:Poly:= roughExquo(x1.alg.poly,x2.alg.poly)$Poly
        [[x1.alg.lev,p0]::CK]::Rep
 

--% DynamicAlgebraicInternalClosure3
 
  ++ Part of DynamicAlgebraicInternalClosure that does not fit 
  ++ into either one of DAICLOS1, DAICLOS2 or DAICLOS

DynamicAlgebraicInternalClosure3(K0,Context): Exports == Implementation where
 
  K0: DynamicFieldCategory
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  I ==> Integer
  Poly ==> DynamicUnivariatePolynomial($)
  SPoly ==> SparseUnivariatePolynomial($)
  Law ==> DynamicAlgebraicLaw($,Context)
  Tower ==> DynamicTower(Law)

  Exports == DynamicAlgebraicInternalClosureCategory(K0,Context) 

  Implementation == DynamicAlgebraicInternalClosure2(K0,Context) add
 
    -- representation

      CK:= Record(lev:PI,poly:Poly)
      Rep:= Union(rat:K0,alg:CK)
 
    -- about reduction

      -- see DAICLOS1 

    -- local operations
 
    -- reduction, no splitting

      refineAndReduceLaws: (PI,Law) -> Void
        ++ refineAndReduceLaws(np,la) refines law of level np to la
        ++ and reduces every law from level np+1 to level(current)

    -- declarations

      x,x1,x2: $
      np: PI
      nn: NNI
      n: I
 
    -- define
 
    -- no reduction, no splitting

      level() ==
        level(current)$Tower

      degree(nn) == 
        tow0:Tower:= current
        lev0:NNI:= level(tow0)$Tower
        nn>lev0 => error "LEVEL IS TOO LARGE IN degree$DAICLOS3"
        deg0:NNI:= 1$NNI
        for i in 1..nn repeat 
          la0:Law:= law(tow0,i pretend PI)
          d0:NNI:= degree(la0)$Law
          deg0:= deg0 ** d0
        deg0

      degree() ==
        tow0:Tower:= current
        lev0:NNI:= level(tow0)$Tower
        degree(lev0)
      
      if K0 has Finite then 

        size0() == size()$K0

        size(nn) == 
          tow0:Tower:= current
          lev0:NNI:= level(tow0)$Tower
          nn>lev0 => error "LEVEL IS TOO LARGE IN size$DAICLOS3"
          siz0:NNI:= size()$K0
          for i in 1..nn repeat 
            la0:Law:= law(tow0,i pretend PI)
            d0:NNI:= degree(la0)$Law
            siz0:= siz0 ** d0
          siz0

        size() ==
          tow0:Tower:= current
          lev0:NNI:= level(tow0)$Tower
          size(lev0)
      
    -- reduction, no splitting
 
      coerce(x):OutputForm ==
        roughCoerce(reduce(x))

      reduce(x) ==
        x case rat => (reduce(x.rat)$K0)::$
        p0:Poly:= x.alg.poly
        np0:PI:= x.alg.lev
        la0:Law:= law(current,np0)$Tower
        linear?(la0)$Law =>
          reduce(eval(p0,constantIfLinear(la0)$Law)$Poly)
        p0:= monicRemainder(p0,polynomialIfNonlinear(la0)$Law)$Poly
        p0:= reduce(p0)$Poly
        roughDegree(p0)$Poly=0 => roughLeadingCoefficient(p0)$Poly
        [[np0,p0]::CK]::$
 
      refineAndReduceLaws(np,la) ==
        tow0:Tower:= refineLaw(current,np,la)$Tower
        setCurrent tow0
        for i in (np+1)..level() repeat
          tow0:= reduceLaw(current,i pretend PI)$Tower
          setCurrent tow0
        if K0 has Finite then 
          la0:Law:= law(current,np)$Tower
          deg1:NNI:= degree(polynomial(la0)$Law)$Poly
          deg2:NNI:= degree(la0)$Law
          if deg1 = deg2 then
            la0:= refineFlag(la0,prime$Flag)$Law
            tow0:= refineLaw(current,np,la0)$Tower
            setCurrent tow0
        void

    -- splitting
 
      representingPolynomial(x) ==
        p0:Poly:=
          if x case rat then makePolynomial(x)$Poly
          else x.alg.poly
        p0::SPoly
 
      definingPolynomial(np) == 
        la0:Law:= law(current,np)$Tower
        (polynomial la0)::SPoly 

      x ** n ==
        n >= 0 => reduce(power(x,n pretend NNI))
        reduce(power(inv(x),(-n) pretend NNI))
 
      recip(x) ==
        x case rat =>
          if (r0:= recip(x.rat)$K0) case "failed" then "failed"
          else [r0::K0]::Rep::$
        x =$$ 0 => "failed"
        x0:$:= reduce(x)
        x0 case rat =>
          if (r0:= recip(x0.rat)$K0) case "failed" then "failed"
          else [r0::K0]::Rep::$
        roughDegree(x0.alg.poly)$Poly = 0 =>
          recip(roughLeadingCoefficient(x0.alg.poly)$Poly)
        -- now law(current,x0.alg.lev) is nonlinear.
        la0:Law:= law(current,x0.alg.lev)$Tower
        p0:Poly:= polynomialIfNonlinear(la0)$Law
        [[x0.alg.lev,inverseModulo(p0,x0.alg.poly)$Poly]::CK]::Rep::$
 
      inv(x) ==
        (ix:=recip(x)) case "failed" => error "DIVISION BY 0 IN inv$DACL"
        ix::$
 
      x1 / x2 ==
        reduce (x1 *$$ inv(x2))
 
      if K0 has Finite then
        degree(x,nn) ==
          ans:NNI:= 1
          x0:$:= power(x,nn)
          while x0 ^= x repeat
            ans:= ans + 1
            x0:= power(x0,nn)
          ans

      x1 = x2 ==
        x0:$:= reduce(x1-x2)
        x0 case rat => ((x0.rat) =$K0 0)
        np0:PI:= x0.alg.lev
        q0:Poly:= x0.alg.poly
        -- q0 is reduced with >0 degree since x0 is reduced
        la0:Law:= law(current,np0)$Tower
        sp0:= split(la0,q0)$Law
        sp0.nzlaw case failed =>
          if not(roughEqual?(la0,sp0.zlaw.trueLaw)$Law) then
            refineAndReduceLaws(np0,sp0.zlaw.trueLaw)
          true
        sp0.zlaw case failed =>
          if not(roughEqual?(la0,sp0.nzlaw.trueLaw)$Law) then
            refineAndReduceLaws(np0,sp0.nzlaw.trueLaw)
          false
        tow0:Tower:= refineLaw(current,np0,sp0.zlaw.trueLaw)$Tower
        addNext(tow0)
        refineAndReduceLaws(np0,sp0.nzlaw.trueLaw)
        false
 
--% DynamicAlgebraicInternalClosure
 
  ++ Part of DynamicAlgebraicInternalClosure that is directly 
  ++ related to rootOf

DynamicAlgebraicInternalClosure(K0,Context): Exports == Implementation where
 
  K0: DynamicFieldCategory
  Context: ContextCategory

  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  I ==> Integer
  Poly ==> DynamicUnivariatePolynomial($)
  SPoly ==> SparseUnivariatePolynomial($)
  Law ==> DynamicAlgebraicLaw($,Context)
  Tower ==> DynamicTower(Law)
   
  Exports == DynamicAlgebraicInternalClosureCategory(K0,Context) 

  Implementation == DynamicAlgebraicInternalClosure3(K0,Context) add
 
    -- representation

      CK:= Record(lev:PI,poly:Poly)
      Rep:= Union(rat:K0,alg:CK)
 
    -- about reduction

      -- see DAICLOS1

    -- local operations
 
    -- splitting

      easyFactorize: ($,Flag,NNI,NNI) -> $
        ++ easyFactorize(x,fl,mu,de) is called inside rootOf,
        ++ it returns the corresponding algebraic kernel
        ++ It makes some equality tests in order to get 
        ++ easy factorization of the corresponding defining polynomial 
        ++ Note: fl, mu and de are those in rootOf
        ++ and x must be the corresponding algebraic kernel
      easySquarefree: (Law,PI,NNI) -> Void
        ++ easySquarefree(la,l,mu) is used in easyFactorize, 
        ++ it makes the defining polynomial squarefree
        ++ Note: only when K0 has characteristic 0 and la is nonlinear
        ++ l is the level of la and mu is the one in rootOf
      easyDistinctDegree: ($,PI) -> Void
        ++ easyDistinctDegree(x,l) is used in easyFactorize, 
        ++ it makes a distinct degree factorization of the 
        ++ defining polynomial
        ++ Note: only when K0 is finite 
        ++ x is the kernel and l is the level of x
      easyDistinctMultiplicity: ($,Law,PI) -> Void
        ++ easyDistinctMultiplicity(x,la,l) is used in easyFactorize, 
        ++ It makes a distinct multiplicity factorization of the 
        ++ defining polynomial
        ++ Note: x is the kernel of la and l is its level

    -- declarations

      x,x1,x2: $
      l: PI
      mu,de: NNI
      n: I
      sup:SPoly
      ct: Context
      sym: Symbol
      fl: Flag
 
    -- define
 
    -- splitting
 
      easyFactorize(x,fl,mu,de) ==
      -- to improve
        lev0:PI:= level(current)$Tower pretend PI
        la0:Law:= law(current,lev0)$Tower
        linear?(la0)$Law =>
          x
        p0:Poly:= polynomialIfNonlinear(la0)$Law
        fl = prime()$Flag => 
          if (de = 0) then 
            la0:= setDegree(la0,roughDegree(p0)$Poly)$Law
          if (mu = 0) and (characteristic()$K0 = 0) then 
            la0:= setMultiplicity(la0,1$NNI)$Law
            tow0:Tower:= refineLaw(current,lev0,la0)$Tower
            setCurrent tow0
          x
        if (fl = noFlag()$Flag) and (characteristic()$K0 = 0) then
          easySquarefree(la0,lev0,mu)
        if (de = 0) and (K0 has Finite) then 
          easyDistinctDegree(x,lev0)
        if (fl = noFlag()$Flag) then
          easyDistinctMultiplicity(x,la0,lev0)
        x

      easySquarefree(la,l,mu) ==
        p0:Poly:= polynomialIfNonlinear(la)$Law 
        p0:= makeSquarefree(p0)$Poly
        la0:Law:= refinePolynomial(la,p0)$Law
        la0:= refineFlag(la0,squarefree$Flag)$Law
        if (mu^=1) then
          la0:= setMultiplicity(la0,1$NNI)$Law
        tow0:Tower:= refineLaw(current,l,la0)$Tower
        setCurrent tow0
        void

      if K0 has Finite then 
        easyDistinctDegree(x,l) ==
          nn0:NNI:= 
            if l=1 then size()$K0
            else size((l - 1) pretend NNI)
          d0:NNI:= degree(x,nn0)
          la0:Law:= law(current,l)$Tower
          la0:= setDegree(la0,d0)$Law
          fl0:Flag:= 
            if d0 = degree(polynomial(la0)$Law)$Poly then 
              prime$Flag
            else
              squarefree$Flag
          la0:= refineFlag(la0,fl0)$Law
          tow0:Tower:= refineLaw(current,l,la0)$Tower
          setCurrent tow0
          void

      easyDistinctMultiplicity(x,la,l) ==
        m0:NNI:= multiplicity(polynomial(la)$Law,x)$Poly
        la0:Law:= law(current,l)$Tower
        la0:= setMultiplicity(la0,m0)$Law
        tow0:Tower:= refineLaw(current,l,la0)$Tower
        setCurrent tow0
        void

      rootOf(sup,sym,fl,mu,de,ct) ==
        if maxLevel?(current)$Tower then
          p0:Poly:= sup::Poly
          p0 =$Poly 0 =>
            error "POLYNOMIAL IN rootOf MUST BE NON-ZERO"
          p0:= makeMonic(p0)$Poly
          errla0:= build(p0,sym,fl,mu,de,ct)$Law
          errla0 case failed =>
            error "POLYNOMIAL IN rootOf MUST BE NON-CONSTANT"
          la0:Law:= errla0.trueLaw
          tow0:Tower:= addLaw(current,la0)$Tower
          lev0:PI:= (level tow0) pretend PI
          setCurrent tow0
        else
          tow0:Tower:= addLevel(current)$Tower
          lev0:PI:= (level tow0) pretend PI
          tow0:= reduceLaw(tow0,lev0)$Tower
          setCurrent tow0
        x0:$:= root(lev0)
        easyFactorize(x0,fl,mu,de)  
 
      rootOf(sup,sym) ==
        rootOf(sup,sym,noFlag()$Flag,0$NNI,0$NNI,empty()$Context)


---- daclos.spad

-- Dynamic algebraic closure

-- ATENCION TERESA : VoidContext en dynev.spad

)abbrev category DACLFC  DynamicAlgebraicClosedFieldCategory
  ++ DynamicAlgebraicClosedFieldCategory

)abbrev category DACLOSC   DynamicAlgebraicClosureCategory
  ++ DynamicAlgebraicClosureCategory(K0:DFIELDC)

)abbrev domain   DACLOS     DynamicAlgebraicClosure
  ++ DynamicAlgebraicClosure(K0:FIELD)

)abbrev domain   DACASE     DynamicAlgebraicCase
  ++ DynamicAlgebraicCase(K0:FIELD) 


--% DynamicAlgebraicClosedFieldCategory
 
DynamicAlgebraicClosedFieldCategory: Category == Exports where

  NNI ==> NonNegativeInteger
  SPoly ==> SparseUnivariatePolynomial($)

  Exports == DynamicAlgebraicInternalClosedFieldCategory(VoidContext) with
 
    -- operations
 
    -- splitting

      rootOf: (SPoly, Symbol, Flag, NNI, NNI) -> $
        ++ ***

--% DynamicAlgebraicClosureCategory

DynamicAlgebraicClosureCategory(K0): Category == Exports where

  K0: Field

  DK0 ==> DynamicBuildField(K0)

  Exports == 
   Join(DynamicAlgebraicClosedFieldCategory, Algebra(K0),
     DynamicAlgebraicInternalClosureCategory(DK0,VoidContext))

--% DynamicAlgebraicClosure
 
DynamicAlgebraicClosure(K0): Exports == Implementation where
 
  K0: Field

  NNI ==> NonNegativeInteger
  DK0 ==> DynamicBuildField(K0)
  SPoly ==> SparseUnivariatePolynomial($)

  Exports == DynamicAlgebraicClosureCategory(K0)

  Implementation == DynamicAlgebraicInternalClosure(DK0,VoidContext) add

    -- representation

      Rep:= DynamicAlgebraicInternalClosure(DK0,VoidContext)

    -- declarations

      x0: K0
      sup: SPoly
      sym: Symbol
      fl: Flag
      mu,de: NNI

    -- define

    -- no reduction, no splitting

      coerce(x0):$ ==
        x0::DK0::Rep

    -- splitting

      rootOf(sup,sym,fl,mu,de) ==
        rootOf(sup,sym,fl,mu,de,empty()$VoidContext)$Rep
 
--% DynamicAlgebraicCase

DynamicAlgebraicCase(K0): Exports == Implementation where

  K0: Field

  K ==> DynamicAlgebraicClosure(K0)
  Law ==> DynamicAlgebraicLaw(K,VoidContext)
  Tower ==> DynamicTower(Law)
  OUTFORM ==> OutputForm

  Exports == DynamicCaseCategory

  Implementation == Tower add

    -- representation

      Rep:= Tower

    -- declarations

      ca: $
      lca: List $

    -- operations

      current() == 
        (current()$K)::$
 
      next() == 
        ltow0:List Tower:= next()$K
        [tow0::$ for tow0 in ltow0]

      setCurrent(ca) == 
        setCurrent(ca::Rep)$K

      setNext(lca) == 
        setNext([ca0::Rep for ca0 in lca])$K

      refreshCurrent() == 
        refreshCurrent()$K

      refreshNext() == 
        refreshNext()$K

      refresh() == 
        refresh()$K

      coerce(ca):OUTFORM ==
        cas0:$:= current()
        setCurrent(ca)
        out0:OUTFORM:= roughCoerce(ca::Rep)
        setCurrent(cas0)
        out0

