% $OpenXM: OpenXM/src/kan96xx/Doc/slope.sm1,v 1.4 2012/06/11 05:23:52 takayama Exp $
(oxasir.sm1.loaded) boundp not {
   [(parse) (oxasir.sm1)  pushfile] extension
} { } ifelse
(cohom.sm1.loaded) boundp not {
   [(parse) (cohom.sm1)  pushfile] extension
} { } ifelse
$slope.sm1, computing the slopes of a D-ideal:  June 15, 2000$ message
$                            (C) N.Takayama, F.Castro-Jimenez$ message
$Imported commands:  slope $ message
/slope.verbose 1 def
/gb.warning 0 def
/slope.geometric 1 def %%Computing the geometric slope. Load cohom.sm1 and oxasir.

/slope.infinity (99999999999999999999).. def
/w_support {
  /arg2 set
  /arg1 set
  [/in-w_support /f /wvec  /ans /g /tt] pushVariables
  [
     /f  arg1 def
     /wvec arg2 def
     /ans [ ] def
     {
        f (0). eq { exit } { } ifelse
        f init /g set
        wvec { g 2 1 roll ord_w (universalNumber) dc } map /tt set
        ans tt append /ans set
        f g sub /f set
     } loop
     /arg1 ans def
  ] pop
  popVariables
  arg1
} def

[(w_support)
 [$f [w1 w2 ...] w_support [ [i1 i2 ...] [j1 j2 ...] [k1 k2 ...] ...]$
  $ i1, ..., j1, ..., k1, ... are universal numbers. $
  $Example: (x Dx+ x ). [ [(x) -1 (Dx) 1] [(Dx) 1]] w_support$
 ]
] putUsages


/w_supports_of_I {
  /arg1 set  
  [/in-w_supports_of_I /ans /v /ff /wvec /gg /gg2] pushVariables
  [
     /ff arg1 0 get def
     /v arg1 1 get def
     /wvec arg1 2 get def
     wvec { [ 2 1 roll ] [ ff v 4 -1 roll ] gb } map /gg set
     gg { 0 get } map /gg set
     gg flatten /gg2 set
     gg2 message
     gg2 0 get (ring) dc ring_def
     gg2 { (string) dc  . } map /gg2 set % reparse
     gg2 { wvec w_support } map /ans set
     /arg1 [ans gg] def
  ] pop
  popVariables
  arg1
} def

[(w_supports_of_I)
 [$[f v [w1 w2 ...]] w_support_of_I [supports gb]$
  $Example 1: [[(x Dx + 2 y Dy) (Dx^2-Dy)] (x,y) [  [(Dx) 1 (Dy) 1] [(y) -1 (Dy) 1]]]$
  $          w_supports_of_I$
  $Example 2: [ [[1 2 3]] [0]] gkz /ff set$
  $          [ ff 0 get ff 1 get [ [(Dx1) 1 (Dx2) 1 (Dx3) 1] [(Dx3) 1 (x3) -1]]$
  $          ] w_supports_of_I $
  $Example 3: [ [[1 2 3]] [0]] gkz /ff set$
  $         [ ff 0 get ff 1 get [ [(x1) 0 (x2) 0 (x3) -3 (Dx1) 6 (Dx2) 6 (Dx3) 9]]] gb /gg set $
  $ gg 1 get { [ [(Dx1) 1 (Dx2) 1 (Dx3) 1] [(Dx3) 1 (x3) -1]] w_support } map /gg2 set $
 ]
] putUsages

/w_supports_of_I_without_gb_computation {
  /arg1 set  
  [/in-w_supports_of_I_without_gb_computation
    /ans /v /ff /wvec /gg2] pushVariables
  [
     /ff arg1 0 get def
     /v arg1 1 get def
     /wvec arg1 2 get def
     /gg2 ff def
%%     gg2 message
     gg2 0 get (ring) dc ring_def
     gg2 { (string) dc  . } map /gg2 set % reparse
     gg2 { wvec w_support } map /ans set
     /arg1 [ans gg2] def
  ] pop
  popVariables
  arg1
} def

