%%% oaku/kan/bfunction.sm1, 1998, 11/5

%%% global variables for bfunction
%%% bfunction.*
/bfunction.version (2.981105) def
bfunction.version [(Version)] system_variable gt 
{ (This package requires the latest version of kan/sm1) message
  (Please get it from http://www.math.kobe-u.ac.jp/KAN) message
  error
} { } ifelse
/bfunction.v [(x) (y) (z)] def   %% default variables of the input polynomial 
/bfunction.s (s) def             %% default variable of the output b-function
/bfunction.vh (v_) def           %% default variable for V-homogenization
/bfunction.t (t_) def            %% default variable for t in delta(t-f)
/bfunction.a  [] def             %% parameters are not available yet
/bfunction.verbose 0 def         %% no messages if 0
/bfunction.strategy 0 def        %% V-homogenization + h-homogenization if 0
                                 %% V-homogenization if 1 (not available yet)
                                 %% h-homogenization if 2 (not available yet)
/bfunction.result 0 def
         
(bfunction.sm1,  11/05,1998  (C) T. Oaku.  bfunction ) message-quiet

[(bfunction)
 [( a bfunction b)
  (array a; poly b;)
  (a :  [f] ;  string f ;)
  (a :  [f] ;  polynomial f ;)
  (a :  [f v] ; string f,v; )
  (a :  [f v] ; polynomial f, string v; )
  (b is the b-function (=Bernstein-Sato polynomial) of a polynomial f)
  (in variables v.)
  (If v is not specified, the variables are assumed to be (x,y,z). )
  (b will be a polynomial in s.  This variable can be changed by typing in)
  (  (variable) /bfunction.s set )
  (For the algorithm, see Duke Math. J. 87 (1997),115-132,)
  (   J. Pure and Applied Algebra 117&118(1997), 495--518.)
  $Example  [(x^3-y^2) (x,y)] bfunction :: $
  (  )
  (See also bfct which implements a new algorithm to compute b-function and is faster. Aug 2002.)
 ]
] putUsages

/bfunction {
  /arg1 set
  [/aa /typev /setarg /f /s /v /bf /bfs /vt ] pushVariables 
  [(CurrentRingp) (KanGBmessage)] pushEnv  %% push current global environment.
  [

    /aa arg1 def
    aa isArray { } { (array bfunction) message error } ifelse
    /setarg 0 def
    aa { tag } map /typev set
    typev [ StringP ] eq
    {  /f aa 0 get def
       /v bfunction.v def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    typev [ PolyP ] eq
    {  /f aa 0 get (string) data_conversion def
       /v bfunction.v def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    typev [StringP StringP] eq
    {  /f aa 0 get def
       /v [ aa 1 get to_records pop ] def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    typev [PolyP StringP] eq
    {  /f aa 0 get (string) data_conversion def
       /v [ aa 1 get to_records pop ] def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    typev [StringP ArrayP] eq
    {  /f aa 0 get def 
       /v aa 1 get def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    typev [PolyP ArrayP] eq
    {  /f aa 0 get (string) data_conversion def 
       /v aa 1 get def
       /s bfunction.s def
       /setarg 1 def
    } { } ifelse
    setarg { } { (Argument mismatch) message error } ifelse

    [(KanGBmessage) bfunction.verbose] system_variable

    v bfunction.t append /vt set

    [f v fw_delta bfunction.t vt] indicial 0 get /bf set
    [bfunction.s ring_of_polynomials 0] define_ring 
    bf . /bf set  
    bfunction.s . /bfs set
    bf [[bfs (-1). bfs sub]] replace /bf set 
    /bfunction.result bf def
    /arg1 bf def  
  ] pop
  popEnv
  popVariables
  arg1
} def

%%  Computing the indicial polynomial (the b-function) of a D-module
/indicial {
  /arg1 set  %% [equations, the variable to be restricted to 0, all variables]
  [/eqs /t /vars /allvars /newvars /x_vars /ans1 /ans2 ] pushVariables 
  [(CurrentRingp)] pushEnv  
  [
    arg1 0 get /eqs set
    arg1 1 get /t set
    arg1 2 get /vars set
    vars bfunction.s append /allvars set 
    [bfunction.t] allvars complement /newvars set 
    [bfunction.t] vars complement /x_vars set 
    [eqs t vars] indicial1 /ans1 set
    [ans1 x_vars newvars] eliminate_Dx /ans2 set
    [ans2 x_vars newvars] eliminate_x /arg1 set
  ] pop
  popEnv
  popVariables
  arg1
} def

%% (-1,0;1,0)-Groebner basis
%% [equations (t) vars] indical1 ---> psi(BFequations) (as a list of strings)
/indicial1 {
 /arg1 set  
 [/bft /bfs /bfh /bf1 /ff /ans /n /i /BFallvarlist /BFDvarlist
  /BFs_weight /BFvarlist ] pushVariables
 [(CurrentRingp)] pushEnv   
 [
  /ff arg1 0 get def
  /bft arg1 1 get def
  /BFvarlist arg1 2 get def
  /BFallvarlist 
    [ bfunction.vh bfunction.s] BFvarlist concat bfunction.a concat 
  def 
  BFvarlist length /n set
  BFvarlist {xtoDx} map /BFDvarlist set
  /BFs_weight 
    [ [ bfunction.vh 1 ]
      [ 0 1 n 1 sub 
          { /i set BFDvarlist i get 1 } 
        for 
        0 1 n 1 sub  
          { /i set BFvarlist i get 1 }
        for ]
    ] def  

  [ BFallvarlist listtostring ring_of_differential_operators
    BFs_weight weight_vector 
  0] define_ring  

  /bfh  (h). def
  /bfs  bfunction.vh . def
  /bf1  (1). def
  ff { bft fw_homogenize . } map /ff set
  ff {[[bfh bf1]] replace} map {homogenize} map /ff set 
  [ff] groebner 0 get {[[bfh bf1]] replace} map /ff set 
  ff reducedBase /ans set 
  ans {bft fw_psi} map /ans set
  ans {(string) data_conversion} map /arg1 set
  ] pop
  popEnv
 popVariables
 arg1
} def

%% eliminates Dx in the ring of differential operators
/eliminate_Dx {
 /arg1 set  %% [operators x variables]
 [/bfh /bf1 /ff /ans /nx /ny /x_varlist /Dx_weight /BFvarlist 
   /allvarlist /Dx_varlist /y_varlist /Dy_varlist /allvarlist /i 
 ] pushVariables 
 [(CurrentRingp)] pushEnv  
 [
  /ff arg1 0 get def
  /x_varlist arg1 1 get def
  /BFvarlist arg1 2 get def
  x_varlist length /nx set
  BFvarlist bfunction.a concat /allvarlist set

  x_varlist {xtoDx} map /Dx_varlist set
  x_varlist BFvarlist complement /y_varlist set
  y_varlist length /ny set
  y_varlist {xtoDx} map /Dy_varlist set

  /Dx_weight 
    [ [ 0 1 nx 1 sub 
          { /i set Dx_varlist i get 1 } 
        for ]
      [ 0 1 nx 1 sub 
          { /i set x_varlist i get 1 } 
        for 
        0 1 ny 1 sub 
          { /i set y_varlist i get 1 } 
        for 
        0 1 ny 1 sub 
          { /i set Dy_varlist i get 1 } 
        for 
      ]
    ] def
  
  [ allvarlist listtostring  ring_of_differential_operators
    Dx_weight weight_vector 
  0] define_ring  

  /bfh (h). def
  /bf1 (1). def
  ff {.} map /ff set
  ff {[[bfh bf1]] replace} map {homogenize} map /ff set
  bfunction.verbose 1 eq 
    {(Eliminating the derivations w.r.t. ) messagen x_varlist ::}
    { } 
  ifelse
  [ff] groebner 0 get {[[bfh bf1]] replace} map /ff set  
  ff reducedBase /ans set
  ans Dx_varlist eliminatev /ans set
  ans {(string) data_conversion} map /arg1 set
 ] pop
 popEnv
 popVariables
 arg1
} def

%% eliminates x in the ring of polynomials
/eliminate_x {
 /arg1 set  %% [operators x variables]
 [/bfh /bfs /bf1 /ff /ans /nx /ny /x_varlist  /BFvarlist 
   /allvarlist /y_varlist /i 
 ] pushVariables 
 [(CurrentRingp)] pushEnv  
 [
  /ff arg1 0 get def
  /x_varlist arg1 1 get def
  /BFvarlist arg1 2 get def
  x_varlist length /nx set
  BFvarlist bfunction.a concat /allvarlist set

  x_varlist BFvarlist complement /y_varlist set
  y_varlist length /ny set

  /x_weight 
    [ [ 0 1 nx 1 sub 
          { /i set x_varlist i get 1 } 
        for ]
      [ 0 1 ny 1 sub 
          { /i set y_varlist i get 1 } 
        for 
      ]
    ] def
  
  [ allvarlist listtostring  ring_of_polynomials x_weight weight_vector 
  0] define_ring  

  /bfh (h). def
  /bf1 (1). def
  ff {.} map /ff set
  ff {[[bfh bf1]] replace} map {homogenize} map /ff set
  bfunction.verbose 1 eq 
    {(Eliminating the variables ) messagen x_varlist ::}
    { } 
  ifelse
  [ff] groebner 0 get {[[bfh bf1]] replace} map /ff set  
  ff reducedBase /ans set
  ans x_varlist eliminatev /ans set
  ans {(string) data_conversion} map /arg1 set
 ] pop
 popEnv
 popVariables
 arg1
} def
%%%%%%%%%%%%%%%%%%%%%%% libraries  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
%% FW-principal part of an operator (FW-homogeneous)
%%  Op (poly) fw_symbol --->  FW-symbol(Op)  (poly)
/fw_symbol {
  [[(h). (1).]] replace bfunction.vh . coefficients 1 get 0 get
} def

%% FW-homogenization 
%% Op (string) (t) fw_homogenize ---> h(Op) (string)
/fw_homogenize {
  /arg2 set  %% bft (string)
  /arg1 set  %% an operator (string)
  [ /bft /bfDt /bfht /bfhDt /Op /degs /m /mn ] pushVariables
  [
    /Op arg1 expand def
    /bft arg2 def                         
    bft xtoDx /bfDt set                   
    bfunction.vh (^(-1)*) bft 3 cat_n /bfht set    
    bfunction.vh (*) bfDt 3 cat_n /bfhDt set       
    Op [[bft expand bfht expand][bfDt expand bfhDt expand]] replace 
      /Op set
    Op bfunction.vh expand coefficients 0 get 
      {(integer) data_conversion} map /degs set 
    degs << degs length 1 sub >> get /m set
    0 m sub /mn set  
    << bfunction.vh expand mn powerZ >> Op mul /Op set 
    Op (string) data_conversion /arg1 set
  ] pop
  popVariables
  arg1
} def

%% setup the ring of differential operators with the variables varlist 
%% and parameters bfunction.a
%% varlist setupBFring
/setupDring {
  /arg1 set  
  [ /varlist /bft /allvarlist /n /dvarlist /D_weight /i
  ] pushVariables
  [  
    arg1 /varlist set
    /allvarlist 
      varlist bfunction.a join 
    def                            
    varlist length /n set
    varlist {xtoDx} map /dvarlist set
    /D_weight 
    [ [ 0 1 n 1 sub 
          { /i set dvarlist i get 1 } 
        for ]
      [ 
        0 1 n 1 sub  
          { /i set varlist i get 1 }
        for ]
    ] def  

    [ allvarlist listtostring ring_of_differential_operators
      D_weight weight_vector 
    0] define_ring 
    
  ] pop
  popVariables
} def
 
%% psi(P)(s) 
%% Op (poly) (t) (string) fw_psi ---> psi(P) (poly)
%% Op should be FW-homogeneous. 
/fw_psi {
 /arg2 set  %% bft (string)
 /arg1 set  %% Op  (polynomial)
 [/bft /bfDt /P /tt /dtt /k /Q /i /m /kk /PPt /PPC /kk /Ss] pushVariables
 [
  arg2 /bft set 
  arg1 fw_symbol /P set 
  /bfDt bft xtoDx def
  /tt bft expand def  /dtt bfDt expand def               
  P bft fw_order /k set    
    << 1 1 k >> 
    {pop tt P mul /P set }
    for
    << -1 -1 k >>
    {pop dtt P mul /P set }
    for 
  (0) expand /Q set
  P dtt coefficients 0 get length /m set
  0 1 << m 1 sub >>  
  {
    /i set 
    P dtt coefficients 0 get i get /kk set 
    kk (integer) data_conversion /kk set
    P dtt coefficients 1 get i get /PPt set
    PPt tt coefficients 1 get 0 get /PPC set
    bfunction.s expand /Ss set
    0 1 << kk 1 sub >> { 
      pop
      PPC Ss mul /PPC set
      Ss (1) expand sub /Ss set
    } for
    Q PPC add /Q set
  } for
  Q  /arg1 set
 ] pop
 popVariables
 arg1
} def

%% get the FW-order
%% Op (poly) (t) fw_order ---> FW-ord(Op) (integer)
%% Op should be FW-homogenized. 
/fw_order {
 /arg2 set  %% bft (string)
 /arg1 set  %% Op (poly)
 [/Op /bft /fws /m /fwsDt /k /tt /dtt] pushVariables
 [
  arg1 /Op set
  arg2 /bft set
  Op fw_symbol /fws set
  /tt bft expand def
  /dtt bft xtoDx  expand def
  fws [[bfunction.s expand  (1).]] replace /fws set
  fws dtt coefficients 0 get 0 get /m set
  fws dtt coefficients 1 get 0 get /fwsDt set
  fwsDt tt coefficients 0 get 0 get /k set
  m k sub (integer) data_conversion /arg1 set
 ] pop
 popVariables
 arg1
} def

/remove0 {
  /arg1 set
  arg1 (0). eq
  { } {arg1} ifelse
} def

%% functions for list operations etc. 

/notidentical {
  /arg2 set
  /arg1 set
  arg1 arg2 eq
  { } {arg1} ifelse
} def

%% [(x1) (x2) (x3)] ---> (x1,x2,x3)
/listtostring {
  /arg1 set
  [/n /j /ary /str] pushVariables 
  [
    /ary arg1 def
    /n ary length def
    arg1 0 get /str set
    n 1 gt  
      { str (,) 2 cat_n /str set }{ }
    ifelse
    1 1 n 1 sub {
      /j set
      j n 1 sub eq 
        {str << ary j get >>  2 cat_n /str set}      
        {str << ary j get >>  (,) 3 cat_n /str set}
      ifelse
    } for
    /arg1 str def
  ] pop
  popVariables
  arg1
} def

%% (x1) --> (Dx1)
/xtoDx {
  /arg1 set
  @@@.Dsymbol arg1 2 cat_n
} def

%% concatenate two lists
/concat {
  /arg2 set
  /arg1 set
  [/n /j /lst1 /lst2 ] pushVariables
  [
    /lst1 arg1 def
    /lst2 arg2 def
    /n lst2 length def
    0 1 n 1 sub {
      /j set
      lst1 lst2 j get append /lst1 set
    } for
    /arg1 lst1 def
  ] pop
  popVariables
  arg1
} def

%% var (poly) m (integer) ---> var^m (poly)
/powerZ {
  /arg2 set %% m
  /arg1 set %% Var
  [ /m /var /varstr /pow /nvar] pushVariables
  [
    arg1 /var set
    arg2 /m set
    var (string) data_conversion /varstr set
    m -1 gt
      { var m npower /pow set}
      { varstr (^(-1)) 2 cat_n expand /nvar set
        nvar << 0 m sub >> npower /pow set 
       }
    ifelse
    pow /arg1 set
  ] pop
  popVariables
  arg1
} def


%% (f) varlist fw_delta ---> [t - f, Dx + f_xDt, ...] 
/fw_delta {
  /arg2 set %% [(x) (y) ...]
  /arg1 set %% (f)
  [ /fstr /f /bft /n /j /varlist /dxvarlist /allvarlist /xi /fxi /dxi /dt
    /delta /BFdt /BFDtx_weight ] pushVariables
  [
    arg1 /fstr set
    arg2 /varlist set
    [bfunction.t] varlist join bfunction.a join /allvarlist set
    bfunction.t xtoDx /BFdt set 
    varlist {xtoDx} map /dxvarlist set
    varlist length /n set
    /BFDtx_weight [ [ BFdt 1
      0 1 n 1 sub {/j set varlist j get 1} for ]
      [ bfunction.t 1
      0 1 n 1 sub {/j set dxvarlist j get 1} for ]
    ] def   

    [ allvarlist listtostring ring_of_differential_operators
     BFDtx_weight weight_vector  0 ] define_ring

    fstr expand /f set 
    bfunction.t expand /bft set 
    BFdt expand /dt set 
    /delta [
      bft f sub
      0 1 n 1 sub {
        /i set
        varlist i get xtoDx expand /dxi set
	<< dxi f mul >> << f dxi mul >> sub [[(h). (1).]] replace /fxi set
        dxi << fxi dt mul >> add 
      } for
    ] def
    delta {(string) data_conversion} map /arg1 set
  ] pop
  popVariables
  arg1
} def