% $OpenXM: OpenXM/src/kan96xx/Doc/complex.sm1,v 1.6 2000/12/29 07:19:37 takayama Exp $
%% lib/complex.sm1  [ functions for complex ], 1999, 9/9
%% cf.  yama:1999/Int/uli.sm1
%%%%%%%%%%%%%%%%%%%   commands %%%%%%%%%%%%%%%%%%%%%%%%%
%%%  res-div, res-solv, res-kernel-image, res-dual
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[(complex.sm1 : 1999, 9/28, res-div, res-solv, res-kernel-image, res-dual )
 (              2000, 6/8,  isExact_h, isExact )
 (In this package, complex is expressed in terms of matrices.) 
] {message-quiet} map
/uli.verbose 0 def
/uli.weight [(x) -1 (y) -1 (Dx) 1 (Dy) 1] def

%%% M = [M_1, ..., M_p],  M_i has the length q
%%%    D^p (row vector) --- M ---> D^q (row vector),   v --> v M
%%% In this package (res-***), all data are expressed by matrices.
/res-nextShift {
  /arg1 set
  [/in-nextShift /f /mm /m /p /q /i /fi] pushVariables
  [
      /f arg1 0  get def
      /mm arg1 1 get def
      %%  D^p[m] ---f---> D^q[mm]   [f mm] nextShift m
      /p f length def
      [1 1 p { pop 0 } for] /m set
      0 1 p 1 sub {
        /i set
        /fi f i get def
        m i  <<  mm  fi { uli.weight ord_w } map add maxInArray >>  put
      } for
      /arg1 m def
  ] pop
  popVariables
  arg1
} def

[(res-nextShift)
[([f mm] nextShift m)
 $Example: [(x,y) ring_of_differential_operators 0] define_ring$
 $ [ [ [ (x). (x^2). (x^3). ] $
 $     [ (Dx). (Dx^2). (Dx^3).]] [5 6 7]] res-nextShift :: $
]] putUsages


%% Input must be a matrix.
/res-init {
   /arg1 set
   [/in-initv /v /n] pushVariables
   [
     /v arg1 def
     /n v length def
     [n [v] fromVectors {init} map] toVectors2
     /arg1 set
   ] pop
   popVariables
   arg1
} def