/decompose_to_w_homogeneous {
  /arg1  set
  [/in-decompose_to_w_homogeneous /f /w /g /ans] pushVariables
  [
     /f arg1 0 get def
     /w arg1 1 get def
     /ans [ ] def
     f (ring) dc ring_def
     /w w weightv def
     {
       f (0). eq { exit } {  } ifelse
       f w init /g set
       ans g append /ans set
       f g sub /f set
     } loop
     /arg1 ans def
  ] pop
  popVariables
  arg1
} def

[(decompose_to_w_homogeneous)
[( [f w] decompose_to_w_homogeneous [f0 f1 f2 ...])
 $Example: [ (x^3+x*h^4+x+1). [(x) 2 (h) 1] ] decompose_to_w_homogeneous $
]] putUsages

%% Check in the polynomial ring.
/w_homogeneousQ {
  /arg1 set
  [/in-w_homogeneousQ  /ii /vv /ww /ans /gg /jj /i] pushVariables
  [
    /ii arg1 0 get def
    /vv arg1 1 get def
    /ww arg1 2 get def
    [ii vv] pgb 0 get /gg set
    gg 0 get (ring) dc ring_def
    gg { (string) dc . } map /ii set
    ii { [ 2 1 roll ww ] decompose_to_w_homogeneous } map /jj set
    jj { dup length 1 eq { pop } { } ifelse } map /jj set
    jj flatten /jj set
    /ans 1 def
    0 1 jj length 1 sub {
       /i set
       jj i get gg reduction-noH 0 get (0). eq {  }
       { jj i get messagen ( does not belong to the ideal ) message
         /ans 0 def
         exit
       } ifelse
    } for
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

[(w_homogeneousQ)
[([ideal variables weight] w_homogeneousQ bool)
 $Example 1: [[(x) (x^2+x) (x^3-y^2)] [(x) (y)] [(x) 1 (y) 1]] w_homogeneousQ$
 $Example 2: [[(x^2+1) (x^3-y^2)] [(x) (y)] [(x) 1 (y) 1]] w_homogeneousQ$
]] putUsages

%% Should move to hol.sm1
/gr_var {
  /arg1 set
  [/in-gr_var /v /ans /i /vec-input] pushVariables
  [
     /v arg1 def
     v isArray {
        /vec-input 1 def
        v { toString } map /v set
     } {
        /vec-input 0 def
        [v to_records pop] /v set
     } ifelse
     /ans v def
     0 1 v length 1 sub {
       /i set
       /ans ans [@@@.Dsymbol v i get] cat append def
     } for
     vec-input not {
       ans from_records /ans set
     } { } ifelse
     /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(gr_var)
[( [v1 ... vn] gr_var [v1 ... vn Dv1 ... Dvn] )
 $ (v1,...,vn) gr_var (v1,...,vn,Dv1,...,Dvn) $
 (cf. wToVW)
]] putUsages

%% Should move to hol.sm1
/reparse {
  /arg1 set
  [/in-reparse /f /ans] pushVariables
  [
     /f arg1 def
     f isArray {
       /ans f { reparse } map def
     }{
       f toString . /ans set
     } ifelse
     /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(reparse)
[(obj reparse obj2)
 (Parse the object in the current ring.)
 (Elements in obj2 belong to the current ring.)
]] putUsages

%% Should move to hol.sm1
/wToVW {
  /arg1 set
  [/in-wToVW /ww /vv /tmp /ans /i] pushVariables
  [
    /tmp arg1 def
    /ww tmp 0 get def
    /vv tmp 1 get def
    /ans [ ] def
    0 1 vv length 1 sub {
      /i set
      ans [ vv i get ww i get (integer) dc] append /ans set
    } for
    /arg1 ans flatten def
  ] pop
  popVariables
  arg1
} def
[(wToVW)
[([ ww vv] wToVW [ v1 w1 ...])
 (cf. gr_var)
 (Example: [ [-1 -2 1 2] [(x) (y) (Dx) (Dy)]] wToVW :: )
]] putUsages

/gr_gb {
  /arg1 set
  [/in-gr_gb /ii /vv /ww /vv_gr /ans /gr_I] pushVariables
  [(CurrentRingp)] pushEnv
  [
    /ii arg1 0 get def
    /vv arg1 1 get def
    /ww arg1 2 get def
    [ii vv ww] gb /ans set
    %% (gr_gb: your gb is) message  ans  message 
    /vv_gr vv gr_var def
    vv_gr isArray { vv_gr from_records /vv_gr set } { } ifelse
    [vv_gr ring_of_polynomials 0] define_ring
    ans 1 get dehomogenize /gr_I set
    gr_I reparse /gr_I set
    /arg1 [ans 0 get gr_I] def
  ] pop
  popEnv
  popVariables
  arg1
} def
[(gr_gb)
[([ii vv ww] gr_gb [ii_gb gr_ii])
 (It computes the Grobner basis ii_gb in D for the weight vector vv.)
 (gr_ii is the initial ideal with respect to ww and is the ideal of)
 (the ring of polynomials with reverse lexicographic order.)
 (The answer is dehomogenized.)
 (cf. gr_var, reparse. Need gb for this function --- load cohom.sm1)
 $Example:   [[(x1*Dx1+2*x2*Dx2+3*x3*Dx3) $
 $             (Dx1^2-Dx2) (-Dx1*Dx2+Dx3)  (Dx2^2-Dx1*Dx3)] $
 $             [ (x1)  (x2)  (x3) ]  ] /ff set $
 $   [ff 0 get ff 1 get [[(x2) -1 (Dx1) 2 (Dx2) 3 (Dx3) 2]]] gr_gb /gg set$
]] putUsages

/firstSlope3 {
  /arg1 set
  [/in-firstSlope3 /ff /gv /gf /wv /wf /vv /vvdd 
   /first-slope /first-weight /first-gb
  ] pushVariables
  [
     /ff arg1 def
     /vv [(x1) (x2) (x3)] def
     /vvdd [(x1) (x2) (x3) (Dx1) (Dx2) (Dx3)] def
     /wf [(Dx1) 1 (Dx2) 1 (Dx3) 1] def %% F-filtration
     /wv [(x2) -1 (Dx2) 1] def         %% V-filtration

     [ff vv [wf]] gb /gf set
     [ff vv [wv]] gb dehomogenize /gv set

     %% determine the first-slope and first-weight here.
     %%  [gf vv [wf]] w_supports_of_I
     %%  [gv vv [wv]] w_supports_of_I
     /firstweight [ (x2) -1 (Dx1) 2 (Dx2) 3 (Dx3) 2] def
     [ff vv [firstweight]] gr_gb 
     /first-gb set
     [ 
      [first-gb 1 get vvdd wf] w_homogeneousQ
      [first-gb 1 get vvdd wv] w_homogeneousQ
      first-gb
     ] /arg1 set
  ] pop
  popVariables
  arg1
} def
%% [ [[1 2 3]] [0]] gkz /ff set ff 0 get firstSlope3 /gg set
%% [ [[1 2 3]] [2]] gkz /ff set ff 0 get firstSlope3 /gg set
%% This program is used to check gr_gb and w_homogeneousQ

/biggest_pq {
  /arg1 set
  [/in-biggest_pq /ex /xmax /ymax /i /ans] pushVariables
  [
     /ex arg1 def
     ex length 0 eq {
       /ans null def
       /LLL.biggest_pq goto
     } {  } ifelse
     /xmax ex 0 get 0 get def
     0 1 ex length 1 sub {
        /i set
        ex i get 0 get xmax ge {
          /xmax ex i get 0 get def
          /ymax ex i get 1 get def
        }{ } ifelse
     } for 
     0 1 ex length 1 sub {
        /i set
        ex i get 0 get xmax eq {
          ex i get 1 get ymax gt {
            /ymax ex i get 1 get def
          } { } ifelse
        }{ } ifelse
     } for 
     /ans [xmax ymax] def
     /LLL.biggest_pq
     /arg1 ans def
  ] pop
  popVariables
  arg1
}def
[(biggest_pq)
[([[i1 j1] [i2 j2] ...] biggest_pq [ik jk])
 (It returns the biggest [i j] with the lexicographic order x > y)
 (Example:  [ [1 2] [1 3] [2 4] [2 -1]] biggest_pq :: )
]] putUsages

/remove_x* {
  /arg1 set
  [/in-remove_x* /ans /i /ex /x] pushVariables
  [
    /ex arg1 0 get def
    /x arg1 1 get def
    /ans [  ] def
    0 1 ex length 1 sub {
      /i set
      ex i get 0 get x eq { 
      }{
        /ans ans ex i get append def
      } ifelse
    } for
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(remove_x*)
[([[[i1 j1] [i2 j2] ...] x] remove_x* [[i1 j1] [i2 j2] ...])
 (It removes [x *] elements from [[i1 j1] ...])
 (Example: [ [ [1 2] [1 3] [2 4] [2 -1]] 2 ] remove_x* :: )
]] putUsages

% f > g ?
/greater_u {
  /arg2 set  /arg1 set
  [/in-greater_u /f /g /tmp /ans] pushVariables
  [
    /f arg1 def /g arg2 def
    f g sub /tmp set
    /ans 0 def
    tmp isInteger {
      tmp 0 gt {
         /ans 1 def
      }{  } ifelse
    }{
      tmp isRational { tmp (numerator) dc /tmp set } { } ifelse
      tmp (0).. gt {
         /ans 1 def
      } {  } ifelse
    } ifelse
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

%% to turn around the a bug of univ-num (universalNumber) dc bug.
/toUniv {
  /arg1 set
  [/in-toUniv /p] pushVariables
  [
    /p arg1 def
     p isInteger {
       /p p (universalNumber) dc def
     }{ } ifelse
    /arg1 p def
  ] pop
  popVariables
  arg1
} def
/smallSlope {
  /arg1 set
  [/in-smallSlope /ex /p /q /tmp /r /s /slope
   /upperBoundOfSlope
  ]  pushVariables
  [
     /ex arg1 0 get def
     /upperBoundOfSlope arg1 1 get def
     (0).. upperBoundOfSlope greater_u {
         (SmallSlope: the upperBoundOfSlope has a negative value.)
         error
     } {  } ifelse
     /slope (0).. def
     /tmp ex biggest_pq def
     /p tmp 0 get def /q tmp 1 get def
     [ex p] remove_x* /ex set
     {
        ex length 0 eq { exit } { } ifelse
        /tmp ex biggest_pq def
        /r tmp 0 get def  %% tmp = (r,s)
        /s tmp 1 get def  %% tmp = (r,s)
        [ex r] remove_x* /ex set
        s q greater_u {
          %% return  (s-q)/(p-r) : positiive
          s q sub toUniv 
          p r sub toUniv div /slope set
          [(cancel) slope] mpzext /slope set
          upperBoundOfSlope slope greater_u {
            exit
          } {
            /p r def
            /q s def
            /slope (0).. def % throw away this slope
          } ifelse
        } {  } ifelse     
     } loop
     /arg1 slope def
  ] pop
  popVariables
  arg1
} def
[(smallSlope)
[([ [[i1 j1] [i2 j2] ...]  upperBound] smallSlope b/a)
 (The absolute value of the smallSlope must be smaller than upperBound.)
 (Example: [ [[1 2] [1 6] [2 4] [2 -1]] slope.infinity] smallSlope :: )
 (Example: [ [[0 7] [1 2] [1 6] [2 4] [2 -1]] (2)..] smallSlope :: )
 (Example: [ [[1 2] [1 3] [2 4] [2 -1]] slope.infinity]smallSlope :: )
 (Example: [ [[1 2] [1 -1]] slope.infinity] smallSlope :: )
 $Example: [ [[1 2 3]] [0]] gkz /ff set$
 $          [ ff 0 get ff 1 get [ [(Dx1) 1 (Dx2) 1 (Dx3) 1] [(Dx3) 1 (x3) -1]]$
 $          ] w_supports_of_I /gg set$
 $        gg 0 get { /pp set [pp slope.infinity] smallSlope } map /hh set $
]] putUsages


/maxSlope {
  /arg1 set
  [/in-maxSlope /ss /ans /i] pushVariables
  [
    /ss arg1 def
    /ans (0).. def
    0 1 ss length 1 sub {
      /i set
      ss i get ans greater_u {
        /ans ss i get def
      } {  } ifelse
    } for
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def

/slope {
  /arg1 set
  [/in-slope /ff /gv /gf /wv /wf /wll /worderf
   /vv /vvdd  /f /v /ll /f-filt
   /w_supp
   /virtualSlope /a /b  /ans /tmp /sslopes
   /pp  /maxSmallSlope
   /first-slope /first-weight /first-gb /first-init
  ] pushVariables
  [
     /ff arg1 0 get def
     /vv arg1 1 get def
     /f  arg1 2 get def
     /v  arg1 3 get def
     vv isArray not { [vv to_records pop] /vv set } { } ifelse
     /f-filt f def
%%   Example:
%%     /ff [ (2 y Dx + 3 x^2 Dy) (3 y^3 Dy - 2 x^4 Dx - 6 x^3 y Dy + 6)] def
%%     /f [ 0 0 1 1] def   %% F-filtration
%%     /v [ -1 0 1 0] def  %% V-filtration
%%     /vv [(x) (y)] def
     %% -3: x=0,   -2 : y =0

     /maxSmallSlope slope.infinity def
     /vvdd vv gr_var def
     vvdd length f length eq {  }
     { (The number of variables <<vvdd>> and the size of weight vector <<f>>do not match.) 
        error } ifelse
     vvdd length v length eq {  }
     { (The number of variables <<vvdd>> and the size of weight vector <<v>>do not match.) 
        error } ifelse
     /ans [ ] def
     /wv [v vvdd] wToVW def

     /worderf [f vvdd] wToVW def

     /wf [f vvdd] wToVW def
     slope.verbose { (Computing gb with ) messagen wf message ( and ) messagen
                      wv message } { } ifelse
     [ff vv [wf wv]] gr_gb /first-gb set
     /firstweight wf def
    {
     /wf [f vvdd] wToVW def

     first-gb 0 get dehomogenize /gf set
     [gf vv [worderf wv]] w_supports_of_I_without_gb_computation
     /w_supp set
      slope.verbose { (w_supp are ) message w_supp 0 get message } { } ifelse
      slope.verbose { (gb is ) message w_supp 1 get message } { } ifelse
      slope.verbose { (weight is ) messagen firstweight message } { } ifelse
      w_supp 0 get { /pp set [pp maxSmallSlope] smallSlope } map /sslopes set
      slope.verbose { (smallSlopes are ) message sslopes message } { } ifelse
      sslopes maxSlope /first-slope set
      first-slope (0).. greater_u {
       (small slope is ) messagen first-slope message
     } {  
       (All the smallSlopes are zero. Exiting...) message
        exit
     } ifelse
     /a first-slope (denominator) dc def
     /b first-slope (numerator) dc def
     %% a v mul  b f mul add /ll set
     a v mul  b f-filt mul add /ll set
     /firstweight [ll vvdd] wToVW def
     (Computing the GB with the weight vector ) messagen firstweight message
     (and ) messagen wv message
     [ff vv [firstweight wv]] gr_gb  % use two weight vectors.
     /first-gb set
     %% (GB is) messagen first-gb message
     first-gb 1 get /first-init set
     slope.geometric {
       (To get the geometric slope, we need to compute the radical.) message
       [ first-init vvdd] radical /first-init set
       [first-init vvdd] pgb 0 get /first-init set
       (Radical is ) messagen first-init message
     } { } ifelse
     [first-init vvdd worderf] w_homogeneousQ
     [first-init vvdd wv] w_homogeneousQ
     and {
       (It is bi-homogeneous!  It is not a slope.) message
       /maxSmallSlope first-slope def %% I think it is necessary.
     } {
       slope.geometric {
          (It is a geometric slope.) message
       }{
         (It is an algebraic slope.) message
       } ifelse
       /maxSmallSlope first-slope def
       /ans ans [first-slope ll] append def
     } ifelse
     (-----------------------------------------------) message
     /f ll def
    } loop
    /arg1 ans def
  ] pop
  popVariables
  arg1
} def
[(slope)
[( [ii vv F-filtration V-filtration] slope [ [-slope1 weight] ...])
 ( ii : ideal, vv : variables, F-filtration : F-filtration by vector)
 ( V-filtration : V-filtration by vector)
 (It computes the algebraic or geometric slopes of ii along the hyperplane)
 (specified by the V-filtration.)
 (When slope.geometric is one, it outputs the geometric slopes.)
 (As to the algorithm, see A.Assi, F.J.Castro-Jimenez and J.M.Granger)
 ( How to calculate the slopes of a D-module, Compositio Math, 104, 1-17, 1996)
 (Note that the signs of the slopes s' are negative, but the absolute values -s')
 (of the slopes are returned.)
 $In other words, when pF+qV is the gap, -s'=q/p is returned.$
 $Note that s=1-1/s' is called the slope in recent literatures. Solutions belongs to O(s).$
 $The number s satisfies 1<= s.$  
 $We have r=s-1=-1/s', and kappa=1/r=-s',$
 $which is used 1/Gamma(1+m*r) factor and exp(-tau^kappa)$
 $in the Borel and Laplace transformations respectively.$
 $Example 1: [ [(x^4 Dx + 3)] (x) [0 1] [-1 1]] slope :: $
 $          The solution is exp(x^(-3)). $
 $Example 2: [ [(x^3 Dx^2 + (x + x^2) Dx + 1)] [(x)] $
 $             [0 1] [-1 1]] slope :: $
 $Example 3: [ [(x^6 Dx^3 + x^3 Dx^2 + (x + x^2) Dx + 1)] [(x)] $
 $             [0 1] [-1 1]] slope :: $
 $Example 4:$ 
 $   /ff [ (2 y Dx + 3 x^2 Dy) (3 y^3 Dy - 2 x^4 Dx - 6 x^3 y Dy + 6)] def$
 $   [ ff (x,y) [ 0 0 1 1]  [ 0 -1 0 1] ] slope :: $
 $   Answer should be  2  ==> -2 $
 $Example 5:$ 
 $   /ff [ [[1 2 3]] [-3]] gkz def $
 $   [ ff 0 get ff 1 get [ 0 0 0 1 1 1]  [ 0 0 -1 0 0 1] ] slope :: $
]] putUsages

/bihomogeneousGrQ {
  /arg1 set
  [/in-checkBihomogeneous  /ff /vv /firstweight /worderf /wv
   /first-gb /ans /vvdd
  ] pushVariables
  [
     arg1 0 get /ff set
     arg1 1 get /vv set
     arg1 2 get /firstweight set
     arg1 3 get 0 get /worderf set
     arg1 3 get 1 get /wv set

     vv isArray not { [vv to_records pop] /vv set} { } ifelse
     vv gr_var /vvdd set
     %%(Computing the GB with the weight vector ) messagen firstweight message
     [ff vv [firstweight]] gr_gb 
     /first-gb set
     %% (GB is) messagen first-gb message
     [first-gb 1 get vvdd worderf] w_homogeneousQ
     [first-gb 1 get vvdd wv] w_homogeneousQ
     and {
       (It is bi-homogeneous!) message  /ans 1 def
     } {
       (It is not bi-homogenous w.r.t ) messagen
        [worderf wv] message
        /ans 0 def
     } ifelse
     /arg1 [ans first-gb firstweight] def
  ] pop
  popVariables
  arg1
} def
[(bihomogeneousGrQ)
[([ ii vv w [vf wv]] bihomogeneousGrQ [ans gb])
 $It checks if in_w(ii) is bihomogeneous w.r.t. vf and wv$
 $Example 1: [ [[1 2 3]] [0]] gkz /ff set $
 $         [ff 0 get ff 1 get [(x3) -2 (Dx1) 1 (Dx2) 1 (Dx3) 3] $
 $              [ [(Dx1) 1 (Dx2) 1 (Dx3) 1] [(x3) -1 (x3) 1]]] $
 $              bihomogeneousGrQ /gg set $
 $ bi-homogeneous $
 $Example 2: [ [[1 2 3]] [0]] gkz /ff set $
 $         [ff 0 get ff 1 get [(x3) -1 (Dx1) 2 (Dx2) 2 (Dx3) 3] $
 $              [ [(Dx1) 1 (Dx2) 1 (Dx3) 1] [(x3) -1 (x3) 1]]] $
 $              bihomogeneousGrQ /gg set $
 $ not bi-homogeneous $
 $Example 3: [ [[1 3]] [0]] gkz /ff set $
 $         [ff 0 get ff 1 get [(x2) -2 (Dx1) 1 (Dx2) 3] $
 $              [ [(Dx1) 1 (Dx2) 1] [(x2) -1 (x2) 1]]] $
 $              bihomogeneousGrQ /gg set $
 $ not bi-homogeneous $
]] putUsages

%% Radical via primary ideal decomposition.
/radical {
  /arg1 set
  [/in-radical /ii /jj /pp0 /n /i /vv /ans] pushVariables
  [
     /ii arg1 def
     ii 1 get /vv set
     ii primadec /jj set
     /n  jj length def
     jj { 1 get } map /pp0 set
     vv isArray {
       /vv vv from_records def
     } { } ifelse
     (Primary components are ) messagen pp0 message
     /ans pp0 0 get def
     pp0 rest /pp0 set
     {
        pp0 length 0 eq { exit } {  } ifelse
        %% [ans pp0 0 get vv] message
        [ans pp0 0 get vv] gr_intersection /ans set
        %%[ans pp0 0 get vv] gr_intersection /ans set
        pp0 rest /pp0 set
     } loop
     ans /arg1 set
  ] pop
  popVariables
  arg1
} def
[(radical)
[([ii vv] radical jj)
 (Computing the radical of ii via primadec.)
 (Example 1: [ [(x^2-1) (x^4-1)] (x)] radical ::)
 (Example 2: [ [(x^2 y) (y^4)  (x y)] (x,y)] radical ::)
]] putUsages

/gr_intersection {
  /arg1 set
  [/in-gr_intersection /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_polynomials
      [[(_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
[(gr_intersection)
[(Ideal intersections in the ring of polynomials.)
 $Example 1: [[(y) (Dx)] [(x) (Dy)] (x,y,Dx,Dy)] gr_intersection ::$
]] putUsages


/tests {

     /ff [ [[1 2 3] ] [0]] gkz 0 get def
     /vv [(x1) (x2) (x3)] def
     /f [ 0 0 0 1 1 1] def   %% F-filtration
     /v [ 0 0 -1  0 0 1] def  %% V-filtration

     /ff [ [[1 2 4] ] [0]] gkz 0 get def
     /vv [(x1) (x2) (x3)] def
     /f [ 0 0 0 1 1 1] def   %% F-filtration
     /v [ 0 0 -1  0 0 1] def  %% V-filtration

     %% [1 2 3]
     /ff  [    $2*(x1-1)*Dx1+4*(x2-2)*Dx2+6*x3*Dx3-1$ , $Dx1^2-Dx2$ , $-Dx1*Dx2+Dx3$ , $Dx2^2-Dx1*Dx3$ ] def
     /vv [(x1) (x2) (x3)] def
     /f [ 0 0 0 1 1 1] def   %% F-filtration
     /v [ 0 0 -1  0 0 1] def  %% V-filtration

     %% [1 2 4]
     /ff  [    $2*(x1-1)*Dx1+4*(x2-2)*Dx2+8*x3*Dx3-1$ , $Dx1^2-Dx2$ , $Dx2^2-Dx3$ ] def
     /vv [(x1) (x2) (x3)] def
     /f [ 0 0 0 1 1 1] def   %% F-filtration
     /v [ 0 0 -1  0 0 1] def  %% V-filtration

     /ff [ (2 y Dx + 3 x^2 Dy) (3 y^3 Dy - 2 x^4 Dx - 6 x^3 y Dy + 6)] def
     /f [ 0 0 1 1] def   %% F-filtration
     /v [ 0 -1 0 1] def  %% V-filtration
     /vv [(x) (y)] def
     %% -3: x=0,   -2 : y =0

     /ff [ [[1 3]] [0]] gkz 0 get def
     /f [ 0 0 1 1] def   %% F-filtration
     /v [ 0 -1 0 1] def  %% V-filtration
     /vv [(x1) (x2)] def


} def

/slope.sm1.loaded 1 def