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

)abbrev package SFRGCD SquareFreeRegularTriangularSetGcdPackage
++ Author: Marc Moreno Maza
++ Date Created: 09/23/1998
++ Date Last Updated: 10/01/1998
++ References :
++  [1] M. MORENO MAZA and R. RIOBOO "Computations of gcd over
++      algebraic towers of simple extensions" In proceedings of AAECC11
++      Paris, 1995.
++  [2] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.
++  [3] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Description: 
++ A internal package for computing gcds and resultants of univariate 
++ polynomials with coefficients in a tower of simple extensions of a field.
++ There is no need to use directly this package since its main operations are
++ available from \spad{TS}. 

SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS) : SIG == CODE where
  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)

  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  S ==> String
  LP ==> List P
  PtoP ==> P -> P
  PS ==> GeneralPolynomialSet(R,E,V,P)
  PWT ==> Record(val : P, tower : TS)
  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Branch ==> Record(eq: List P, tower: TS, ineq: List P)
  UBF ==> Union(Branch,"failed")
  Split ==> List TS
  KeyGcd ==> Record(arg1: P, arg2: P, arg3: TS, arg4: B)
  EntryGcd ==> List PWT
  HGcd ==> TabulatedComputationPackage(KeyGcd, EntryGcd)
  KeyInvSet ==> Record(arg1: P, arg3: TS)
  EntryInvSet ==> List TS
  HInvSet ==> TabulatedComputationPackage(KeyInvSet, EntryInvSet)
  iprintpack ==> InternalPrintPackage()
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  quasicomppack ==> SquareFreeQuasiComponentPackage(R,E,V,P,TS)

  SQUAREFREE ==> SquareFreeRegularTriangularSetCategory(R,E,V,P)

  SIG ==> with

     startTableGcd! : (S,S,S) -> Void

     stopTableGcd! : () -> Void

     startTableInvSet! : (S,S,S) -> Void

     stopTableInvSet! : () -> Void

     stosePrepareSubResAlgo : (P,P,TS) -> List LpWT 

     stoseInternalLastSubResultant : (P,P,TS,B,B) -> List PWT 

     stoseInternalLastSubResultant : (List LpWT,V,B) -> List PWT 

     stoseIntegralLastSubResultant : (P,P,TS) -> List PWT 

     stoseLastSubResultant : (P,P,TS) -> List PWT 

     stoseInvertible? : (P,TS) -> B

     stoseInvertible?_sqfreg : (P,TS) -> List BWT

     stoseInvertibleSet_sqfreg : (P,TS) -> Split

     stoseInvertible?_reg : (P,TS) -> List BWT

     stoseInvertibleSet_reg : (P,TS) -> Split

     stoseInvertible? : (P,TS) -> List BWT

     stoseInvertibleSet : (P,TS) -> Split

     stoseSquareFreePart : (P,TS) -> List PWT 

  CODE ==> add

     startTableGcd!(ok: S, ko: S, domainName: S): Void == 
       initTable!()$HGcd
       printInfo!(ok,ko)$HGcd
       startStats!(domainName)$HGcd
       void()

     stopTableGcd!(): Void ==   
       if makingStats?()$HGcd then printStats!()$HGcd
       clearTable!()$HGcd

     startTableInvSet!(ok: S, ko: S, domainName: S): Void == 
       initTable!()$HInvSet
       printInfo!(ok,ko)$HInvSet
       startStats!(domainName)$HInvSet
       void()

     stopTableInvSet!(): Void ==   
       if makingStats?()$HInvSet then printStats!()$HInvSet
       clearTable!()$HInvSet

     stoseInvertible?(p:P,ts:TS): Boolean == 
       q := primitivePart initiallyReduce(p,ts)
       zero? q => false
       normalized?(q,ts) => true
       v := mvar(q)
       not algebraic?(v,ts) => 
         toCheck: List BWT := stoseInvertible?(p,ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val = false => return false
         return true
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       lgwt := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,true)
       for gwt in lgwt repeat
         g := gwt.val; 
         (not ground? g) and (mvar(g) = v) => 
           return false
       true

     stosePrepareSubResAlgo(p1:P,p2:P,ts:TS): List LpWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- ASSUME init(p1) invertible modulo ts !!!
       toSee: List LpWT := [[[p1,p2],ts]$LpWT]
       toSave: List LpWT := []
       v := mvar(p1)
       while (not empty? toSee) repeat
         lpwt := first toSee; toSee := rest toSee
         p1 := lpwt.val.1; p2 := lpwt.val.2
         ts := lpwt.tower
         lbwt := stoseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT)
         for bwt in lbwt repeat
           (bwt.val = true) and (degree(p2,v) > 0) =>
             p3 := prem(p1, -p2)
             s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N
             toSave := cons([[p2,p3,s],bwt.tower]$LpWT,toSave)
           -- p2 := initiallyReduce(p2,bwt.tower)
           newp2 := primitivePart initiallyReduce(p2,bwt.tower)
           (bwt.val = true) =>
             -- toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave)
             toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave)
           -- zero? p2 => 
           zero? newp2 => 
             toSave := cons([[p1,0,1],bwt.tower]$LpWT,toSave)
           -- toSee := cons([[p1,p2],bwt.tower]$LpWT,toSee)
           toSee := cons([[p1,newp2],bwt.tower]$LpWT,toSee)
       toSave

     stoseIntegralLastSubResultant(p1:P,p2:P,ts:TS): List PWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- ASSUME p1 and p2 have no algebraic coefficients
       lsr := lastSubResultant(p1, p2)
       ground?(lsr) => [[lsr,ts]$PWT]
       mvar(lsr) < mvar(p1) => [[lsr,ts]$PWT]
       gi1i2 := gcd(init(p1),init(p2))
       ex: Union(P,"failed") := (gi1i2 * lsr) exquo$P init(lsr)
       ex case "failed" => [[lsr,ts]$PWT]
       [[ex::P,ts]$PWT]
            
     stoseInternalLastSubResultant(p1:P,p2:P,ts:TS,b1:B,b2:B): List PWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- if b1 ASSUME init(p2) invertible w.r.t. ts
       -- if b2 BREAK with the first non-trivial gcd 
       k: KeyGcd := [p1,p2,ts,b2]
       e := extractIfCan(k)$HGcd
       e case EntryGcd => e::EntryGcd
       toSave: List PWT 
       empty? ts => 
         toSave := stoseIntegralLastSubResultant(p1,p2,ts)
         insert!(k,toSave)$HGcd
         return toSave
       toSee: List LpWT 
       if b1
         then
           p3 := prem(p1, -p2)
           s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N
           toSee := [[[p2,p3,s],ts]$LpWT]
         else
           toSee := stosePrepareSubResAlgo(p1,p2,ts)
       toSave := stoseInternalLastSubResultant(toSee,mvar(p1),b2)
       insert!(k,toSave)$HGcd
       toSave

     stoseInternalLastSubResultant(llpwt: List LpWT,v:V,b2:B): List PWT ==
       toReturn: List PWT := []; toSee: List LpWT; 
       while (not empty? llpwt) repeat
         toSee := llpwt; llpwt := []
         -- CONSIDER FIRST the vanishing current last subresultant
         for lpwt in toSee repeat 
           p1 := lpwt.val.1; 
           p2 := lpwt.val.2; 
           s := lpwt.val.3; 
           ts := lpwt.tower
           lbwt := stoseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT)
           for bwt in lbwt repeat
             bwt.val = false => 
               toReturn := cons([p1,bwt.tower]$PWT, toReturn)
               b2 and positive?(degree(p1,v)) => return toReturn
             llpwt := cons([[p1,p2,s],bwt.tower]$LpWT, llpwt)
         empty? llpwt => "leave"
         -- CONSIDER NOW the branches where the computations continue
         toSee := llpwt; llpwt := []
         lpwt := first toSee; toSee := rest toSee
         p1 := lpwt.val.1; p2 := lpwt.val.2; s := lpwt.val.3
         delta: N := (mdeg(p1) - degree(p2,v))::N
         p3: P := LazardQuotient2(p2, leadingCoefficient(p2,v), s, delta)
         zero?(degree(p3,v)) =>
           toReturn := cons([p3,lpwt.tower]$PWT, toReturn)
           for lpwt in toSee repeat 
             toReturn := cons([p3,lpwt.tower]$PWT, toReturn)
         (p1, p2) := (p3, next_subResultant2(p1, p2, p3, s))
         s := leadingCoefficient(p1,v)
         llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt)
         for lpwt in toSee repeat 
           llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt)
       toReturn

     stoseLastSubResultant(p1:P,p2:P,ts:TS): List PWT ==
       ground? p1 => 
         error"in stoseLastSubResultantElseSplit$SFRGCD  : bad #1"
       ground? p2 => 
         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
       not (mvar(p2) = mvar(p1)) => 
         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
       algebraic?(mvar(p1),ts) =>
         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1"
       not initiallyReduced?(p1,ts) => 
         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1"
       not initiallyReduced?(p2,ts) => 
         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
       purelyTranscendental?(p1,ts) and purelyTranscendental?(p2,ts) =>
         stoseIntegralLastSubResultant(p1,p2,ts)
       if mdeg(p1) < mdeg(p2) then 
          (p1, p2) := (p2, p1)
          if odd?(mdeg(p1)) and odd?(mdeg(p2)) then p2 := - p2
       stoseInternalLastSubResultant(p1,p2,ts,false,false)

     stoseSquareFreePart_wip(p:P, ts: TS): List PWT ==
     -- ASSUME p is not constant and mvar(p) > mvar(ts)
     -- ASSUME init(p) is invertible w.r.t. ts
     -- ASSUME p is mainly primitive
       mdeg(p) = 1 => [[p,ts]$PWT]
       v := mvar(p)$P
       q: P := mainPrimitivePart D(p,v)
       lgwt: List PWT := stoseInternalLastSubResultant(p,q,ts,true,false)
       lpwt : List PWT := []
       sfp : P
       for gwt in lgwt repeat
         g := gwt.val; us := gwt.tower
         (ground? g) or (mvar(g) < v) =>
           lpwt := cons([p,us],lpwt)
         g := mainPrimitivePart g
         sfp := lazyPquo(p,g)
         sfp := mainPrimitivePart stronglyReduce(sfp,us)
         lpwt := cons([sfp,us],lpwt)
       lpwt

     stoseSquareFreePart_base(p:P, ts: TS): List PWT == [[p,ts]$PWT]

     stoseSquareFreePart(p:P, ts:TS): List PWT == stoseSquareFreePart_wip(p,ts)
       
     stoseInvertible?_sqfreg(p:P,ts:TS): List BWT ==
       --iprint("+")$iprintpack
       q := primitivePart initiallyReduce(p,ts)
       zero? q => [[false,ts]$BWT]
       normalized?(q,ts) => [[true,ts]$BWT]
       v := mvar(q)
       not algebraic?(v,ts) => 
         lbwt: List BWT := []
         toCheck: List BWT := stoseInvertible?_sqfreg(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => lbwt := cons(bwt,lbwt)
           newq := removeZero(q,bwt.tower)
           zero? newq => lbwt := cons(bwt,lbwt)
           lbwt := 
            concat(stoseInvertible?_sqfreg(newq,bwt.tower)@(List BWT), lbwt)
         return lbwt
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
       lbwt: List BWT := []
       lts, lts_g, lts_h: Split
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           lts := augment(ts_v,ts)$TS
           lts := augment(members(ts_v_+),lts)$TS
           for ts in lts repeat
             lbwt := cons([true, ts]$BWT,lbwt)
         g := mainPrimitivePart g
         lts_g := augment(g,ts)$TS
         lts_g := augment(members(ts_v_+),lts_g)$TS
         -- USE stoseInternalAugment with parameters ??
         for ts_g in lts_g repeat
           lbwt := cons([false, ts_g]$BWT,lbwt)
         h := lazyPquo(ts_v,g)
         (ground? h) or (mvar(h) < v) => "leave"
         h := mainPrimitivePart h
         lts_h := augment(h,ts)$TS
         lts_h := augment(members(ts_v_+),lts_h)$TS
         -- USE stoseInternalAugment with parameters ??
         for ts_h in lts_h repeat
           lbwt := cons([true, ts_h]$BWT,lbwt)
       sort((x,y) +-> x.val < y.val,lbwt)

     stoseInvertibleSet_sqfreg(p:P,ts:TS): Split ==
       --iprint("*")$iprintpack
       k: KeyInvSet := [p,ts]
       e := extractIfCan(k)$HInvSet
       e case EntryInvSet => e::EntryInvSet
       q := primitivePart initiallyReduce(p,ts)
       zero? q => []
       normalized?(q,ts) => [ts]
       v := mvar(q)
       toSave: Split := []
       not algebraic?(v,ts) => 
         toCheck: List BWT := stoseInvertible?_sqfreg(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => toSave := cons(bwt.tower,toSave)
           newq := removeZero(q,bwt.tower)
           zero? newq => "leave"
           toSave := concat(stoseInvertibleSet_sqfreg(newq,bwt.tower), toSave)
         toSave := removeDuplicates toSave
         return algebraicSort(toSave)$quasicomppack
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
       lts, lts_h: Split
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           lts := augment(ts_v,ts)$TS
           lts := augment(members(ts_v_+),lts)$TS
           toSave := concat(lts,toSave)
         g := mainPrimitivePart g
         h := lazyPquo(ts_v,g)
         h := mainPrimitivePart h
         (ground? h) or (mvar(h) < v) => "leave"
         lts_h := augment(h,ts)$TS
         lts_h := augment(members(ts_v_+),lts_h)$TS
         toSave := concat(lts_h,toSave)
       toSave := algebraicSort(toSave)$quasicomppack
       insert!(k,toSave)$HInvSet
       toSave
       
     stoseInvertible?_reg(p:P,ts:TS): List BWT ==
       --iprint("-")$iprintpack
       q := primitivePart initiallyReduce(p,ts)
       zero? q => [[false,ts]$BWT]
       normalized?(q,ts) => [[true,ts]$BWT]
       v := mvar(q)
       not algebraic?(v,ts) => 
         lbwt: List BWT := []
         toCheck: List BWT := stoseInvertible?_reg(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => lbwt := cons(bwt,lbwt)
           newq := removeZero(q,bwt.tower)
           zero? newq => lbwt := cons(bwt,lbwt)
           lbwt := 
            concat(stoseInvertible?_reg(newq,bwt.tower)@(List BWT), lbwt)
         return lbwt
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
       lbwt: List BWT := []
       lts, lts_g, lts_h: Split
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           lts := augment(ts_v,ts)$TS
           lts := augment(members(ts_v_+),lts)$TS
           for ts in lts repeat
             lbwt := cons([true, ts]$BWT,lbwt)
         g := mainPrimitivePart g
         lts_g := augment(g,ts)$TS
         lts_g := augment(members(ts_v_+),lts_g)$TS
         -- USE internalAugment with parameters ??
         for ts_g in lts_g repeat
           lbwt := cons([false, ts_g]$BWT,lbwt)
         h := lazyPquo(ts_v,g)
         (ground? h) or (mvar(h) < v) => "leave"
         h := mainPrimitivePart h
         lts_h := augment(h,ts)$TS
         lts_h := augment(members(ts_v_+),lts_h)$TS
         -- USE internalAugment with parameters ??
         for ts_h in lts_h repeat
           inv := stoseInvertible?_reg(q,ts_h)@(List BWT)
           lbwt := concat([bwt for bwt in inv | bwt.val],lbwt)
       sort((x,y) +-> x.val < y.val,lbwt)

     stoseInvertibleSet_reg(p:P,ts:TS): Split ==
       --iprint("/")$iprintpack
       k: KeyInvSet := [p,ts]
       e := extractIfCan(k)$HInvSet
       e case EntryInvSet => e::EntryInvSet
       q := primitivePart initiallyReduce(p,ts)
       zero? q => []
       normalized?(q,ts) => [ts]
       v := mvar(q)
       toSave: Split := []
       not algebraic?(v,ts) =>
         toCheck: List BWT := stoseInvertible?_reg(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => toSave := cons(bwt.tower,toSave)
           newq := removeZero(q,bwt.tower)
           zero? newq => "leave"
           toSave := concat(stoseInvertibleSet_reg(newq,bwt.tower), toSave)
         toSave := removeDuplicates toSave
         return algebraicSort(toSave)$quasicomppack
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
       lts, lts_h: Split
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           lts := augment(ts_v,ts)$TS
           lts := augment(members(ts_v_+),lts)$TS
           toSave := concat(lts,toSave)
         g := mainPrimitivePart g
         h := lazyPquo(ts_v,g)
         h := mainPrimitivePart h
         (ground? h) or (mvar(h) < v) => "leave"
         lts_h := augment(h,ts)$TS
         lts_h := augment(members(ts_v_+),lts_h)$TS
         for ts_h in lts_h repeat
           inv := stoseInvertibleSet_reg(q,ts_h)
           toSave := removeDuplicates concat(inv,toSave)
       toSave := algebraicSort(toSave)$quasicomppack
       insert!(k,toSave)$HInvSet
       toSave

     if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P)
     then

       stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_sqfreg(p,ts)

       stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_sqfreg(p,ts)

     else
       
       stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_reg(p,ts)
 
       stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_reg(p,ts)