% $OpenXM: OpenXM/src/kan96xx/Doc/hol.sm1,v 1.4 2000/03/14 13:01:28 takayama Exp $
%% hol.sm1, 1998, 11/8, 11/10, 11/14, 11/25, 1999, 5/18, 6/5.
%% rank, rrank, characteristic
%% This file is error clean.
/hol.version (2.990515) def
hol.version [(Version)] system_variable gt
{ [(This package hol.sm1 requires the latest version of kan/sm1) nl
   (Please get it from http://www.math.kobe-u.ac.jp/KAN) 
  ] cat
  error
} { } ifelse

$hol.sm1, basic package for holonomic systems (C) N.Takayama, 1999, 12/07 $ 
message-quiet

/rank.v [(x) (y) (z)] def   %% default value of v (variables).
/rank.ch [ ] def  %% characteristic variety.
/rank.verbose 0 def
/rank {
  /arg1 set
  [/in-rank /aa /typev /setarg /f /v /vsss /vddd
   /gg /wv /vd /vdweight /chv
   /one
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> rank) error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v rank.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    setarg { } { (rank : Argument mismatch) error } ifelse
    
    [(KanGBmessage) rank.verbose ] system_variable

    f { toString } map /f set  
    v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map  
    /vddd set   %% vddd = [(Dx) 1 (Dy) 1 (Dz) 1]
    v { @@@.Dsymbol 2 1 roll 2 cat_n } map 
    /vd set     %% vd = [(Dx) (Dy) (Dz)]
    /vdweight 
       vd { [ 2 1 roll -1 ] } map  %% vdweight=[[(Dx) -1] [(Dy) -1] [(Dz) -1]] 
    def   
    
    [v from_records
     ring_of_differential_operators [vddd] weight_vector 0] define_ring
    f { . dehomogenize } map /f set
    [f] groebner_sugar 0 get /gg set

    /wv vddd weightv def
    gg { wv init } map /chv set  %%obtained the characteristic variety.
    /rank.ch chv def
    chv { toString } map /chv set
   
    [ v vd join from_records
      ring_of_polynomials
      [vddd]  vdweight join weight_vector
      0
    ] define_ring
    [chv {.} map] groebner_sugar 0 get { init } map /chii set

    /rank.chii chii def
    rank.verbose { chii message } {  } ifelse
    v {[ 2 1 roll . (1).]} map /one set 
    %% [[(x). (1).] [(y). (1).] [(z). (1).]]
    %% chii { one replace } map  %% buggy code.
    %% Arg of hilb should be a reduced GB.
    [chii { one replace } map] groebner 0 get
    vd hilb /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def


[(rank)
 [( a rank b)
  ( array a;  number b)
  (Example 1 : )
  $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] rank :: $
  (Example 2 : )
  $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] rank :: $
 ]
] putUsages
(rank ) messagen-quiet

/characteristic.verbose 0 def
/characteristic.v [(x) (y) (z)] def
/characteristic.ch [ ] def
/ch { characteristic } def
/characteristic {
  /arg1 set
  [/in-rank /aa /typev /setarg /f /v /vsss /vddd
   /gg /wv /vd  /chv
   /one
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> characteristic) error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v characteristic.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    setarg { } { (rank : Argument mismatch) error } ifelse
    
    [(KanGBmessage) characteristic.verbose ] system_variable

    f { toString } map /f set  
    v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map  
    /vddd set   %% vddd = [(Dx) 1 (Dy) 1 (Dz) 1]
    v { @@@.Dsymbol 2 1 roll 2 cat_n } map 
    /vd set     %% vd = [(Dx) (Dy) (Dz)]
    
    [v from_records
     ring_of_differential_operators [vddd] weight_vector 0] define_ring
    f { . dehomogenize } map /f set
    [f] groebner_sugar 0 get /gg set

    /wv vddd weightv def
    gg { wv init } map /chv set
    /characteristic.ch [chv] def
%%    gg { wv init toString} map /chv set  %%obtained the characteristic variety.
%%    /characteristic.ch chv def
   
%%    [ v vd join from_records
%%      ring_of_polynomials
%%      [vddd] weight_vector
%%      0
%%    ] define_ring
%%    [chv {.} map] groebner_sugar 0 get /characteristic.ch set

    characteristic.ch /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

[(characteristic)
 [( a characteristic b)
  ( array a;  number b)
  (b is the generator of the characteristic variety of a.)
  (For the algorithm, see Japan J. of Industrial and Applied Math., 1994, 485--497.)
  (Example 1 : )
  $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] characteristic :: $
  (Example 2 : )
  $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] characteristic :: $
 ]
] putUsages
(characteristic ) messagen-quiet
[(ch)
 [(ch is the abbreviation of characteristic.)
  ( a ch b)
  ( array a;  number b)
  (b is the generator of the characteristic variety of a.)
  (For the algorithm, see, Japan J. of Industrial and Applied Math., 1994, 485--497.)
  (Example 1 : )
  $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] ch :: $
  (Example 2 : )
  $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] ch :: $
 ]
] putUsages
(ch ) messagen-quiet

