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

)abbrev domain PACRAT PseudoAlgebraicClosureOfRationalNumber
++ Authors: Gaetan Hache
++ Date Created: feb 1997
++ Date Last Updated: May 2010 by Tim Daly
++ Description: 
++ This domain implements dynamic extension using the simple notion of 
++ tower extensions. A tower extension T of the ground field K is any 
++ sequence of field extension (T : K_0, K_1, ..., K_i...,K_n) where K_0 = K 
++ and for i =1,2,...,n, K_i is an extension of K_{i-1} of degree > 1 and 
++ defined by an irreducible polynomial p(Z) in K_{i-1}.
++ Two towers (T_1: K_01, K_11,...,K_i1,...,K_n1) and 
++ (T_2: K_02, K_12,...,K_i2,...,K_n2) are said to be related if T_1 <= T_2 
++ (or T_1 >= T_2), that is if K_i1 = K_i2 for i=1,2,...,n1 
++ (or i=1,2,...,n2). Any algebraic operations defined for several elements 
++ are only defined if all of the concerned elements are comming from 
++ a set of related tour extensions. 

PseudoAlgebraicClosureOfRationalNumber() : SIG == CODE where
  
  INT     ==> Integer
  K       ==> Fraction Integer
  NNI     ==> NonNegativeInteger
  SUP     ==> SparseUnivariatePolynomial
  BOOLEAN ==> Boolean
  PI      ==> PositiveInteger
  FACTRN  ==> FactorisationOverPseudoAlgebraicClosureOfRationalNumber

  recRep  ==> Record(recEl:SUP(%),_
                    recTower:SUP(%),_
                    recDeg:PI,_
                    recPrevTower:%,_
                    recName:Symbol) 
    
  SIG ==> PseudoAlgebraicClosureOfRationalNumberCategory with 
  
    fullOutput : % -> OutputForm

    newElement : (SUP(%), SUP(%), PI,  %, Symbol) -> %
  
  CODE ==> add

    Rep := Union(recRep,K)
    
    -- signature of local function  
    replaceRecEl: (%,SUP(%)) -> %
    down: % -> %
    
    down(a:%) == 
      a case K => a
      aa:=(a pretend recRep)
      elel := aa.recEl
      ^ground?(elel)$SUP(%) => a
      gel:%:=ground(elel)
      down(gel)
      
    coerce(a:Integer):%== (a :: K)

    n:INT * a:% ==
      one?(n) => a
      zero?(a) or zero?(n) => 0
      (n < 0) => - ((-n)*a)
      mm:PositiveInteger:=(n pretend PositiveInteger)
      double(mm,a)$RepeatedDoubling(%)

    replaceRecEl(a,el)==
      a case K => a
      aa:=copy a
      aa.recEl := el
      aa

    -- local variable    
    localTower :% := 1$K
        
    -- implemetation of exported function

    lift(a) == 
      a case K => monomial(a,0)
      (a pretend recRep).recEl

    lift(a,b)==
      extDegree a > extDegree b => _
          error "Cannot lift something at lower level !!!!!"
      extDegree a < extDegree b => monomial(a,0)$SUP(%)
      lift a

    reduce(a)==
      localTower case K => 
        coefficient(a,0)
      ar:= a rem (localTower pretend recRep).recTower
      replaceRecEl(localTower,ar)
       
    maxTower(la)==
        --return an element from the list la which is in the largest
        --extension of the ground field
        --PRECONDITION: all elements in same tower, else no meaning?
      m:="max"/[extDegree(a)$% for a in la]
      first [b for b in la | extDegree(b)=m]

    ground?(a)== a case K

    vectorise(a,lev)==
      da:=extDegree a
      dlev:=extDegree lev
      dlev < da => _
        error "Cannot vectorise at a lower level than the element to vectorise"
      lev case K => [a]
      pa:SUP(%)
      na:%
      ^(da = dlev) =>
        pa:=  monomial(a,0)$SUP(%)
        na:=  replaceRecEl(lev,pa)
        vectorise(na,lev)$%
      prevLev:=previousTower(lev)
      a case K => error "At this point a is not suppose to be in K"
      aEl:=(a pretend recRep).recEl
      daEl:=degree definingPolynomial(a)$%
      lv:=[vectorise(c,prevLev)$% for c in entries(vectorise(aEl,daEl)$SUP(%))]
      concat lv        

    setTower!(a) ==
      localTower:=a
      void()
      
    definingPolynomial == definingPolynomial(localTower)
         
    a:% + b:% ==
      (a case K) and (b case K) => a +$K b
      extDegree(a) > extDegree(b) => b + a
      res1:SUP(%)
      res2:%
      if extDegree(a) = extDegree(b) then
        res1:=   b.recEl +$SUP(%) a.recEl
        res2:=   replaceRecEl(b,res1)
      else
        res1:=   b.recEl +$SUP(%) monomial(a,0)$SUP(%)
        res2:= replaceRecEl(b,res1)
      down(res2)
         
    a:% * b:% ==
      (a case K) and (b case K) => a *$K b
      extDegree(a) > extDegree(b) => b * a
      res1:SUP(%)
      res2:%
      if extDegree(a) = extDegree(b) then
        res1:=   b.recEl *$SUP(%) a.recEl rem b.recTower
        res2:=   replaceRecEl(b,res1)
      else
        res1:=   b.recEl *$SUP(%) monomial(a,0)$SUP(%)
        res2:=  replaceRecEl(b,res1)
      down(res2)
      
    distinguishedRootsOf(polyZero,ee) ==
        setTower!(ee)
        zero?(polyZero) => error "to lazy to give you all the roots of 0 !!!"
        factorf: Factored SUP % :=  factor(polyZero,ee)$FACTRN(%)
        listFact:List SUP %  := [pol.fctr for pol in factorList(factorf)]
        listOfZeros:List(%):=empty()
        for p in listFact repeat
          root:=newElement(p, new(D::Symbol)$Symbol)
          listOfZeros:List(%):=concat([ root ], listOfZeros)
        listOfZeros

    1 == 1$K

    0 == 0$K

    newElement(pol:SUP(%),subF:%,inName:Symbol): % ==
        -- pol is an irreducible polynomial over the field extension
        -- given by subF. 
        -- The output of this function is a root of pol.
      dp:=degree pol
      one?(dp) =>
        listCoef:=coefficients(pol)
        one?(#listCoef) => 0
        - last(listCoef) / first(listCoef)
      ground?(pol) => error "Cannot create a new element with a constant"
      d:PI := (dp pretend PI) * extDegree(subF)
      [monomial(1$%,1),pol,d,subF,inName] :: Rep
      
    newElement(poll:SUP(%),inName:Symbol)==
      newElement(poll,localTower,inName)

    newElement(elPol:SUP(%),pol:SUP(%),d:PI,subF:%,inName:Symbol): % ==
      [elPol, pol,d,subF,inName] :: Rep

    --Field operations 
    inv(a)==
      a case K => inv(a)$K
      aRecEl:= (a pretend recRep).recEl
      aDefPoly:= (a pretend recRep).recTower
      aInv := extendedEuclidean( aRecEl , aDefPoly, 1 )
      aInv  case "failed" => error "PACOFF : division by zero"
      -- On doit retourner un Record représentant l'inverse de a.
      -- Ce Record est exactement le même que celui de a sauf
      -- qu'il faut remplacer le polynôme du selecteur recEl
      -- par le polynôme représentant l'inverse de a :
      -- C'est ce que fait la fonction replaceRecEl.
      replaceRecEl( a , aInv.coef1 )     

    a:% / b:% == a * inv(b)
    
    a:K * b:%==
      (a :: %) * b
      
    b:% * a:K == a*b

    a:% - b:% ==
      a + (-b)
    
    a:% * b:Fraction(Integer) ==
      bn:=numer b
      bd:=denom b
      ebn:%:= bn * 1$%
      ebd:%:= bd * 1$%
      a * ebn * inv(ebd)

    -a:% ==
       a case K => -$K a 
       [-$SUP(%) (a pretend recRep).recEl,_
        (a pretend recRep).recTower,_
        (a pretend recRep).recDeg,_
        (a pretend recRep).recPrevTower,_
        (a pretend recRep).recName ]
        
    bb:% = aa:% ==
      b:=down bb
      a:=down aa
      ^( extDegree(b) =$NNI extDegree(a) ) => false
      (b case K)  =>  ( (retract a)@K  =$K (retract b)@K )
      rda := a :: recRep
      rdb := b :: recRep
      not (rda.recTower =$SUP(%) rdb.recTower) => false
      rdb.recEl =$SUP(%) rda.recEl
        
    zero?(a:%) == 
      da:=down a  -- just to be sure !!!
      ^(da case K) => false
      zero?(da)$K
    
    one?(a:%) ==
      da:= down a  -- just to be sure !!!
      ^(da case K) => false
      one?(da)$K
    
    --Coerce Functions
      
    coerce(a:K):% == a       

    retractIfCan(a:%):Union(Integer,"failed")==
      a case K => retractIfCan(a)$K
      "failed"
         
    retractIfCan(a:%):Union(K,"failed")==
      a case K => a
      "failed"
               
    coerce(a:%):OutputForm ==
      a case K => ((retract a)@K) ::OutputForm
      outputForm((a pretend recRep).recEl,_
                ((a pretend recRep).recName)::OutputForm) $SUP(%)

    fullOutput(a:%):OutputForm==
      a case K => ((retract a)@K) ::OutputForm
      (a pretend recRep)::OutputForm
      
    definingPolynomial(a:%): SUP % ==
      a case K => monomial(1,1)$SUP(%)
      (a pretend recRep).recTower

    extDegree(a:%): PI ==
      a case K => 1
      (a pretend recRep).recDeg
      
    previousTower(a:%):% ==
      a case K => error "No previous extension for ground field element"
      (a pretend recRep).recPrevTower
 
    name(a:%):Symbol ==
      a case K => error "No name for ground field element"
      (a pretend recRep).recName

    -- function related to the ground field 
    
    characteristic == characteristic()$K