/res-isVadapted {
  /arg1 set
  [/in-res-isVstrict /f /m /mm /ans] pushVariables
  [
    /f arg1 0 get def
    /m arg1 1 get def
    /mm arg1 2 get def
    %%  D^p[m] ---f---> D^q[mm]   [f m mm] res-isVadapted
    f [ [ ] ] eq { 
      /ans 1 def
    } {
      [f mm] res-nextShift m eq {/ans 1 def} { /ans 0 def} ifelse
    } ifelse
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

/res-gb {
  /arg1 set
  [/in-res-gb /aa /gg /qq /ans] pushVariables
  [(KanGBmessage)] pushEnv
  [
    /aa arg1 def  %% Input is a matrix.
    aa [ ] eq { /arg1 [ ] def /res-gb.LLL goto } {  } ifelse
    aa 0 get isArray {
    }{ aa { [ 2 1 roll ] } map /aa} ifelse
    /qq aa 0 get length def
    aa { dehomogenize homogenize } map /aa set
    uli.verbose { } { [(KanGBmessage) 0] system_variable} ifelse
    [aa] groebner 0 get /ans set
    ans 0 get isArray { }
    { [qq ans] toVectors2 /ans set } ifelse
    /arg1 ans def
    /res-gb.LLL
  ] pop
  popEnv
  popVariables
  arg1
} def

%% Utility functions res-setRing and res-toString
/res-toString {
  /arg1 set
  [/in-res-toString /s /ans] pushVariables
  [
    /s arg1 def
    s isArray {
      s {res-toString} map /ans set
    }{
      s isPolynomial {
        /ans s toString def
      } { 
        /ans s def 
      } ifelse
    } ifelse
    ans /arg1 set
  ] pop
  popVariables
  arg1
} def

%% res-setRing.v  res-setRing.vlist are global variables that contain,
%% for example, (x,y) and [(x) (y)].
/res-setRing {
  /arg1 set
  [/in-res-setRing /R /v] pushVariables
  [
     /v arg1 def
     v isArray {
       /v v res-toString from_records def
     }{
       v isString {
       }{
         [(res-setRing: ) v toString 
          ( is not a set of variables to define a ring.)] cat
         error
       }ifelse
     }ifelse
     /res-setRing.v v def
     /res-setRing.vlist [v to_records pop] def
     [v ring_of_differential_operators 0] define_ring /R set
     /arg1 R def
   ] pop
   popVariables
   arg1
} def
     

%% [M N] res-div  It returns ker(M/N)  i.e. D^*/ [M N] res-div = M/N
%% First size(M) part of the syzygy of M and N.
/res-div {
  /arg1 set
  [/in-res-div /M /N /ss /m /n /ss2 /ans] pushVariables
  [(KanGBmessage)] pushEnv
  [
    /M arg1 0 get def
    /N arg1 1 get def
    /m M length def
    /n N length def
    M 0 get isArray {
    }{ M { [ 2 1 roll ] } map /M } ifelse
    M { dehomogenize homogenize } map /M set

    n 0 eq not {
      N 0 get isArray {
      }{ N { [ 2 1 roll ] } map /N } ifelse
      N { dehomogenize homogenize } map /N set
    } { } ifelse

    uli.verbose { } { [(KanGBmessage) 0] system_variable} ifelse
    [M N join [(needSyz)]] groebner 2 get /ss set
    ss dehomogenize /ss set
    ss { [ 2 1 roll  aload pop 1 1 n { pop pop } for ] } map
    /ss2 set
    ss2 {homogenize} map /ss2 set
    ss2 [ ] eq {
      [ m res-newpvec ] /ans set
    }{
      [ss2 0 get length [ss2] groebner 0 get dehomogenize ] toVectors2
      /ans set
    } ifelse

    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
[(res-div)
[( [M N] res-div K )
 ( matrix M, N, K ; Each element of M and N must be an element of a ring.)
 ( coker(K) is isomorphic to M/N. )
 (Example: [(x,y) ring_of_differential_operators 0] define_ring )
 (   [[[(x+x^2+y^2).] [(x y).]] [[(x+x^2+y^2).] [(x y).]]] res-div)
 (  )
 $res*div accepts string inputs, too. For example,$
 $ [[[[(x+x^2+y^2)] [(x y)]] [[(x+x^2+y^2)] [(x y)]]]$
 $   [(x) (y)]]  res*div ::$
 (See also res-toString, res-setRing.)
]] putUsages

/res*div {
  /arg1 set
  [/in-res*div /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get res-toString expand res-div /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

/res-syz {
  /arg1 set
  [/in-res-syz /M /m] pushVariables
  [
    /M arg1 def

    M 0 get isArray {
    }{ M { [ 2 1 roll ] } map /M } ifelse

    M { dehomogenize homogenize } map /M set
    [M [(needSyz)]] groebner 2 get dehomogenize /arg1 set
  ] pop
  popVariables
  arg1
} def
[(res-syz) 
[( M res-syz N)
 ( matrix M, N ; each element of M and N must be an element of a ring.)
 ( N is a set of generators of the syzygy module of M.)
 (res*syz is also provided. It accepts string inputs.)
]] putUsages
/res*syz {
  /arg1 set
  [/in-res*syz /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get res-toString expand res-syz /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

/res-getx {
  /arg1 set
  [/in-res-getx /xx /nn /ff] pushVariables
  [
    /ff arg1 def
    /xx ff getvNamesCR def  
    [(N)] system_variable /nn set
    [ xx aload pop 1 1 nn { pop pop } for pop ] rest
    /arg1 set
  ] pop
  popVariables
  arg1
} def

%% Solving \sum c_i M_i = d
%% [M d] res-solv c'/r  ;   M : matrix,  d, c' : vectors, r : scalar, c'/r =c
/res-solv {
  /arg1 set
  [/in-res-solv /M /d /ans /B /vv /G /rr /rng /nn] pushVariables
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [
     /nn arg1 length def
     /M arg1  0 get def
     /d arg1  1 get def
     nn 3 eq {
       /rng arg1 2 get def
     }{
       M getRing /rng set
       rng tag RingP eq {  }
       { d getRing /rng set } ifelse
     }ifelse
     rng res-getx /vv set
     uli.verbose { (res-solv : vv = ) messagen vv message } { } ifelse
     uli.verbose { } { [(KanGBmessage) 0] system_variable } ifelse
     M dehomogenize /M set
     [vv from_records ring_of_differential_operators 0] define_ring
     M 0 get isArray {
       M { { toString . } map } map /M set
     } {
       M { toString . } map /M set
     } ifelse
     [M [(needBack)]] groebner_sugar /G set
     G 1 get /B set

     d isArray {
       d 0 get isArray { [d] fromVectors 0 get /d set } { } ifelse
       [d] fromVectors 0 get /d set
     } {  } ifelse
     d toString . dehomogenize /d set

     /res-solv.d d def
     /res-solv.G G def

     d G 0 get reduction-noH  /rr set
     rr 0 get (0). eq {
       [rr 2 get] B mul 0 get /ans set
       /ans [ ans { toString rng ,, (-1) rng ,, mul} map  
              rr 1 get toString .. ] def
     } {
       /ans null def
     } ifelse
     /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
[(res-solv)
[$[M d] res-solv [c' r] $ 
 $ M : matrix,  d, c' : vectors, r : scalar(integer) $ 
 $ c:=c'/r is a solutions of Sum[c_i M_i] = d where c_i is the i-th element $
 $ of the vector c and M_i is the i-th row vector of M.$
 $If there is no solution, then res-solv returns null. $
 (Note that M and d are not treated as an element of the homogenized Weyl)
 (algebra. If M or d contains the homogenization variable h, it automatically)
 (set to 1. If you need to use h, use the command res-solv-h)
 $[M d rng] res-solv [c' r] $ 
 $ rng is a ring object. $
 $ res-solv extracts variables names from rng, but defines a new ring. $
 $Example 1:  [(x,y) ring_of_differential_operators [[(x) -1 (Dx) 1]] weight_vector 0] $
 $              define_ring $
 $ [ [ [(x Dx + 2).] [ (Dx (x Dx + 3) - (x Dx + 2) (x Dx -4)).]]   [(1).]] $
 $  res-solv :: $
 $Example 2: $
 $ [ [ (x Dx + 2).  (Dx (x Dx + 3) - (x Dx + 2) (x Dx -4)).]   (1).] $
 $  res-solv :: $
 $Example 3: $
 $ [ [[(x Dx + 2). (0).] $
 $    [(Dx+3).     (x^3).]$
 $    [(3).        (x).]$
 $    [(Dx (x Dx + 3) - (x Dx + 2) (x Dx -4)). (0).]]   [(1). (0).]] $
 $  res-solv :: $
 $Example 4: $
 $ [[ (x*Dx+h^2). (Dx^2+x*h).] [(x^2+h^2). (h Dx + x^2).]] /ff set $
 $ [[ (x^2 Dx + x h^2). (Dx^3).]] /gg set  $
 $ [ff gg ff mul 0 get ] res-solv-h :: $
 $   $
 $res*solv and res*solv*h accept string inputs, too. For example,$
 $ [[ [ [(x Dx + 2)] [ (Dx (x Dx + 3) - (x Dx + 2) (x Dx -4))]]   [(1)]] $
 $  (x)]  res*solv :: $
]] putUsages
/res*solv {
  /arg1 set
  [/in-res*solv /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get res-toString expand res-solv /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

%% Solving \sum c_i M_i = d
%% [M d] res-solv-h c'/r  ;   
%% M : matrix,  d, c' : vectors, r : scalar, c'/r =c
/res-solv-h {
  /arg1 set
  [/in-res-solv-h /M /d /ans /B /vv /G /rr /rng /nn] pushVariables
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [
     /nn arg1 length def
     /M arg1  0 get def
     /d arg1  1 get def
     nn 3 eq {
       /rng arg1 2 get def
     }{
       M getRing /rng set
       rng tag RingP eq {  }
       { d getRing /rng set } ifelse
     }ifelse
     rng res-getx /vv set
     uli.verbose { (res-solv-h : vv = ) messagen vv message } { } ifelse
     uli.verbose { } { [(KanGBmessage) 0] system_variable } ifelse
     [vv from_records ring_of_differential_operators 0] define_ring
     M 0 get isArray {
       M { { toString . } map } map /M set
     } {
       M { toString . } map /M set
     } ifelse

     getOptions /options set
     (grade) (module1v) switch_function
     [M [(needBack)]] groebner /G set
     options restoreOptions

     G 1 get /B set

     d isArray {
       d 0 get isArray { [d] fromVectors 0 get /d set } { } ifelse
       [d] fromVectors 0 get /d set
     } {  } ifelse
     d toString . /d set

     /res-solv.d d def
     /res-solv.G G def

     d G 0 get reduction  /rr set
     rr 0 get (0). eq {
       [rr 2 get] B mul 0 get /ans set
       /ans [ ans { toString rng ,, (-1) rng ,, mul} map  
              rr 1 get toString .. ] def
     } {
       /ans null def
     } ifelse
     /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
/res*solv*h {
  /arg1 set
  [/in-res*solv*h /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get res-toString expand res-solv-h /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

%% See also xm, sm1_mul, sm1_mul_d, sm1_mul_h
/res*mul {
  /arg1 set
  [/in-res*mul /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get 0 get res-toString expand 
    A 0 get 1 get res-toString expand
    mul dehomogenize
    /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def
/res*mul*h {
  /arg1 set
  [/in-res*mul*h /A] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get 0 get res-toString expand 
    A 0 get 1 get res-toString expand
    mul 
    /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

%% cf. sm1_adjoint
/res*adjoint {
  /arg1 set
  [/in-res*adjoint /A /p /v /p0 /ans] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /A arg1 def
    A 1 get res-setRing pop
    A 0 get res-toString expand dehomogenize /p set
    /v res-setRing.v def
    p isArray {
      p { /p0 set [p0 v] res*adjoint } map /ans set
    }{
      p v adjoint dehomogenize /ans set
    }ifelse
    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
 
/res-init-m {
  /arg1 set
  [/in-res-init-m /A /ans] pushVariables
  [
    /A arg1 def
    A isArray {
       A { res-init-m } map /ans set
    }{
       A init /ans set
    }ifelse
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

/res-ord_w-m {
  /arg2 set
  /arg1 set
  [/in-ord_w-m /A /ans /w] pushVariables
  [
    /A arg1 def
    /w arg2 def
    A isArray {
       A { w res-ord_w-m } map /ans set
    }{
       A w ord_w /ans set
    }ifelse
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

%% cf. sm1_resol1
/res*resol1 {
  /arg1 set
  [/in-res*resol1 /A /ans /w /ans1 /ans2] pushVariables
  [
    /A arg1 def
    A length 3 ge {
     /w A 2 get def  %% weight vector
    } {
     /w null def
    }ifelse
    A resol1 /ans set
    /ans1 ans res-init-m def
    w tag 0 eq {
      /ans [ans ans1] def
    }{
      ans w 0 get res-ord_w-m /ans2 set
      /ans [ans ans1 ans2] def
    }ifelse
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

%% @@@

%% submodule to quotient module
%% M res-sub2Q  ==> J, where M \simeq D^m/J
/res-sub2Q {
  /arg1 set
  [/in-res-sub2Q /M /m] pushVariables
  [
    /M arg1 def
    M 0 get isArray {
    }{ M { [ 2 1 roll ] } map /M } ifelse
    M { dehomogenize homogenize } map /M set
    [M [(needSyz)]] groebner 2 get dehomogenize /arg1 set
  ] pop
  popVariables
  arg1
} def
[(res-sub2Q)
[(M res-sub2Q J)
 (matrix M, J; )
 (The submodule generated by M is isomorphic to D^m/J.)
]] putUsages


%% submodules to quotient module
%% [M N] res-subsub2Q  ==> J, where M \simeq D^m/J
/res-subsub2Q {
  /arg1 set
  [/in-res-subsub2Q /M /N /ss /m /n /ss2] pushVariables
  [
    /M arg1 0 get def
    /N arg1 1 get def
    /m M length def
    /n N length def
    M 0 get isArray {
    }{ M { [ 2 1 roll ] } map /M } ifelse
    N 0 get isArray {
    }{ N { [ 2 1 roll ] } map /N } ifelse
    M { dehomogenize homogenize } map /M set
    N { dehomogenize homogenize } map /N set
    [M N join [(needSyz)]] groebner 2 get /ss set
    ss dehomogenize /ss set
    ss { [ 2 1 roll  aload pop 1 1 n { pop pop } for ] } map
    /ss2 set
    ss2 {homogenize} map /ss2 set
    [ss2 0 get length [ss2] groebner 0 get dehomogenize ] toVectors2
    /arg1 set
  ] pop
  popVariables
  arg1
} def

/res-newpvec {
  /arg1 set
  [/in-res-newpvec /n ] pushVariables
  [
    /n arg1 def
    [1 1 n { pop (0). } for] /arg1 set
  ] pop
  popVariables
  arg1
} def

%% ki.sm1   kernel/image,  1999, 2/4
%% ki.sm1 is now moved to gbhg3/Int.
%% It is included in lib/complex.sm1
/kernel-image.v 1 def
/kernel-image.p 0 def % characteristic
%%
%%  D^p <-- m --- D^q <-- n -- D^r
%%       ker(m)/im(n)
%%
/res-kernel-image {
  /arg1 set
  [/in-res-kernel-image /p /q /r /m /n /t 
   /vlist  /s0 /s1 /ans
  ] pushVariables
  [
    /m arg1 0 get def
    /n arg1 1 get def
    /vlist arg1 2 get def
    vlist isArray {
      vlist from_records /vlist 
    } { } ifelse
    [vlist ring_of_differential_operators kernel-image.p] define_ring
    m { {toString . dehomogenize toString} map } map /m set  
    m length /q set
    n { {toString . dehomogenize toString} map } map /n set  
    n length /r set

    [m vlist] syz  0 get {{toString} map} map /s0 set
    /t s0 length def
    [ s0 n join vlist ] syz 0 get /s1 set 
    s1 { t carN } map /ans set

    /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(res-kernel-image)
[( [m n vlist] res-kernel-image c )
 (When, D^p <-- m --- D^q <-- n -- D^r )
 (D^q/c is isomorhic to ker(m)/im(n).)
 (vlist is a list of variables.)
]] putUsages


/res-dual {
  /arg1 set
  [/in-res-dual ] pushVariables
  [
    arg1 0 get /input set
    arg1 1 get /vlist set
    /n vlist length def
    /vv vlist from_records def

    %% preprocess to input resol0. Future version of resol1 should do them.
    input 0 get isArray {
      /kernel-image.unknowns input 0 get length def
    } { /kernel-image.unknowns 1 def } ifelse
    [vv ring_of_differential_operators 
     kernel-image.p ] define_ring
    input 0 get isArray {
       input { {toString . dehomogenize toString} map 
       } map /input set
    }{ input { toString . dehomogenize toString} map /input set } ifelse

    [input  vv]
    resol0 /rr set
    
    %% Postprocess of resol0
    [vv ring_of_differential_operators 
     kernel-image.p ] define_ring
    [ [kernel-image.unknowns rr 0 get { toString . dehomogenize } map]
       toVectors2 { {toString} map } map ]
    rr 1 get join /rr-syz set
    %%% end. The result is in rr-syz.

    /M rr-syz << n       >> get def
    /N rr-syz << n 1 sub >> get def
    M [ ] eq {
     /q N length def
     /M [ [0 1 q 1 sub { pop (0). } for] ] def
    } {  } ifelse

    %% regard them as a map from row vector v to row vector w; v M --> w
    uli.verbose {
      (M = ) messagen M pmat 
      (N = ) messagen N pmat
    } { } ifelse
    M transpose { { toString . dehomogenize vv adjoint} map } map /M set
    N transpose { { toString . dehomogenize vv adjoint} map } map /N set
    uli.verbose {
      $We are now computing ker (*N)/im (*M).$ message
      (*N = ) messagen N pmat
      (*M = ) messagen M pmat
      ( *N *M = ) messagen N M mul dehomogenize message
      (  ) message 
    }{  } ifelse
    /M M {{toString} map } map def
    /N N {{toString} map } map def
    [M N vv] res-kernel-image {{toString} map}map /ans1 set
    [ans1 vv] gb 0 get /arg1 set
  ] pop
  popVariables
  arg1
} def

[(res-dual)
[$[F V] res-dual G$
 $G is the dual D-module of F. V is a list of variables.$
 $Example 1:  [ [( x^3-y^2 )  ( 2 x Dx + 3 y Dy + 6 )  ( 2 y Dx + 3 x^2 Dy) ] $
 $              [(x) (y)]] res-dual $
 $Example 2:  [[1 3 4 5]] appell1 res-dual  $
 $Example 3:  [ [(-x1 Dx1 + x1 + 2) (x2 Dx2 - Dx2 -3)] [(x1) (x2)]] res-dual $
 $Example 4:  [ [(x2 Dx2 - Dx2 + 4) (x1 Dx1 + x1 +3)] [(x1) (x2)]] res-dual $
 $            3 and 4 are res-dual each other. $
 $Example 5:  [ [[1 1 1][0 1 2]] [0 0]] gkz res-dual $
 $Example 6:  [ [[1 1 1][0 1 2]] [-2 -1]] gkz res-dual $
 $    $
 $Example 7:  [ [(x Dx -1) (Dx^2)]     [(x)]] res-dual $
 $Example 8:  [ [[(1) (0)] [(0) (Dx)]] [(x)]] res-dual $
 $Example 9:  [ [((x Dx + x +1) (Dx-1))] [(x)]] res-dual $
]] putUsages

%%% From 1999/Int/sst.sm1
/saturation1 {
  /arg1 set
  [/in-saturation1 /ff /vlist /ulist /mm /hlist /iii
   /i  /uweight /aaa 
  ] pushVariables
  [(KanGBmessage) (CurrentRingp)] pushEnv
  [
    /ff arg1 def
    /iii ff 0 get {toString} map def  %% ideal
    /hlist ff 1 get {toString} map def %% saturation polynomials 
    /vlist [ff 2 get to_records pop] def
    /mm hlist length def

    [(KanGBmessage) 0] system_variable    
    /ulist [ 0 1 mm 1 sub { /i set [(_u) i] cat } for ] def
    /uweight ulist { 1 } map def
    [vlist ulist join from_records ring_of_polynomials
     [uweight] weight_vector 0] define_ring
    [0 1 mm 1 sub { /i set hlist i get .
                           ulist i get . mul (1). sub } for]
    /hlist set
    %%hlist pmat
    [iii {.} map hlist join] groebner_sugar 0 get /aaa set
    %%[aaa ulist] pmat
    aaa ulist eliminatev /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

[(saturation1)
[([ideal saturation-poly vlist] saturation jjj)
 $It returns(((ideal:f_1^\infty):f_2^\infty) ...) where$
 $saturation-poly is [f_1, f_2, ...]$
 $Example 1:   $
 $           [[(x1 y1 + x2 y2 + x3 y3 + x4 y4) $
 $             (x2 y2 + x4 y4) (x3 y3 + x4 y4) (y1 y4 - y2 y3)]$
 $            [(y1) (y2) (y3) (y4)] (x1,x2,x3,x4,y1,y2,y3,y4)] saturation1$
 $            /ff set [ff (x1,x2,x3,x4,y1,y2,y3,y4) $
 $                     [[(y1) 1 (y2) 1 (y3) 1 (y4) 1]]] pgb $
 $            0 get [(y1) (y2) (y3) (y4)] eliminatev ::$
]] putUsages


/intersection {
  /arg1 set
  [/in-intersection2 /ii /jj /rr /vlist /ii2 /jj2 ] pushVariables
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [
     /ii arg1 0 get def
     /jj arg1 1 get def
     /vlist arg1 2 get def
  
    [(KanGBmessage) 0] system_variable    

     [vlist to_records pop] /vlist set
     [vlist [(_t)] join from_records ring_of_differential_operators
      [[(_t) 1]] weight_vector 0] define_ring
     ii { toString . (_t). mul } map /ii2 set
     jj { toString . (1-_t). mul } map /jj2 set
     [ii2 jj2 join] groebner_sugar 0 get
     [(_t)] eliminatev /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

[(intersection)
[(Ideal intersections in the ring of differential operators.)
 $Example 1: [[[(x1) (x2)] [(x2) (x4)] (x1,x2,x3,x4)] intersection$
 $             [(x2) (x4^2)] (x1,x2,x3,x4)] intersection :: $
 $Example 2: [[[(x1) (x2)] [(x2) (x4)] (x1,x2,x3,x4)] intersection$
 $             [(x2) (x4^2)] (x1,x2,x3,x4)] intersection /ff set ff message$
 $           [ ff [(x2^2) (x3) (x4)] (x1,x2,x3,x4)] intersection :: $
 $Example 3: [[[(x1) (x2)] [(x2) (x4^2)] (x1,x2,x3,x4)] intersection$
 $             [(x2^2) (x3) (x4)] (x1,x2,x3,x4)] intersection :: $
]] putUsages


/saturation2 {
  /arg1 set
  [/in-saturation2 /ff /vlist /mm /slist /iii
   /i  /aaa 
  ] pushVariables
  [(KanGBmessage) (CurrentRingp)] pushEnv
  [
    /ff arg1 def
    /iii ff 0 get {toString} map def  %% ideal
    /slist ff 1 get {toString} map def %% saturation polynomials 
    /vlist ff 2 get  def
    /mm slist length def

    /aaa [iii [slist 0 get] vlist] saturation1 def
    1 1 mm 1 sub {
      /i set
      [[iii [slist i get] vlist] saturation1
       aaa vlist] intersection /aaa set
    } for
    /arg1 aaa def
  ] pop
  popEnv
  popVariables
  arg1
} def

[(saturation2)
[([ideal saturation-poly vlist] saturations jjj)
 $It returns (ideal:f_1^infty) \cap (ideal:f_2^\infty) \cap ... where$
 $saturation-poly is [f_1, f_2, ...]$
 $Example 1:   $
 $           [[(x1 y1 + x2 y2 + x3 y3 + x4 y4) $
 $             (x2 y2 + x4 y4) (x3 y3 + x4 y4) (y1 y4 - y2 y3)]$
 $            [(y1) (y2) (y3) (y4)] (x1,x2,x3,x4,y1,y2,y3,y4)] saturation2$
 $            /ff set [ff (x1,x2,x3,x4,y1,y2,y3,y4) $
 $                     [[(y1) 1 (y2) 1 (y3) 1 (y4) 1]]] pgb $
 $            0 get [(y1) (y2) (y3) (y4)] eliminatev ::$
 $Example 2: [[(x2^2) (x2 x4) (x2) (x4^2)] [(x2) (x4)] (x2,x4)] saturation2$
]] putUsages

/innerProduct {
  { [ 2 1 roll ] } map /innerProduct.tmp2 set
  /innerProduct.tmp1 set
  [innerProduct.tmp1] innerProduct.tmp2 mul 
  0 get 0 get
} def

/saturation {
  /arg1 set
  [/in-saturation /ff /vlist /mm /slist /iii
   /i  /aaa  /vlist2
  ] pushVariables
  [(KanGBmessage) (CurrentRingp)] pushEnv
  [
    /ff arg1 def
    /iii ff 0 get {toString} map def  %% ideal
    /slist ff 1 get {toString} map def %% saturation polynomials 
    /vlist ff 2 get  def
    /mm slist length def

    [vlist to_records pop] [(_z) (_y)] join /vlist2 set
    [vlist2 from_records ring_of_polynomials 
     [[(_z) 1 (_y) 1]] weight_vector
    0] define_ring 

    [
     [
      [0 1 mm 1 sub { /i set (_y). i npower } for ] 
      slist {.} map innerProduct  (_z). sub
     ]
     iii {.} map join

     [(_z)]
     vlist2 from_records
    ] saturation1 /aaa set

    [(KanGBmessage) 0] system_variable
    aaa {toString .} map /aaa set
    [aaa] groebner_sugar 0 get 
    [(_z) (_y)] eliminatev
    /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

[(saturation)
[([ideal J vlist] saturations jjj)
 $It returns (ideal : J^\infty) $
 (Saturation is computed in the ring of polynomials.) 
 $When J=[f_1, f_2, ...], it is equal to $
 $((ideal, z-(f_1 + y f_2 + y^2 f_3 +...)) : z^\infty) \cap k[x].$
 $Example 1:   $
 $           [[(x1 y1 + x2 y2 + x3 y3 + x4 y4) $
 $             (x2 y2 + x4 y4) (x3 y3 + x4 y4) (y1 y4 - y2 y3)]$
 $            [(y1) (y2) (y3) (y4)] (x1,x2,x3,x4,y1,y2,y3,y4)] saturation$
 $            /ff set [ff (x1,x2,x3,x4,y1,y2,y3,y4) $
 $                     [[(y1) 1 (y2) 1 (y3) 1 (y4) 1]]] pgb $
 $            0 get [(y1) (y2) (y3) (y4)] eliminatev ::$
 $Example 2: [[(x2^2) (x2 x4) (x2) (x4^2)] [(x2) (x4)] (x2,x4)] saturation$
]] putUsages


%% 2000, 6/8,  at Fernando Colon, 319,  Sevilla


/isExact.verbose 1 def  %% should be changed to gb.verbose
/isExact_h {
  /arg1 set
  [/in-isExact_h  /vv /comp /i /j /n /kernel.i /ans] pushVariables
  [
    /comp arg1 0 get def
    /vv arg1 1 get def
    /n comp length def
    /ans 1 def
    0 1 n 2 sub {
      /i set
      /j i 1 add def
      isExact.verbose { (Checking ker ) messagen i messagen ( = im of ) messagen
                     j message } {   } ifelse
      [comp i get vv] syz_h 0 get /kernel.i set
      [ kernel.i comp j get vv] isSameIdeal_h /ans set
      ans 0 eq {
        (image != kernel at ) messagen i messagen ( and ) messagen j message
         /LLL.isExact_h goto
      } {  } ifelse
      isExact.verbose { (OK) message } {  } ifelse
    } for
    /LLL.isExact_h
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

[(isExact_h)
[( complex isExact_h bool )
 (It returns 1 when the given complex is exact. All computations are done)
 (in D<h>, the ring of homogenized differential operators.)
 (cf. syz_h, isSameIdeal_h )
 $Example1: [ [[1 2 3]] [0]] gkz /ff set $
 $         [ff 0 get (x1,x2,x3) [[(x2) -1 (Dx2) 1]]] resol1 /gg set $
 $         [gg (x1,x2,x3)] isExact_h :: $
 $         gg 1 get 0 get /pp set $
 $         gg [1 1] pp put $
 $         [gg (x1,x2,x3)] isExact_h :: $
]] putUsages

/isExact {
  /arg1 set
  [/in-isExact  /vv /comp /i /j /n /kernel.i /ans] pushVariables
  [
    /comp arg1 0 get def
    /vv arg1 1 get def
    /n comp length def
    /ans 1 def
    0 1 n 2 sub {
      /i set
      /j i 1 add def
      isExact.verbose { (Checking ker ) messagen i messagen ( = im of ) messagen
                     j message } {   } ifelse
      [comp i get vv] syz 0 get /kernel.i set
      [ kernel.i comp j get vv] isSameIdeal /ans set
      ans 0 eq {
        (image != kernel at ) messagen i messagen ( and ) messagen j message
         /LLL.isExact goto
      } {  } ifelse
      isExact.verbose { (OK) message } {  } ifelse
    } for
    /LLL.isExact
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

[(isExact)
[( complex isExact bool )
 (It returns 1 when the given complex is exact. All computations are done)
 (in D, the ring of differentialoperators. Inputs are dehomogenized.)
 (cf. syz, isSameIdeal )
 $Example1: [ [[1 2 3]] [0]] gkz /ff set $
 $         [ff 0 get (x1,x2,x3) [[(x2) -1 (Dx2) 1]]] resol1 /gg set $
 $         [gg (x1,x2,x3)] isExact :: $
 $         gg 1 get 0 get /pp set $
 $         gg [1 1] pp put $
 $         [gg (x1,x2,x3)] isExact :: $
 $Example2: [ [[1 2 3]] [0]] gkz /ff set $
 $         [ff 0 get (x1,x2,x3) [[(x2) -1 (Dx2) 1]]] resol1 /gg set $
 $         gg dehomogenize /gg set $
 $         [gg (x1,x2,x3)] isExact :: $
 (       The syzygies of f_i^h in D<h> do not always give generators of )
 (       the corresponding syzygy of f_i in D.)
]] putUsages