%%%% developing rrank.sm1
/rrank.v [(x) (y) (z)] def   %% default value of v (variables).
/rrank.init [ ] def  %% initial ideal.
/rrank.verbose 0 def
/rrank {
  /arg1 set
  [/in-rrank /aa /typev /setarg /f /v /vsss /vddd
   /gg /wv /vd /vdweight
   /one /i /chv
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> rrank) error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v rrank.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    setarg { } { (rrank : Argument mismatch) error } ifelse
    
    [(KanGBmessage) rrank.verbose ] system_variable

    f { toString } map /f set  
    v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map  

    v { @@@.Dsymbol 2 1 roll 2 cat_n } map 
    /vd set     %% vd = [(Dx) (Dy) (Dz)] , v = [(x) (y) (z)]
    /vdweight 
      [ 0 1 v length 1 sub { /i set v i get << 0 i sub >>
                                    vd i get << i >> } for ]
    def   
    rrank.verbose { vdweight message } { } ifelse
    
    [v from_records
     ring_of_differential_operators [vdweight] weight_vector 0] define_ring
    f { . dehomogenize homogenize } map /f set
    [f] groebner 0 get {dehomogenize} map /gg set

    /wv vdweight weightv def
    gg { wv init } map /rrank.init set  %%obtained the initial ideal
    rrank.init {toString} map /chv set 
    /arg1 [chv v] rank def
  ] pop
  popEnv
  popVariables
  arg1
} def


[(rrank)
 [( a rrank b)
  ( array a;  number b)
  (It computes the holonomic rank for regular holonomic system.)
  (For the algorithm, see Grobner deformations of hypergeometric differential equations, 1999, Springer.)
  (Chapter 2.)
  (Example 1 : )
  $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] rrank :: $
 ]
] putUsages
(rrank ) messagen-quiet

/gb.v 1 def
/gb.verbose 0 def
/gb.options [ ] def
/gb {
  /arg1 set
  [/in-gb /aa /typev /setarg /f /v 
   /gg /wv /termorder /vec /ans /rr /mm
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> gb) error } ifelse
    /setarg 0 def
    /wv 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v gb.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse

    setarg { } { (gb : Argument mismatch) error } ifelse
    
    [(KanGBmessage) gb.verbose ] system_variable

    %%% Start of the preprocess
    f getRing /rr set
    %% To the normal form : matrix expression.
    f gb.toMatrixOfString /f set
    /mm gb.itWasMatrix def

    rr tag 0 eq {
      %% Define our own ring
      v isInteger {
        (Error in gb: Specify variables) error
      } {  } ifelse
      wv isInteger {        
        [v ring_of_differential_operators 
        0] define_ring
        /termorder 1 def
      }{
        [v ring_of_differential_operators 
         wv weight_vector
        0] define_ring
        wv gb.isTermOrder /termorder set
      } ifelse
    } {
      %% Use the ring structre given by the input.
      v isInteger not {
        (Warning : the given ring definition is not used.) message
      } {  } ifelse
      rr ring_def
      /wv rr gb.getWeight def
      wv gb.isTermOrder /termorder set
    } ifelse
    %%% Enf of the preprocess

    gb.verbose { (gb.options = ) messagen gb.options message } { } ifelse
    termorder {
      f { {. dehomogenize} map } map /f set
      [f gb.options] groebner_sugar 0 get /gg set
    }{
      f { {. dehomogenize} map} map /f set
      f fromVectors { homogenize } map /f set
      [f gb.options] groebner 0 get /gg set
    }ifelse
    wv isInteger {
      /ans [gg gg {init} map] def
    }{
      /ans [gg gg {wv 0 get weightv init} map] def
    }ifelse

    %% Postprocess : recover the matrix expression.
    mm {
      ans { /tmp set [mm tmp] toVectors } map
      /ans set
    }{ }
    ifelse
    %% 

    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
(gb ) messagen-quiet 

