Return to slope.sm1 CVS log

Up to [local] / OpenXM / src / kan96xx / Doc

The function slope computes the slopes of a given D-module.
It uses the ACG algorithm (Assi,Castro,Granger, How to calculate
the slopes of a D-module, Compositio Math (1996), 107--123) to
get the algebraic slopes.  In order to get the geometric slopes,
the ACG algorithm should be modified and repeated use of radical
computation is necessary.  Radical computation is done by ox_asir
module via OpenXM connection.
This module is a joint work with Castro-Jimenez.

% $OpenXM: OpenXM/src/kan96xx/Doc/slope.sm1,v 1.1 2000/11/01 01:57:55 takayama Exp $
(cohom.sm1) run (oxasir.sm1) run
$slope.sm1, computing the slopes of a D-ideal:  June 15, 2000$ message
$                            (C) F.Castro-Jimenez, N.Takayama$ 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 are negative, but the absolute values)
 (of the slopes are returned.)
 $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