===================================================================
RCS file: /home/cvs/OpenXM/src/kan96xx/Doc/gfan.sm1,v
retrieving revision 1.9
retrieving revision 1.12
diff -u -p -r1.9 -r1.12
--- OpenXM/src/kan96xx/Doc/gfan.sm1	2005/06/30 08:39:39	1.9
+++ OpenXM/src/kan96xx/Doc/gfan.sm1	2005/07/07 07:53:37	1.12
@@ -1,6 +1,6 @@
-%  $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.8 2004/10/13 23:36:52 takayama Exp $
+%  $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.11 2005/07/07 06:07:46 takayama Exp $
 % cp cone.sm1 $OpenXM_HOME/src/kan96xx/Doc/gfan.sm1
-% $Id: gfan.sm1,v 1.9 2005/06/30 08:39:39 takayama Exp $
+% $Id: gfan.sm1,v 1.12 2005/07/07 07:53:37 takayama Exp $
 % iso-2022-jp
 %%Ref:  @s/2004/08/21-note.pdf
 
@@ -17,6 +17,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 /cone.sample {
   cone.load.cohom
+  /cone.ckmFlip  1 def
 % write a comment about the problem.  "nl" means new line.  
 /cone.comment [
   (Toric ideal for 1-simplex x 2-simplex, in k[x]) nl
@@ -76,6 +77,7 @@ def   
   ]
 def
 
+/cone.DhH  0 def
 % Set a function to compute Grobner basis.
 %  cone.gb_Dh   : For computing in Homogenized Weyl algebra h[1,1](D).
 %  cone.gb_DhH  : For computing in doubly homogenized Weyl algebra. 
@@ -128,6 +130,7 @@ printGrobnerFan 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 /cone.sample2 {
   cone.load.cohom
+  /cone.ckmFlip  1 def
 % write a comment about the problem.  "nl" means new line.  
 /cone.comment [
   (BS for y and y-(x-1)^2, t1, t2 space, in doubly homogenized Weyl algebra.) nl
@@ -194,6 +197,7 @@ def
   cone.input { . homogenize toString } map /cone.input set
   dh.end
 
+/cone.DhH  1 def
 % Set a function to compute Grobner basis.
 %  cone.gb_Dh   : For computing in Homogenized Weyl algebra h[1,1](D).
 %  cone.gb_DhH  : For computing in doubly homogenized Weyl algebra. 
@@ -398,8 +402,36 @@ cone.comment message
 %         cone.d (pointed cones lies in this space. cf. cone.Lp)
 % These are set during getting the cone.startingCone
 
+%<
+% global
+%cone.ckmFlip. Collar-Kalkbrener-Mall の flip アルゴリズムを使わない 0. 使う 1.
+%  Default は 0.
+%> 
+/cone.ckmFlip 0 def
 
 %<
+% global
+% cone.DhH  dx x = x dx + h H なら 1. dx x = x dx + h^2 なら 0. Default 0.
+%>
+/cone.DhH  0 def
+
+%<
+% Global
+% gbCheck をするか? しないと結果はあやふや. しかしメモリ exhaust は防げる.
+% 使うときは /cone.epsilon,  /cone.epsilon.limit を十分小さくしておく.
+%>
+/cone.do_gbCheck 1 def
+
+% Default の cone.gb の定義. 各プログラムで再度定義してもよい.
+/cone.gb {
+  cone.DhH { 
+     cone.gb_DhH
+  } {
+     cone.gb_Dh
+  } ifelse
+} def
+
+%<
 % Usage:  wv g coneEq1
 % in(f) が monomial 専用.  in_w(f) = LT(f) となる weight w の満たす
 % 不等式制約を求める.
@@ -2023,15 +2055,20 @@ def
 %>
 /cone.gb_Dh {
   /arg2 set /arg1 set
-  [/ff /ww /gg] pushVariables
+  [/ff /ww /gg /gbopt] pushVariables
   [
     /ff arg1 def
     /ww arg2 def
     [(AutoReduce) 1] system_variable
     [cone.vv ring_of_differential_operators
      [ww] weight_vector 0] define_ring
-    [ff {toString .} map] ff getAttributeList setAttributeList 
-    groebner 0 get /gg set
+    %(---) messagen ff getAttributeList message
+    ff getAttributeList tag 0 eq {/gbopt [ ] def }
+    {
+       /gbopt ff getAttributeList def
+    } ifelse
+   [ff {toString .} map gbopt] 
+    groebner 0 get /gg set   %% groenber は attribute を受け付けない.
     /cone.gb_Dh.g gg def
     /arg1 gg def
   ] pop
@@ -2361,7 +2398,7 @@ def   
 %  Usages: result_getNextFlip getNextCone ncone
 %  flip して新しい ncone を得る.
 %>
-/getNextCone {
+/getNextCone.orig {
  /arg1 set
  [/ncone /ccone /kk /w /next_weight_w_wv] pushVariables
  [ 
@@ -3430,6 +3467,7 @@ def   
 %  gb の check をやるので, それに失敗したら null を戻す.
 %  weight はすべて vw 形式で. vw 形式 = variable weight の繰り返しの形式
 %  reducedGb は文字列のリストではなく多項式の形式のこと.
+%   理由は reducedGb より ring の構造を読むため.
 %>
 /ckmFlip {
   /arg1 set
@@ -3476,6 +3514,12 @@ def   
 % 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2]
     facetWeight { dup isString { . rTable replace toString }
                                { } ifelse } map /facetWeight_gr set
+
+% newWeight も 新しい環 gr_ww  の weight に変換.
+% 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2]
+    newWeight { dup isString { . rTable replace toString }
+                               { } ifelse } map /newWeight_gr set
+
 % Dx x = x Dx + h H  or Dx x = x Dx + h^2 で計算.  
 % どちらをとるかは cone.gb_gr で区別するしかなし
     %% [ch1 vlist_gr oldWeight_gr] /ttt set 
@@ -3484,10 +3528,6 @@ def   
     ch1 {toString .} map  /ch1 set
 %% ここまででとりあえずテストをしよう.
 %%    ch1 /arg1 set
-% newWeight も 新しい環 gr_ww  の weight に変換.
-% 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2]
-    newWeight { dup isString { . rTable replace toString }
-                               { } ifelse } map /newWeight_gr set
     [ch1 { toString } map vlist_gr newWeight_gr] cone.gb_gr /ch2 set
 
 % Dx x = x Dx + h H  or Dx x = x Dx + h^2 で計算.  
@@ -3505,15 +3545,23 @@ def   
     ccf { 2 get {toString  . rTable2 replace toString} map } map /ccf2 set
     %% ccf2 は gr でない ring の元.
     gOld getRing ring_def
-    cone.beginH  % Hh か h^2 か.
+    cone.DhH { cone.begin_DhH } {  } ifelse % Hh か h^2 か.
     ccf2 { {.} map gOld mul } map /gNew set
     gNew { toString } map /gNew set
-    cone.endH
+    cone.DhH { cone.end_DhH } {  } ifelse % Hh か h^2 か.
     % gNew /arg1 set
     %gNew が newWeight での GB か check. Yes なら reduced basis へ.
     %No なら null を戻す. 
-    gNew [(gbCheck) 1] setAttributeList newWeight 
-       cone.gb (gb) getAttribute 
+%%Ref: note @s/2005/06/30-note-gfan.pdf
+    cone.do_gbCheck not { 
+       (Warning! gbCheck is skipped.) message
+    } { 
+       (Doing gbCheck.) message
+    } ifelse
+    cone.do_gbCheck {
+     gNew [(gbCheck) 1] setAttributeList newWeight 
+        cone.gb (gb) getAttribute 
+    } { 1 } ifelse
     1 eq {
      gNew [(reduceOnly) 1] setAttributeList newWeight cone.gb /arg1 set
     }{ /arg1 null def } ifelse
@@ -3524,30 +3572,122 @@ def   
 
 %<
 % Usages: f gbasis cone.reduction_DhH
+%       dx x = x dx + h H での reduction. 
 %>
 /cone.reduction_DhH {
   /arg2 set /arg1 set
   [/ff /ggbasis /eenv /ans] pushVariables
   [  
     /ff arg1 def /ggbasis arg2 def
-    cone.beginH
+    cone.begin_DhH
     ff ggbasis reduction /ans set
-    cone.endH
+    cone.end_DhH
     /arg1 ans def
   ] pop
   popVariables
   arg1
 } def
 
+%<
+% Usages: f gbasis cone.reduction_Dh
+%       dx x = x dx + h^2 での reduction. 
+%>
+/cone.reduction_Dh {
+  /arg2 set /arg1 set
+  [/ff /ggbasis /eenv /ans] pushVariables
+  [  
+    /ff arg1 def /ggbasis arg2 def
+    ff ggbasis reduction /ans set
+    /arg1 ans def
+  ] pop
+  popVariables
+  arg1
+} def
+
+%<
+% Usages: cone.begin_DhH   dx x = x dx + h H を開始.
+%>
 /cone.begin_DhH {
   [(Homogenize) (AutoReduce) (KanGBmessage)] pushEnv /cone.eenv set
   [(Homogenize) 3] system_variable
 } def
 
+%<
+% Usages: cone.begin_DhH   dx x = x dx + h H を終了.
+%>
 /cone.end_DhH {
   cone.eenv popEnv
 } def
 
+%<
+% Usages: ff vv ww cone.gb_gr_DhH   dx x = x dx + h H で計算.
+%   dh.gb は dhecart.sm1 で定義されており, dx x = x dx + h H での計算.
+%   gr をとっても, -w,w の場合は 微分作用素環のままであり, これが必要.
+%   bug? cone.gb で十分?
+%>
+/cone.gb_gr_DhH {
+  /arg1 set
+  [/ff /ww /vv] pushVariables
+  [
+     /ff arg1 0 get def
+     /vv arg1 1 get def
+     /ww arg1 2 get def
+     /dh.gb.verbose 1 def
+     /dh.autoHomogenize 0 def
+     [(AutoReduce) 1] system_variable
+     [ff { toString } map vv
+      [ww vv generateD1_1]] dh.gb 0 get /arg1 set
+  ] pop
+  popVariables
+  arg1 
+} def
+%<
+% Usages: ff vv ww cone.gb_gr_Dh   dx x = x dx + h^2 で計算.
+%   gb は dhecart.sm1 で定義されており, dx x = x dx + h^2 での計算.
+%   gr をとっても, -w,w の場合は 微分作用素環のままであり, これが必要.
+%   bug? cone.gb で十分?
+%>
+/cone.gb_gr_Dh {
+  /arg1 set
+  [/ff /ww /vv /gg /envtmp] pushVariables
+  [
+     /ff arg1 0 get def
+     /vv arg1 1 get def
+     /ww arg1 2 get def
+
+     [(AutoReduce) (KanGBmessage)] pushEnv /envtmp set
+     [(AutoReduce) 1] system_variable
+     [(KanGBmessage) 1] system_variable
+     [vv ring_of_differential_operators
+     [ww] weight_vector 0] define_ring
+     [ff {toString .} map] ff getAttributeList setAttributeList 
+     groebner 0 get /gg set
+     envtmp popEnv
+
+     /arg1 gg def
+  ] pop
+  popVariables
+  arg1 
+} def
+
+
+% これらは cone.ckmFlip 1 の時しか使わず.
+/cone.reduction {
+  cone.DhH {
+    cone.reduction_DhH
+  }{
+    cone.reduction_Dh
+  } ifelse
+} def
+/cone.gb_gr {
+  cone.DhH {
+    cone.gb_gr_DhH
+  }{
+    cone.gb_gr_Dh
+  } ifelse
+} def
+
+
 /test1.ckmFlip {
  % cf. cone.sample2
    cone.load.cohom
@@ -3561,6 +3701,10 @@ def   
  /cone.parametrizeWeightSpace {
    4 2 parametrizeSmallFan
  } def
+
+ /cone.DhH 1 def
+ /cone.ckmFlip 1 def
+
  /cone.local 1 def	
  /cone.w_start  null def   
  /cone.h0 1 def
@@ -3579,38 +3723,8 @@ def   
    cone.input { . homogenize toString } map /cone.input set
    dh.end
 
- /cone.gb {
-   cone.gb_DhH
- } def
 
- /cone.reduction {
-   cone.reduction_DhH
- } def
-
- /cone.beginH {
-   cone.begin_DhH
- } def
- /cone.endH {
-   cone.end_DhH
- } def
 % テストを開始する.
- /cone.gb_gr {
-  /arg1 set
-  [/ff /ww /vv] pushVariables
-  [
-     /ff arg1 0 get def
-     /vv arg1 1 get def
-     /ww arg1 2 get def
-     /dh.gb.verbose 1 def
-     /dh.autoHomogenize 0 def
-     [(AutoReduce) 1] system_variable
-     [ff { toString } map vv
-      [ww vv generateD1_1]] dh.gb 0 get /arg1 set
-  ] pop
-  popVariables
-  arg1 
- } def
-
 % getStartingCone /cone.ncone set
 % cone.ncone updateFan 
 % cone.gblist 0 get message
@@ -3625,4 +3739,108 @@ def   
  [ff (t1,t2,x,y) wOld wFacet wNew] ckmFlip /ff2 set
  (See ff and ff2) message  
 
+} def
+
+%<
+% Usages: cone i getaVectorOnFacet
+% cone の i 番目の facet の上の vector を求める.
+% cf. liftWeight
+%> 
+/getaVectorOnFacet {
+  /arg2 set /arg1 set
+  [/cone /facet_i /ep /vp /v /v /ii] pushVariables
+  [
+    /cone arg1 def /facet_i arg2 def
+    facet_i to_int32 /facet_i set
+  
+    cone (facetsv) getNode 2 get facet_i get /v set
+    /vp v 0 get def
+    1 1 v length 1 sub {
+      /ii set
+      vp v ii get  add /vp set
+    } for
+    vp nnormalize_vec /vp set
+    /arg1 vp def
+  ] pop 
+  popVariables
+  arg1 
+} def
+
+/getNextCone {
+  getNextCone_ckm
+} def
+
+%<
+%  Usages: result_getNextFlip getNextCone_ckm ncone
+%  flip して新しい ncone を得る.  Collar-Kalkbrener-Moll のアルゴリズムを使う
+%  if (cone.ckmFlip == 0) 普通の計算 else CKM.
+%>
+/getNextCone_ckm {
+ /arg1 set
+ [/ncone /ccone /kk /w /next_weight_w_wv /cid /ttt] pushVariables
+ [ 
+  /ccone arg1 def
+  /ncone null def
+  /kk ccone 1 get def  % kk は cid 番目の cone の kk 番目の facet を表す.
+  /cid ccone 2 get def % cid は cone の 番号.
+  ccone 0 get /ccone set
+  {
+   ccone tag 0 eq { exit } {  } ifelse
+
+% ccone の kk 番目の facet について flip する.
+   ccone kk cone.epsilon flipWeight  /w set
+   (Trying new weight is ) messagen w message
+   w liftWeight /next_weight_w_wv set
+   (Trying new weight [w,wv] is ) messagen next_weight_w_wv message
+
+   cone.ckmFlip {
+    [ 
+     cone.gblist cid get (grobnerBasis) getNode 2 get % reduce gb
+     cone.vv
+     cone.gblist cid get (weight) getNode [2 0 2] get % weight
+     ccone kk getaVectorOnFacet liftWeight 1 get  % weight on facet
+     next_weight_w_wv 1 get  % new weight  
+    ] /ttt set
+     ttt message
+     ttt ckmFlip /cone.cgb set
+   }{ 
+     cone.input next_weight_w_wv 1 get cone.gb /cone.cgb set
+   } ifelse
+
+  cone.cgb tag 0 eq not {
+   [w] next_weight_w_wv join /cone.cgb_weight set
+   next_weight_w_wv 1 get cone.cgb coneEq /cone.g_ineq set
+   cone.g_ineq cone.w_ineq join cone.Wt mul cone.Lpt mul 
+   pruneZeroVector /cone.gw_ineq_projectedWtLpt set
+
+   (cone.gw_ineq_projectedWtLpt is obtained.) message
+
+   cone.gw_ineq_projectedWtLpt getConeInfo /cone.nextConeInfo set
+% 次元を調べる.  だめなら retry
+   cone.nextConeInfo 0 get 0 get to_int32 cone.d eq {
+     cone.nextConeInfo 1 get newCone /ncone set
+     ccone ncone getCommonFacet 0 get {
+       (Flip succeeded.) message
+       exit
+     } { } ifelse
+   } { } ifelse
+% common face がなければ やはり epsilon を小さく.
+   cone.nextConeInfo 0 get 0 get to_int32 cone.d eq {
+    (ccone and ncone do not have a common facet.) message
+   } {  
+    (ncone is not maximal dimensional. ) message
+   } ifelse
+  }{ } ifelse
+
+   (Decreasing epsilon to ) messagen
+   cone.epsilon (1).. (2).. div mul /cone.epsilon set
+     cone.epsilon cone.epsilon.limit sub numerator (0).. lt {
+       (Too small cone.epsilon ) error
+     }  {  } ifelse
+   cone.epsilon message
+  } loop
+  /arg1 ncone def
+ ] pop
+ popVariables
+ arg1
 } def