/pgb {
  /arg1 set
  [/in-pgb /aa /typev /setarg /f /v 
   /gg /wv /termorder /vec /ans /rr /mm
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage) (UseCriterion1)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { (<< array >> pgb) error } ifelse
    /setarg 0 def
    /wv 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v gb.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse

    setarg { } { (pgb : Argument mismatch) error } ifelse
    
    [(KanGBmessage) gb.verbose ] system_variable

    %%% Start of the preprocess
    f getRing /rr set
    %% To the normal form : matrix expression.
    f gb.toMatrixOfString /f set
    /mm gb.itWasMatrix def

    rr tag 0 eq {
      %% Define our own ring
      v isInteger {
        (Error in pgb: Specify variables) error
      } {  } ifelse
      wv isInteger {        
        [v ring_of_polynomials
        0] define_ring
        /termorder 1 def
      }{
        [v ring_of_polynomials
         wv weight_vector
        0] define_ring
        wv gb.isTermOrder /termorder set
      } ifelse
    } {
      %% Use the ring structre given by the input.
      v isInteger not {
        (Warning : the given ring definition is not used.) message
      } {  } ifelse
      rr ring_def
      /wv rr gb.getWeight def
      wv gb.isTermOrder /termorder set
    } ifelse
    %%% Enf of the preprocess

    gb.verbose { (gb.options = ) messagen gb.options message } { } ifelse
    termorder {
      f { {. dehomogenize} map } map /f set
      [(UseCriterion1) 1] system_variable
      [f gb.options] groebner_sugar 0 get /gg set
      [(UseCriterion1) 0] system_variable
    }{
      f { {. dehomogenize} map} map /f set
      f fromVectors { homogenize } map /f set
      [(UseCriterion1) 1] system_variable
      [f gb.options] groebner 0 get /gg set
      [(UseCriterion1) 0] system_variable
    }ifelse
    wv isInteger {
      /ans [gg gg {init} map] def
    }{
      /ans [gg gg {wv 0 get weightv init} map] def
    }ifelse

    %% Postprocess : recover the matrix expression.
    mm {
      ans { /tmp set [mm tmp] toVectors } map
      /ans set
    }{ }
    ifelse
    %% 

    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def

/pgb.old {
  /arg1 set
  [/in-pgb /aa /typev /setarg /f /v 
   /gg /wv /termorder /vec /ans
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage) (UseCriterion1)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { (array pgb) message (pgb) usage error } ifelse
    /setarg 0 def
    /wv 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v gb.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse

    setarg { } { (pgb : Argument mismatch) message error } ifelse
    
    [(KanGBmessage) gb.verbose ] system_variable

    %% Input must not be vectors.
    f { toString } map /f set  

    wv isInteger {        
      [v ring_of_polynomials
      0] define_ring
      /termorder 1 def
    }{
      [v ring_of_polynomials
       wv weight_vector
      0] define_ring
      wv gb.isTermOrder /termorder set
    } ifelse
    termorder {
      f { . dehomogenize } map /f set
      [(UseCriterion1) 1] system_variable
      [f] groebner_sugar 0 get /gg set
      [(UseCriterion1) 0] system_variable
    }{
      f { . dehomogenize homogenize} map /f set
      [(UseCriterion1) 1] system_variable
      [f] groebner 0 get /gg set
      [(UseCriterion1) 0] system_variable
    }ifelse
    wv isInteger {
      /ans [gg gg {init} map] def
    }{
      /ans [gg gg {wv 0 get weightv init} map] def
    }ifelse
    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1
} def
(pgb ) messagen-quiet 

/gb.toMatrixOfString {
  /arg1 set
  [/in-gb.toMatrixOfString /ff /aa /ans] pushVariables
  [
     /aa arg1 def
     aa length 0 eq { /ans [ ] def /gb.toMatrixOfString.LLL goto }{ } ifelse
     aa 0 get isArray {
       /gb.itWasMatrix aa 0 get length def
     }{
       /gb.itWasMatrix 0 def 
     } ifelse
     aa {
       /ff set
       ff isArray {
         ff {toString} map /ff set
       }{
         [ff toString] /ff set
       } ifelse
       ff
     } map /ans set
    /gb.toMatrixOfString.LLL
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(gb.toMatrixOfString)
[(It translates given input into a matrix form which is a data structure)
 (for computations of kernel, image, cokernel, etc.)
 (gb.itWasMatrix is set to the length of the input vector.)
 $Example 1: $
 $  [ (x). (y).] gb.toMatrixOfString ==> [[(x)] [(y)]] $
 $  gb.itWasMatrix is 0.$
 $Example 2: $
 $  [ [(x). (1).] [(y). (0).]] gb.toMatrixOfString ==>  [ [(x) (1)] [(y) (0)]] $ 
 $  gb.itWasMatrix is 2.$
]] putUsages

/gb.toMatrixOfPoly {
  /arg1 set
  [/in-gb.toMatrixOfPoly /ff /aa /ans] pushVariables
  [
     /aa arg1 def
     aa length 0 eq { /ans [ ] def /gb.toMatrixOfPoly.LLL goto }{ } ifelse
     aa 0 get isArray {
       /gb.itWasMatrix aa 0 get length def
     }{
       /gb.itWasMatrix 0 def 
     } ifelse
     aa {
       /ff set
       ff isArray {
       }{
         [ff] /ff set
       } ifelse
       ff
     } map /ans set
    /gb.toMatrixOfPoly.LLL
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(gb.toMatrixOfPoly)
[(It translates given input into a matrix form which is a data structure)
 (for computations of kernel, image, cokernel, etc.)
 (gb.itWasMatrix is set to the length of the input vector.)
 $Example 1: $
 $  [ (x). (y).] gb.toMatrixOfPoly ==> [[(x)] [(y)]] $
 $  gb.itWasMatrix is 0.$
 $Example 2: $
 $  [ [(x). (1).] [(y). (0).]] gb.toMatrixOfPoly ==>  [ [(x) (1)] [(y) (0)]] $ 
 $  gb.itWasMatrix is 2.$
]] putUsages

/gb.getWeight {
  /arg1 set
  [/in-gb.getWeight /rr /ww /vv /ans /nn /ii] pushVariables
  [(CurrentRingp)] pushEnv
  [
     /rr arg1 def
     rr ring_def
     getVariableNames /vv set
     [(orderMatrix)] system_variable 0 get /ww set
     /nn vv length 1 sub def
     [0 1 nn {
        /ii set
        ww ii get 0 eq {
        } {
          vv ii get
          ww ii get
        } ifelse
      } for
     ] /ans set
     /arg1 [ans] def
  ] pop
  popEnv
  popVariables
  arg1
} def
[(gb.getWeight)
[(ring gb.getWeight wv)
 (It gets the weight vector field of the ring ring.)
]] putUsages


/gb.isTermOrder {
  /arg1 set
  [/in-gb.isTermOrder /vv /ww /yes /i /j] pushVariables
  [
     /vv arg1 def
     /yes 1 def
     0 1 vv length 1 sub {
       /i set
       /ww vv i get def
       0 1 ww length 1 sub {
          /j set
          ww j get isInteger {
            ww j get 0 lt { /yes 0 def } { } ifelse
          }{ } ifelse
       }for
     }for
     /arg1 yes def
  ] pop
  popVariables
  arg1
} def
[(gb)
 [(a gb b)
  (array a; array b;)
  (b : [g ii];  array g; array in; g is a Grobner basis of f)
  (             in the ring of differential operators.)
   $            ii is the initial ideal in case of w is given or <<a>> belongs$
   $            to a ring. In the other cases, it returns the initial monominal.$
  (a : [f ];    array f;  f is a set of generators of an ideal in a ring.)
  (a : [f v];   array f; string v;  v is the variables. )
  (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
  (  )
  $Example 1: [ [( (x Dx)^2 + (y Dy)^2 -1) ( x y Dx Dy -1)] (x,y) $
  $             [ [ (Dx) 1 ] ] ] gb pmat ; $
  (Example 2: )
  (To put h=1, type in, e.g., )
  $ [ [(2 x Dx + 3 y Dy+6) (2 y Dx + 3 x^2 Dy)] (x,y) $
  $   [[(x) -1 (Dx) 1 (y) -1 (Dy) 1]]] gb /gg set gg dehomogenize pmat ;$
  (  )
  $Example 3: [ [( (x Dx)^2 + (y Dy)^2 -1) (  x y Dx Dy -1)] (x,y) $
  $             [ [ (Dx) 1 (Dy) 1] ] ] gb pmat ; $
  (  )
  $Example 4: [[ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] (x,y) $
  $             [ [ (x) -1 (y) -1] ] ] gb pmat ; $
  (  )
  (cf. gb, groebner, groebner_sugar, syz. )
]] putUsages

[(pgb)
 [(a pgb b)
  (array a; array b;)
  (b : [g ii];  array g; array in; g is a Grobner basis of f)
  (             in the ring of polynomials.)
  $             ii is the initial ideal in case of w is given or <<a>>belongs$
  $             to a ring. In the other cases, it returns the initial monominal.$
  (a : [f ];    array f;  f is a set of generators of an ideal in a ring.)
  (a : [f v];   array f; string v;  v is the variables.)
  (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
  $Example 1: [(x,y) ring_of_polynomials 0] define_ring $
  $           [ [(x^2+y^2-4). (x y -1).] ] pgb :: $
  $Example 2: [ [(x^2+y^2) (x y)]   (x,y)  [ [(x) -1 (y) -1] ] ] pgb :: $
  (cf. gb, groebner, groebner_sugar, syz. )
]] putUsages


%/syz.v 1 def
/syz.v 1 def
/syz.verbose 0 def
/syz {
  /arg1 set
  [/in-syz /aa /typev /setarg /f /v 
   /gg /wv /termorder /vec /ans /ggall /vectorInput /vsize /gtmp /gtmp2
   /rr /mm
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { (<< array >> syz) error } ifelse
    /setarg 0 def
    /wv 0 def
    aa { tag } map /typev set
    typev [ ArrayP ] eq
    {  /f aa 0 get def
       /v syz.v def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse
    typev [ArrayP ArrayP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get from_records def
       /wv aa 2 get def
       /setarg 1 def
    } { } ifelse

    setarg { } { (syz : Argument mismatch) error } ifelse
    
    [(KanGBmessage) syz.verbose ] system_variable



    %%% Start of the preprocess
    f getRing /rr set
    %% To the normal form : matrix expression.
    f gb.toMatrixOfString /f set
    /mm gb.itWasMatrix def
    mm 0 gt {
      /vectorInput 1 def
    }{
      /vectorInput 1 def
    } ifelse

    rr tag 0 eq {
      %% Define our own ring
      v isInteger {
        (Error in syz: Specify variables) error
      } {  } ifelse
      wv isInteger {        
        [v ring_of_differential_operators 
        0] define_ring
        /termorder 1 def
      }{
        [v ring_of_differential_operators 
         wv weight_vector
        0] define_ring
        wv gb.isTermOrder /termorder set
      } ifelse
    }{
      %% Use the ring structre given by the input.
      v isInteger not {
        (Warning : the given ring definition is not used.) message
      } {  } ifelse
      rr ring_def
      /wv rr gb.getWeight def
      wv gb.isTermOrder /termorder set
    } ifelse
    %%% Enf of the preprocess

    termorder {
      f { {. dehomogenize} map } map /f set
      [f [(needBack) (needSyz)]] groebner_sugar /ggall set
      ggall 2 get /gg set
    }{
      f { {. dehomogenize } map homogenize } map /f set
      [f [(needBack) (needSyz)]] groebner /ggall set 
      ggall 2 get /gg set
    }ifelse
    vectorInput {
      /vsize f 0 get length def  %% input vector size.
      /gtmp ggall 0 get def
       [vsize gtmp] toVectors /gtmp set
       ggall 0 gtmp put
    }{  } ifelse
      /arg1 [gg dehomogenize ggall] def
  ] pop
  popEnv
  popVariables
  arg1
} def
(syz ) messagen-quiet 

[(syz)
 [(a syz [b c])
  (array a; array b; array c)
  (b is a set of generators of the syzygies of f.)
  (c = [gb, backward transformation, syzygy without dehomogenization].)
  (See groebner.)
  (a : [f ];    array f;  f is a set of generators of an ideal in a ring.)
  (a : [f v];   array f; string v;  v is the variables.)
  (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
  $Example 1: [(x,y) ring_of_polynomials 0] define_ring $
  $           [ [(x^2+y^2-4). (x y -1).] ] syz :: $
  $Example 2: [ [(x^2+y^2) (x y)]   (x,y)  [ [(x) -1 (y) -1] ] ] syz :: $
  $Example 3: [ [( (x Dx)^2 + (y Dy)^2 -1) ( x y Dx Dy -1)] (x,y) $
  $             [ [ (Dx) 1 ] ] ] syz pmat ; $
  $Example 4:  [ [(2 x Dx + 3 y Dy+6) (2 y Dx + 3 x^2 Dy)] (x,y) $
  $             [[(x) -1 (Dx) 1 (y) -1 (Dy) 1]]] syz pmat ;$
  $Example 5:  [ [ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] $
  $              (x,y) ] syz pmat ;$
  $Example 6:  [ [ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] $
  $              (x,y) [[(x) -1 (y) -2]] ] syz pmat ;$
  $Example 7:  [ [ [(0) (0)] [(0) (0)] [(x) (y)]] $
  $              [(x) (y)]] syz pmat ;$
]] putUsages


%%%%%%%%%%%%%%%%%% package fs  %%%%%%%%%%%%%%%%%%%%%%%
[(genericAnn)
 [ (f [s v1 v2 ... vn] genericAnn [L1 ... Lm])
   (L1, ..., Lm are annihilating ideal for f^s.)
   (f is a polynomial of v1, ..., vn)
   (<string> | <poly>  f, s, v1, ..., vn ; <poly>  L1, ..., Lm )
   $Example: (x^3+y^3+z^3) [(s) (x) (y) (z)] genericAnn$
 ]
] putUsages ( genericAnn ) messagen-quiet
/fs.verbose 0 def
/genericAnn {
  /arg2 set /arg1 set
  [/in-genericAnn /f /vlist /s  /vvv /nnn /rrr
   /v1  /ops /ggg /ggg0 
   ] pushVariables
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [
    /f arg1 def  /vlist arg2 def
    f toString /f set
    vlist { toString } map /vlist set
    [(KanGBmessage) fs.verbose] system_variable
    /s vlist 0 get def
    /vvv (_u,_v,_t,) vlist rest { (,) 2 cat_n } map aload length /nnn set
         s nnn 2 add cat_n def
    fs.verbose { vvv message } {  }ifelse
    [vvv ring_of_differential_operators
     [[(_u) 1 (_v) 1]] weight_vector 0] define_ring /rrr set
   
    [ (_u*_t). f . sub  (_u*_v-1). ]
    vlist rest { /v1 set  
%%D-clean   f . (D) v1 2 cat_n . 1 diff0 (_v*D_t). mul
        f . @@@.Dsymbol v1 2 cat_n . 1 diff0 [(_v*) @@@.Dsymbol (_t)] cat . mul
        @@@.Dsymbol v1 2 cat_n . add } map
    join
    /ops set
    ops {[[(h). (1).]] replace } map /ops set
    fs.verbose { ops message  } {  }ifelse
    [ops] groebner_sugar 0 get /ggg0 set
    fs.verbose { ggg0 message } { } ifelse
    ggg0 [(_u) (_v)] eliminatev
%%D-clean        { [(_t).] [ (D_t).] [s .] distraction 
        { [(_t).] [ [@@@.Dsymbol (_t)] cat .] [s .] distraction 
          [[s . << (0). s . sub (1). sub >>]] replace
         } map /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

%% Find differential equations for  f^(m), r0 the minimal integral root.
[(annfs)
 [( [ f v m r0] annfs g )
  (It returns the annihilating ideal of f^m where r0 must be smaller)
  (or equal to the minimal integral root of the b-function.)
  (Or, it returns the annihilating ideal of f^r0, r0 and the b-function)
  (where r0 is the minial integral root of b.)
  (For the algorithm, see J. Pure and Applied Algebra 117&118(1997), 495--518.)
  (Example 1: [(x^2+y^2+z^2+t^2) (x,y,z,t) -1 -2] annfs :: )
  $           It returns the annihilating ideal of (x^2+y^2+z^2+t^2)^(-1).$
  (Example 2: [(x^2+y^2+z^2+t^2) (x,y,z,t)] annfs :: )
  $           It returns the annihilating ideal of f^r0 and [r0, b-function]$
  $           where r0 is the minimal integral root of the b-function.$
  (Example 3: [(x^2+y^2+z^2) (x,y,z) -1 -1] annfs :: )
  (Example 4: [(x^3+y^3+z^3) (x,y,z)] annfs :: )
  (Example 5: [((x1+x2+x3)(x1 x2 + x2 x3 + x1 x3) - t x1 x2 x3 ) )
  (            (t,x1,x2,x3) -1 -2] annfs :: )
  (           Note that the example 4 uses huge memory space.)
]] putUsages
( annfs ) messagen-quiet
/annfs.verbose fs.verbose def
/annfs.v [(x) (y) (z)] def
/annfs.s (_s) def
%% The first variable must be s.
/annfs {
  /arg1 set
  [/in-annfs /aa /typev /setarg  /v /m /r0 /gg /ss /fs /gg2
   /ans /vtmp /w2 /velim /bbb /rrr /r0
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> annfs) error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    /r0 [ ] def 
    /m  [ ]  def
    /v annfs.v def
    aa 0 << aa 0 get toString >> put
    typev [ StringP ] eq
    {  /f aa 0 get def
       /setarg 1 def
    } { } ifelse
    typev [StringP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /setarg 1 def
    } { } ifelse
    typev [StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    typev [StringP ArrayP IntegerP IntegerP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /m aa 2 get def
       /r0 aa 3 get def
       /setarg 1 def
    } { } ifelse
    typev [StringP StringP IntegerP IntegerP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /m aa 2 get def
       /r0 aa 3 get def
       /setarg 1 def
    } { } ifelse
    setarg 1 eq { } { (annfs : wrong argument) error } ifelse

    [annfs.s] v join /v set

    /ss v 0 get def
    annfs.verbose {
       (f, v, s, f^{m}, m+r0 = ) messagen
       [ f  (, ) v  (, )  ss  (, ) 
         (f^) m (,) m (+)  r0 ] {messagen} map ( ) message
    } { } ifelse

    f v genericAnn /fs set

    annfs.verbose {
      (genericAnn is ) messagen fs message
    } { } ifelse
    [(KanGBmessage) annfs.verbose] system_variable

    m isArray {
      %% Now, let us find  the b-function. /vtmp /w2 /velim /bbb /rrr /r0
      v rest { /vtmp set vtmp  @@@.Dsymbol vtmp 2 cat_n } map /velim set
      velim { 1 } map /w2 set
      annfs.verbose { w2 message } {  } ifelse
      [v from_records ring_of_differential_operators 
       [w2] weight_vector 0] define_ring
      [ fs { toString . } map [ f toString . ] join ]
      groebner_sugar 0 get velim eliminatev 0 get /bbb set
      [[(s) annfs.s] from_records ring_of_polynomials 0] define_ring
      bbb toString . [[annfs.s . (s).]] replace /bbb set
      annfs.verbose { bbb message } {  } ifelse
      bbb findIntegralRoots  /rrr set
      rrr 0 get /r0 set  %% minimal integral root.
      annfs.verbose { rrr message } {  } ifelse
      fs 0 get (ring) dc ring_def
      fs { [[annfs.s . r0 toString .]] replace } map /ans set
      /ans [ans [r0 bbb]] def
      /annfs.label1 goto
    } { } ifelse
    m 0 ge {
      (annfs works only for getting annihilating ideal for f^(negative))
      error
    } { } ifelse
    r0 isArray {
      [(Need to compute the minimal root of b-function) nl
       (It has not been implemented.) ] cat
      error
    } {  } ifelse

    [v from_records ring_of_differential_operators 0] define_ring
    fs {toString . dehomogenize [[ss . r0 (poly) dc]] replace} 
       map /gg set
    annfs.verbose { gg message } { } ifelse

    [ [f . << m r0 sub >> npower ] gg join  
      [(needBack) (needSyz)]] groebner_sugar 2 get /gg2 set
    
    gg2 { 0 get } map /ans set
    /ans ans { dup (0). eq {pop} { } ifelse } map def

    /annfs.label1
    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1 
} def
  
/genericAnnWithL.s (s) def
/annfs.verify 0 def
/genericAnnWithL {
  /arg1 set
  [/in-genericAnnWithL /aa /typev /setarg  /v /m /r0 /gg /ss /fs /gg2
   /ans /vtmp /w2 /velim /bbb /rrr /r0  /myL /mygb /jj
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage) (Homogenize)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> annfs) error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    /r0 [ ] def 
    /m  [ ]  def
    /v annfs.v def
    aa 0 << aa 0 get toString >> put
    typev [ StringP ] eq
    {  /f aa 0 get def
       /setarg 1 def
    } { } ifelse
    typev [StringP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /setarg 1 def
    } { } ifelse
    typev [StringP ArrayP] eq
    {  /f aa 0 get def
       /v aa 1 get def
       /setarg 1 def
    } { } ifelse
    setarg 1 eq { } { (genericAnnWithL : wrong argument) error } ifelse

    [genericAnnWithL.s] v join /v set

    /ss v 0 get def
    annfs.verbose {
       (f, v, s, f^{m}, m+r0 = ) messagen
       [ f  (, ) v  (, )  ss  (, ) 
         (f^) m (,) m (+)  r0 ] {messagen} map ( ) message
    } { } ifelse

    f v genericAnn /fs set

    annfs.verbose {
      (genericAnn is ) messagen fs message
    } { } ifelse
    [(KanGBmessage) annfs.verbose] system_variable

    m isArray {
      %% Now, let us find  the b-function. /vtmp /w2 /velim /bbb /rrr /r0
      v rest { /vtmp set vtmp  @@@.Dsymbol vtmp 2 cat_n } map /velim set
      velim { 1 } map /w2 set
      annfs.verbose { w2 message } {  } ifelse
      [v from_records ring_of_differential_operators 
       [w2] weight_vector 0] define_ring

      [ [ f toString . ] fs { toString . } map join [(needBack)]]
      groebner_sugar /mygb set
      mygb 0 get velim eliminatev 0 get /bbb set
      mygb 0 get bbb position /jj set
      mygb 1 get jj get 0 get /myL set

      annfs.verbose { bbb message } {  } ifelse

      annfs.verify {
        (Verifying L f - b belongs to genericAnn(f)) message
        [(Homogenize) 0] system_variable
        << myL f . mul bbb sub >> 
        [fs { toString . } map] groebner_sugar 0 get
        reduction 0 get message
        (Is it zero? Then it's fine.) message
      } { } ifelse

      /ans [bbb [myL fs] ] def
      /annfs.label1 goto
    } { } ifelse

    /annfs.label1
    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1 
} def


[(genericAnnWithL)
[$[f v] genericAnnWithL [b [L I]]$
 $String f,v; poly b,L; array of poly I;$
 $f is a polynomial given by a string. v is the variables.$
 $ v must not contain names  s, e.$
 $b is the b-function (Bernstein-Sato polynomial) for f and$
 $ L is the operator satisfying L f^{s+1} = b(s) f^s $
 $ I is the annihilating ideal of f^s.$
 $cf. bfunction, annfs, genericAnn.$
 $Example 1:  [(x^2+y^2) (x,y)] genericAnnWithL ::$
 $Example 2:  [(x^2+y^2+z^2) (x,y,z)] genericAnnWithL ::$
 $Example 3:  [(x^3-y^2 z^2) (x,y,z)] genericAnnWithL ::$
]] putUsages

/reduction*.noH 0 def
/reduction* {
  /arg1 set
  [/in-reduction* /aa /typev /setarg /f /v 
   /gg /wv /termorder /vec /ans /rr /mm /h /size /a0 /a3
   /opt
  ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv
  [

    /aa arg1 def
    aa isArray { } { ( << array >> reduction*) error } ifelse
    /setarg 0 def
    /wv 0 def
    aa { tag } map /typev set
    typev [StringP ArrayP ArrayP] eq
    typev [ArrayP ArrayP ArrayP] eq or
    typev [PolyP ArrayP ArrayP] eq or
    {  /h aa 0 get def
       /f aa 1 get def
       /v aa 2 get from_records def
       /setarg 1 def
    } { } ifelse
    typev [StringP ArrayP ArrayP ArrayP] eq
    typev [ArrayP ArrayP ArrayP ArrayP] eq or
    typev [PolyP ArrayP ArrayP ArrayP] eq or
    {  /h aa 0 get def
       /f aa 1 get def
       /v aa 2 get from_records def
       /wv aa 3 get def
       /setarg 1 def
    } { } ifelse

    setarg { } { (reduction* : Argument mismatch) error } ifelse
    
    [(KanGBmessage) gb.verbose ] system_variable

    %%% Start of the preprocess
    f getRing /rr set
    

    rr tag 0 eq {
      %% Define our own ring
      v isInteger {
        (Error in reduction*: Specify variables) error
      } {  } ifelse
      wv isInteger {        
        [v ring_of_differential_operators 
        0] define_ring
        /termorder 1 def
      }{
        [v ring_of_differential_operators 
         wv weight_vector
        0] define_ring
        wv gb.isTermOrder /termorder set
      } ifelse
    } {
      %% Use the ring structre given by the input.
      v isInteger not {
        (Warning : the given ring definition is not used.) message
      } {  } ifelse
      rr ring_def
      /wv rr gb.getWeight def
      wv gb.isTermOrder /termorder set
    } ifelse
    %%% Enf of the preprocess

    f 0 get isArray {
      /size f 0 get length def
      f { { toString . } map } map /f set
      f fromVectors /f set
    }{
      /size -1 def
      f { toString . } map /f set
    } ifelse 

    h isArray {
      h { toString . } map /h set
      [h] fromVectors 0 get /h set
    }{ 
      h toString . /h set
    } ifelse
    f { toString . } map /f set
    getOptions /opt set
    [(ReduceLowerTerms) 1] system_variable
    reduction*.noH {
      h f reduction-noH /ans set
    } {
      h f reduction /ans set
    } ifelse
    opt restoreOptions
    size -1 eq not {
      [size ans 0 get] toVectors /a0 set
      [size ans 3 get] toVectors /a3 set
      /ans [a0 ans 1 get ans 2 get a3] def
    } { } ifelse
    /arg1 ans def
  ] pop
  popEnv
  popVariables
  arg1  
} def 


[(reduction*)
[([f base v] reduction* [h c0 syz input])
 ([f base v weight] reduction* [h c0 syz input])
 (reduction* is an user interface for reduction and reduction-noH.)
 (If reduction*.noH is one, then reduction-noH will be called.)
 (Example 1: [(x^2) [(x^2+y^2-4) (x y-1)] [(x) (y)]] reduction* )
 (Example 2: [[(1) (y^2-1)] [ [(0) (y-1)] [(1) (y+1)]] [(x) (y)]] reduction*)
 (Example 3: [(x^2) [(x^2+y^2-4) (x y-1)] [(x) (y)] [[(x) 10]] ] reduction* )
]] putUsages

( ) message-quiet ;