=================================================================== RCS file: /home/cvs/OpenXM/src/kan96xx/Doc/gfan.sm1,v retrieving revision 1.1 retrieving revision 1.13 diff -u -p -r1.1 -r1.13 --- OpenXM/src/kan96xx/Doc/gfan.sm1 2004/09/05 10:19:29 1.1 +++ OpenXM/src/kan96xx/Doc/gfan.sm1 2009/08/26 04:54:17 1.13 @@ -1,8 +1,264 @@ -% $OpenXM$ +% $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.12 2005/07/07 07:53:37 takayama Exp $ % cp cone.sm1 $OpenXM_HOME/src/kan96xx/Doc/gfan.sm1 -% $Id: gfan.sm1,v 1.1 2004/09/05 10:19:29 takayama Exp $ +% $Id: gfan.sm1,v 1.13 2009/08/26 04:54:17 takayama Exp $ % iso-2022-jp +%%Ref: @s/2004/08/21-note.pdf +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Two examples are given below to get a global Grobner fan and +%% a local Grobner fan ; cone.sample and cone.sample2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Global Grobner Fan +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% How to input data? An example. (cf. test13.sm1) +%% Modify the following or copy the /cone.sample { ... } def +%% to your own file, +%% edit it, and execute it by " cone.sample ; " +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +/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 +] cat def + +% List of variables +% If cone.type=1, then (H) should be added. +/cone.vlist [(x11) (x12) (x13) (x21) (x22) (x23) + (Dx11) (Dx12) (Dx13) (Dx21) (Dx22) (Dx23) (h)] def + +% List of variables in the form for define_ring. +/cone.vv (x11,x12,x13,x21,x22,x23) def + +% If cone.type=0, then x,Dx, +% If cone.type=1, then x,Dx,h,H (Doubly homogenized) +% If cone.type=2, then x,Dx,h +/cone.type 2 def + +% Set how to parametrize the weight space. +% In the example below, 6 means the number of variables x11,x12,x13,x21,x22,x33 +% p q parametrizeSmallFan (p >= q) : Enumerate Grobner cones in the Small +% Grobner fan. +% The weights for the last p-q variables +% are 0. +% Example. 6 2 parametrizeSmallFan weights for x12,x21,x22,x23 are 0. +% +% p q parametrizeTotalFan (p = q = number of variables in cone.vv) +% p > q has not yet been implemented. +% +/cone.parametrizeWeightSpace { + 6 6 parametrizeSmallFan +} def + +% If you want to enumerate Grobner cones in local order (i.e., x^e <= 0), +% then cone.local = 1 else cone.local = 0. +/cone.local 0 def + +% Initial value of the weight in the weight space of which dimension is +% cone.m +% If it is null, then a random weight is used. +/cone.w_start + null +def + +% If cone.h0=1, then the weight for h is 0. +% It is usally set to 1. +/cone.h0 1 def + +% Set input polynomials which generate the ideal. +% Input must be homogenized. +% (see also data/test14.sm1 for double homogenization.) +/cone.input + [ + (x11 x22 - x12 x21) + (x12 x23 - x13 x22) + (x11 x23 - x13 x21) + ] +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. +% ( Computation in ^O and h[0,1](^D) need this +% as the first step. /cone.local 1 def ) +/cone.gb { + cone.gb_Dh +} def + + +cone.comment message +(cone.input = ) message +cone.input message +%%%% Step 1. Enumerating the Grobner Cones in a global ring. +%%%% The result is stored in cone.fan +getGrobnerFan + +%%%% If you want to print the output, then uncomment. +printGrobnerFan + +%%%% If you want to save the data to the file sm1out.txt, then uncomment. +% /cone.withGblist 1 def saveGrobnerFan /ff set ff output + +%%%% Step 2. Dehomogenize the Grobner Cones +%%%% by the equivalence relation in a local ring (uncomment). +% dhCones_h + +%%%% Generate the final data dhcone2.fan (a list of local Grobner cones.) +% dhcone.rtable + +%%%% Output dhcone2.fan with explanations +% dhcone.printGrobnerFan + +} def +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% End of " How to input data? An example. " +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Local Grobner Fan +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% How to input data? The example 2 (cf. test14.sm1). +%% Modify the following or copy the /cone.sample2 { ... } def +%% to your own file, +%% edit it, and execute if by " cone.sample2 ; " +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +/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 + (The Grobner cones are dehomogenized to get local Grobner fan.) nl +] cat def + +% List of variables +% If cone.type=1, then (H) should be added. +/cone.vlist [(t1) (t2) (x) (y) (Dt1) (Dt2) (Dx) (Dy) (h) (H)] def + +% List of variables in the form for define_ring. +/cone.vv (t1,t2,x,y) def + +% If cone.type=0, then x,Dx, +% If cone.type=1, then x,Dx,h,H (Doubly homogenized) +% If cone.type=2, then x,Dx,h +/cone.type 1 def + +% Set how to parametrize the weight space. +% In the example below, 6 means the number of variables x11,x12,x13,x21,x22,x33 +% p q parametrizeSmallFan (p >= q) : Enumerate Grobner cones in the Small +% Grobner fan. +% The weights for the last p-q variables +% are 0. +% Example. 6 2 parametrizeSmallFan weights for x12,x21,x22,x23 are 0. +% +% p q parametrizeTotalFan (p = q = number of variables in cone.vv) +% p > q has not yet been implemented. +% +/cone.parametrizeWeightSpace { + 4 2 parametrizeSmallFan +} def + +% If you want to enumerate Grobner cones in local order (i.e., x^e <= 0), +% then cone.local = 1 else cone.local = 0. +/cone.local 1 def + +% Initial value of the weight in the weight space of which dimension is +% cone.m +% If it is null, then a random weight is used. +/cone.w_start + null +def + +% If cone.h0=1, then the weight for h is 0. +% It is usally set to 1. +/cone.h0 1 def + +% Set input polynomials which generate the ideal. +% Input must be homogenized. +% (see also data/test14.sm1 for double homogenization.) +/cone.input + [ + (t1-y) (t2 - (y-(x-1)^2)) + ((-2 x + 2)*Dt2+Dx) + (Dt1+Dt2+Dy) + ] +def +% homogenize + [cone.vv ring_of_differential_operators + [[(t1) -1 (t2) -1 (Dt1) 1 (Dt2) 1]] ecart.weight_vector + 0] define_ring + dh.begin + 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. +% ( Computation in ^O and h[0,1](^D) need this +% as the first step. /cone.local 1 def ) +/cone.gb { + cone.gb_DhH +} def + +cone.comment message +(cone.input = ) message +cone.input message +%%%% Step 1. Enumerating the Grobner Cones in a global ring. +%%%% The result is stored in cone.fan +getGrobnerFan + +%%%% If you want to print the output, then uncomment. +printGrobnerFan + +%%%% If you want to save the data to the file sm1out.txt, then uncomment. +% /cone.withGblist 1 def saveGrobnerFan /ff set ff output + +%%%% Step 2. Dehomogenize the Grobner Cones +%%%% by the equivalence relation in a local ring (uncomment). +dhCones_h + +%%%% Generate the final data dhcone2.fan (a list of local Grobner cones.) +dhcone.rtable + +%%%% Output dhcone2.fan with explanations +dhcone.printGrobnerFan + +} def +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% End of " How to input data? The example 2. " +%% Do not touch below. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +[(parse) (cgi.sm1) pushfile] extension + +% If you use local polymake, then comment out. +% If you use the cgi/polymake on the net, then uncomment out. +%/doPolymake {doPolymake.OoHG} def (Using doPolymake.OoHG ) message +%/polymake.start {polymake.start.OoHG} def (Using polymake.start.OoHG ) message +%% Choose it automatically. +[(which) (polymake)] oxshell tag 0 eq { + (Polymake is not installed in this system.) message + /doPolymake {doPolymake.OoHG} def + (Using doPolymake.OoHG ) message + /polymake.start {polymake.start.OoHG} def + (Using polymake.start.OoHG ) message +} { (Local polymake will be used.) message } ifelse + /cone.debug 1 def /ox.k0.loaded boundp { @@ -10,6 +266,57 @@ [(parse) (ox.sm1) pushfile] extension } ifelse +/cone.load.cohom { + /cone.loaded boundp { } + { + [(parse) (cohom.sm1) pushfile] extension +% [(parse) (cone.sm1) pushfile] extension % BUG? cone.sm1 overrides a global + % in cohom.sm1? + [(parse) (dhecart.sm1) pushfile] extension + /cone.loaded 1 def + oxNoX + polymake.start ( ) message + } ifelse +} def + +%% Usages: cone.gb_DhH. h H (double homogenized) 用の GB. +%% dhecart.sm1 を load してあること. 入力は同次でないといけない. +%% [cone.vv ring_of_differential_operators +%% [[(t1) -1 (t2) -1 (Dt1) 1 (Dt2) 1]] ecart.weight_vector +%% 0] define_ring +%% dh.begin homogenize dh.end などの方法で同次化できる. +/cone.gb_DhH { + /arg2 set /arg1 set + [/ff /ww] pushVariables + [ + /ff arg1 def + /ww arg2 def + /dh.gb.verbose 1 def + /dh.autoHomogenize 0 def + [(AutoReduce) 1] system_variable + [ff { toString } map cone.vv + [ww cone.vv generateD1_1]] ff getAttributeList setAttributeList + dh.gb 0 get /arg1 set + ] pop + popVariables + arg1 +} def + +% +% cone.fan, cone.gblist に fan のデータがはいる. +% +%%%%<<<< 初期データの設定例. 日本語版 data/test13 より. <<<<<<<<<<<<<< +/cone.sample.test13.ja { + /cone.loaded boundp { } + { + [(parse) (cohom.sm1) pushfile] extension + [(parse) (cone.sm1) pushfile] extension + /cone.loaded 1 def + } ifelse +/cone.comment [ + (Toric ideal for 1-simplex x 2-simplex, in k[x]) nl +] cat def +%------------------Globals---------------------------------------- % Global: cone.type % どの exponents を取り出すのか指定する. % cf. exponents, gbext h や H も見るか? @@ -20,6 +327,70 @@ % Global: cone.local % cone.local: Local か? 1 なら local +/cone.local 0 def + + +% Global: cone.h0 +% cone.h0: 1 なら h の weight 0 での Grobner fan を計算する. +/cone.h0 1 def + +% --------------- 入力データ用大域変数の設定 -------------------------- +% +% cone.input : 入力多項式系 +/cone.input + [ + (x11 x22 - x12 x21) (x12 x23 - x13 x22) + (x11 x23 - x13 x21) + ] +def + +% cone.vlist : 全変数のリスト +/cone.vlist [(x11) (x12) (x13) (x21) (x22) (x23) + (Dx11) (Dx12) (Dx13) (Dx21) (Dx22) (Dx23) (h)] def + +% cone.vv : define_ring 形式の変数リスト. +/cone.vv (x11,x12,x13,x21,x22,x23) def + +% cone.parametrizeWeightSpace : weight 空間を parametrize する関数. +% 大域変数 cone.W , cone.Wpos もきまる. +/cone.parametrizeWeightSpace { + 6 6 parametrizeSmallFan +} def + +% cone.w_start : weight空間における weight の初期値. +% この値で max dim cone が得られないと random weight による サーチが始まる. +% random にやるときは null にしておく. +/cone.w_start + [9 8 5 4 5 6] +def + +% cone.gb : gb を計算する関数. +/cone.gb { + cone.gb_Dh +} def + + + +( ) message +cone.comment message +(cone.input = ) messagen cone.input message +(Type in getGrobnerFan) message +(Do clearGlobals if necessary) message +(printGrobnerFan ; saveGrobnerFan /ff set ff output ) message + +} def +%%%%%%>>>>> 初期データの設定例おわり >>>>>>>>>>>>>>>>>>>>>> + +% Global: cone.type +% どの exponents を取り出すのか指定する. +% cf. exponents, gbext h や H も見るか? +% 0 : x,y,Dx,Dy +% 1 : x,y,Dx,Dy,h,H +% 2 : x,y,Dx,Dy,h +/cone.type 2 def + +% Global: cone.local +% cone.local: Local か? 1 なら local /cone.local 1 def % Global: cone.h0 @@ -31,8 +402,36 @@ % 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 の満たす % 不等式制約を求める. @@ -88,7 +487,7 @@ ww to_int32 /ww set % univNum があれば int32 に直しておく. /ww2 ww weightv def % v-w 形式を 数字のベクトルに. (init 用) - /eqs [ ] def % 不等式系の係数 + /eqs null def % 不等式系の係数 /gsize g length def 0 1 gsize 1 sub { /i set @@ -101,20 +500,21 @@ % in_ww(f) > f_j となる項の処理. iterms 1 exps length 1 sub { /j set - eqs [expsTop exps j get sub] join /eqs set + expsTop exps j get sub eqs cons /eqs set % exps[0]-exps[j] を eqs へ格納していく. } for % in_ww(f) = f_j となる項の処理. [(exponents) f ww2 init cone.type] gbext /exps set % exps は in(f) 1 1 iterms 1 sub { /j set - eqs [exps j get expsTop sub] join /eqs set - eqs [expsTop exps j get sub] join /eqs set + exps j get expsTop sub eqs cons /eqs set + expsTop exps j get sub eqs cons /eqs set % exps[j]-exps[0], exps[0]-exps[j] を格納. % 結果的に (exps[j]-exps[0]).w = 0 となる. } for } { } ifelse } for + eqs listToArray reverse /eqs set /arg1 eqs def ] pop popVariables @@ -439,6 +839,15 @@ $It translates null to (0)..$ ]] putUsages +%< +% Usages: newVector.with-1 +% (-1).. で埋めたベクトルを作る. +%> +/newVector.with-1 { + newVector { pop (-1).. } map +} def + + % [2 0] lcm は 0 をもどすがいいか? --> OK. %< @@ -615,6 +1024,11 @@ def [ /ceq arg1 def ceq pruneZeroVector /ceq set + + ceq length 0 eq { + (Monomial ideal is not accepted as an input.) cone_ir_input + } { } ifelse + ceq genPo2 /ceq2 set % ceq2 は polymake.data(polymake.INEQUALITIES(...)) 形式 % polymake で ceq2 の次元の計算. @@ -720,6 +1134,8 @@ def polydata (FACETS) getNode 2 get 0 get to_univNum { nnormalize_vec} map /facets set [[ ] ] facets join shell rest removeFirstFromPolymake /facets set + facets length 0 eq + {(Internal error. Facet data is not obtained. See OpenXM_tmp.) error} { } ifelse % vertices は cone の上にあるので整数倍 OK. 正規かする. polydata (VERTICES) getNode 2 get 0 get to_univNum { nnormalize_vec} map /vertices set @@ -729,11 +1145,16 @@ def { nnormalize_vec } map /ineq set [[ ] ] ineq join shell rest removeFirstFromPolymake /ineq set +% nextcid, nextfid を加える. nextcid は nextConeId の略. となりの cone 番号. +% nextfid は nextFacetId の略. となりの cone の facet +% 番号. [(cone) [ ] [ [(facets) [ ] facets] arrayToTree [(flipped) [ ] facets length newVector null_to_zero] arrayToTree [(facetsv) [ ] facets vertices newCone_facetsv] arrayToTree + [(nextcid) [ ] facets length newVector.with-1 ] arrayToTree + [(nextfid) [ ] facets length newVector.with-1 ] arrayToTree [(vertices) [ ] vertices] arrayToTree [(inequalities) [ ] ineq] arrayToTree ] @@ -787,6 +1208,33 @@ def } def %< +% Usages: [gb weight] newConeGB +% gb と weight を tree 形式にして格納する. +%> +/newConeGB { + /arg1 set + [/gbdata /gg /ww /rr] pushVariables + [ + /gbdata arg1 def +% gb + gbdata 0 get /gg set +% weight + gbdata 1 get /ww set +% + [(coneGB) [ ] + [ + [(grobnerBasis) [ ] gg] arrayToTree + [(weight) [ ] [ww]] arrayToTree + [(initial) [ ] gg { ww 2 get weightv init } map ] arrayToTree + ] + ] arrayToTree /rr set + /arg1 rr def + ] pop + popVariables + arg1 +} def + +%< % Usages: cone_random %> /cone_random.start (2).. def @@ -1196,6 +1644,8 @@ def %< % Usages: pruneZeroVector % genPo, getConeInfo 等の前に使う. 0 ベクトルは意味のない制約なので除く. +% 同じ制約条件ものぞく. polymake FACET が正しく動かない場合があるので. +% cf. pear/OpenXM_tmp/x3y2.poly, x^3+y^2, x^2+y^3 data/test15.sm1 %> /pruneZeroVector { /arg1 set @@ -1203,6 +1653,7 @@ def [ /mm arg1 def mm to_univNum /mm set + [ [ ] ] mm join shell rest uniq /mm set [ 0 1 mm length 1 sub { /ii set @@ -1367,8 +1818,30 @@ def popVariables } def +%< +% Usages: cone i [cid fid] markNext +% cone の i 番目の facet のとなりの cone id (cid) と face id (fid) を設定する. +% cone の nextcid[i] = cid; nextfid[i] = fid となる. +% cone 自体が変更される. +% cone は class-tree. +%> +/markNext { + /arg3 set /arg2 set /arg1 set + [/cone /facet_i /vv /nextid] pushVariables + [ + /cone arg1 def /facet_i arg2 def /nextid arg3 def + facet_i to_int32 /facet_i set + cone (nextcid) getNode 2 get /vv set + vv facet_i , nextid 0 get to_univNum , put + cone (nextfid) getNode 2 get /vv set + vv facet_i , nextid 1 get to_univNum , put + ] pop + popVariables +} def + + %< % Usages: cone getNextFacet i % flipped の mark のない facet の index facet_i を戻す. @@ -1410,6 +1883,11 @@ def % note: 2004.9.2 cone (facetsv) getNode 2 get facet_i get /v set cone (facets) getNode 2 get facet_i get /f set + + v length 0 eq { + (The codimension of the linarity space of the Grobner cone seems to be 1 or 0.) cone_ir_input + } { } ifelse + /vp v 0 get def 1 1 v length 1 sub { /ii set @@ -1587,14 +2065,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] 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 @@ -1711,8 +2195,10 @@ def wv_start pmat %[3] reduced GB の計算. cone.input wv_start cone.gb /reduced_G set - (Reduced GB : ) message - reduced_G pmat + (Reduced GB is obtained: ) message + %reduced_G pmat + /cone.cgb reduced_G def + [cone.w_start w_start wv_start] /cone.cgb_weight set %[4] 射影してから polytope のデータを計算. wv_start reduced_G coneEq /cone.g_ineq set @@ -1731,7 +2217,8 @@ def cone.cinit 0 get 0 get to_int32 cone.m eq { exit } { (Failed to get the max dim cone. Updating the weight ...) messagen - /w_start cone.m cone_random_vec cone.W mul def + cone.m cone_random_vec /cone.w_start set + /w_start cone.w_start cone.W mul def % cone.cinit を再度計算するために clear する. /cone.cinit null def } ifelse @@ -1803,17 +2290,22 @@ def %> /markBorder { /arg1 set - [/cone /facets_t /flipped_t /kk] pushVariables + [/cone /facets_t /flipped_t /kk /nextcid_t /nextfid_t] pushVariables [ /cone arg1 def cone (facets) getNode 2 get /facets_t set cone (flipped) getNode 2 get /flipped_t set + cone (nextcid) getNode 2 get /nextcid_t set + cone (nextfid) getNode 2 get /nextfid_t set 0 1 flipped_t length 1 sub { /kk set flipped_t kk get (0).. eq { cone kk isOnWeightBorder { % Border の上にあるので flip 済のマークをつける. flipped_t kk (2).. put +% となりの cone の id (nextcid, nextfid) は -2 とする. + nextcid_t kk (-2).. put + nextfid_t kk (-2).. put } { } ifelse } { } ifelse } for @@ -1833,6 +2325,8 @@ def /cone.fan [ ] def % global: cone.incidence /cone.incidence [ ] def +% global: cone.gblist gb's standing for each cones in cone.fan. +/cone.gblist [ ] def /updateFan { /arg1 set @@ -1840,6 +2334,9 @@ def [ /ncone arg1 def /cone.fan.n cone.fan length def +% -1. cone.cgb (直前に計算された gb) と cone.cgb_weight (直前の計算の weight) +% を cone.gblist へ格納する. + cone.gblist [ [cone.cgb cone.cgb_weight] newConeGB ] join /cone.gblist set % 0. ncone が cone.fan にすでにあればエラー 0 1 cone.fan.n 1 sub { /kk set @@ -1865,6 +2362,9 @@ def ncone ii markFlipped cone.fan kk get /tcone set tcone jj markFlipped +% nextcid, nextfid を設定する. + ncone ii [kk jj] markNext + tcone jj [cone.fan.n ii] markNext } { } ifelse } for % 3. ncone を加える. @@ -1874,22 +2374,25 @@ def } def %< -% usages: getNextFlip [cone, k] +% usages: getNextFlip [cone, k, cid] % cone.fan を検索して まだ flip してない cone と facet の組を戻す. % もうないときには null を戻す. +% cid は cone が cone.fan の 何番目であるかの index. cone.gblist の検索等に +% 用いる. %> /getNextFlip { - [/tcone /ans /ii ] pushVariables + [/tcone /ans /ii /cid] pushVariables [ - /ans null def + /ans null def /cid -1 def 0 1 cone.fan length 1 sub { /ii set cone.fan ii get /tcone set + /cid ii def tcone getNextFacet /ans set ans tag 0 eq { } { exit } ifelse } for ans tag 0 eq { /arg1 null def } - { /arg1 [tcone ans] def } ifelse + { /arg1 [tcone ans cid] def } ifelse ] pop popVariables arg1 @@ -1899,12 +2402,13 @@ def % flip の時の epsilon /cone.epsilon (1).. (10).. div def /cone.epsilon.limit (1).. (100).. div def +% cone.epsilon.limit を負にすれば停止しない. %< % Usages: result_getNextFlip getNextCone ncone % flip して新しい ncone を得る. %> -/getNextCone { +/getNextCone.orig { /arg1 set [/ncone /ccone /kk /w /next_weight_w_wv] pushVariables [ @@ -1922,6 +2426,7 @@ def (Trying new weight [w,wv] is ) messagen next_weight_w_wv message cone.input next_weight_w_wv 1 get cone.gb /cone.cgb set + [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 @@ -1975,5 +2480,1394 @@ def cone.nextflip tag 0 eq { exit } { } ifelse cone.nextflip getNextCone /cone.ncone set } loop - (Construction is completed. See cone.fan and cone.incidence.) message -} def \ No newline at end of file + (Construction is completed. See cone.fan, cone.incidence and cone.gblist.) + message +} def + +%< +% Usages: vlist generateD1_1 +% -1,1 weight を生成する. +% vlist は (t,x,y) か [(t) (x) (y)] +% +%> +/generateD1_1 { + /arg1 set + [/vlist /rr /rr /ii /vv] pushVariables + [ + /vlist arg1 def + vlist isString { + [vlist to_records pop] /vlist set + } { } ifelse + [ + 0 1 vlist length 1 sub { + /ii set + vlist ii get /vv set + vv -1 + [@@@.Dsymbol vv] cat 1 + } for + ] /rr set + /arg1 rr def + ] pop + popVariables + arg1 +} def + +/listNodes { + /arg1 set + [/in-listNodes /ob /rr /rr /ii] pushVariables + [ + /ob arg1 def + /rr [ ] def + { + ob isClass { + ob (array) dc /ob set + } { exit } ifelse + rr [ob 0 get] join /rr set + ob 2 get /ob set + 0 1 ob length 1 sub { + /ii set + rr ob ii get listNodes join /rr set + } for + exit + } loop + /arg1 rr def + ] pop + popVariables + arg1 +} def +[(listNodes) +[(ob listNodes) + (cf. getNode) + (Example:) + ( /dog [(dog) [[(legs) 4] ] [ ]] [(class) (tree)] dc def) + ( /man [(man) [[(legs) 2] ] [ ]] [(class) (tree)] dc def) + ( /ma [(mammal) [ ] [man dog]] [(class) (tree)] dc def) + ( ma listNodes ) +]] putUsages + +%< +% Usages: obj printTree +%> +/printTree { + /arg1 set + [/ob /rr /rr /ii /keys /tt] pushVariables + [ + /ob arg1 def + /rr [ ] def + /keys ob listNodes def + keys 0 get /tt set + keys rest /keys set + keys { ob 2 1 roll getNode } map /rr set + (begin ) messagen tt messagen + ( ---------------------------------------) message + 0 1 rr length 1 sub { + /ii set + keys ii get messagen (=) message + rr ii get 2 get pmat + } for + (--------------------------------------- end ) messagen + tt message + /arg1 rr def + ] pop + popVariables + arg1 +} def + +%< +% Usages は (inputForm) usages をみよ. +%> +/inputForm { + /arg1 set + [/ob /rr /i ] pushVariables + [ + /ob arg1 def + /rr [ ] def + { + ob isArray { + rr [ ([) ] join /rr set + 0 1 ob length 1 sub { + /i set + i ob length 1 sub lt { + rr [ob i get inputForm $ , $] join /rr set + } { + rr [ob i get inputForm] join /rr set + } ifelse + } for + rr [ (]) ] join cat /rr set + exit + } { } ifelse + ob isClass { + ob etag 263 eq { % tree + /rr ob inputForm.tree def exit + } { /rr [( $ this etag is not implemented $ )] cat def exit } ifelse + } { } ifelse + ob isUniversalNumber { + [$($ ob toString $)..$] cat /rr set + exit + } { } ifelse + ob isPolynomial { + [$($ ob toString $).$] cat /rr set + exit + } { } ifelse + ob isRational { + [$ $ ob (numerator) dc inputForm $ $ + ob (denominator) dc inputForm $ div $ ] cat /rr set + exit + } { } ifelse + ob isString { + [$($ ob $)$ ] cat /rr set + exit + } { } ifelse + ob toString /rr set + exit + } loop + rr /arg1 set + ] pop + popVariables + arg1 +} def +[(inputForm) + [(obj inputForm str) +]] putUsages +% should be moved to dr.sm1 + +/inputForm.tree { + /arg1 set + [/ob /key /rr /rr /ii] pushVariables + [ + /ob arg1 def + /rr [ ] def + { + ob (array) dc /ob set + /rr [ $[$ ob 0 get inputForm $ , $ + ob 1 get inputForm $ , $ + ] def + rr [ob 2 get inputForm ] join /rr set + rr [$ ] $] join /rr set + rr [ $ [(class) (tree)] dc $ ] join /rr set + rr cat /rr set + exit + } loop + /arg1 rr def + ] pop + popVariables + arg1 +} def + +%< +% Usages: str inputForm.value str +%> +/inputForm.value { + /arg1 set + [/key /val /valstr /rr] pushVariables + [ + arg1 /key set + key isString { } {(inputForm.value: argument must be a string) error } ifelse + key boundp { + [(parse) key] extension pop + /val set + val inputForm /valstr set + [( ) valstr ( /) key ( set )] cat /rr set + } { + /valstr [] cat /rr set + } ifelse + rr /arg1 set + ] pop + popVariables + arg1 +} def + +% global: cone.withGblist +/cone.withGblist 0 def +%< +% Usages: saveGrobnerFan str +% GrobnerFan のデータを inputForm に変更して文字列に変える. +% このデータを parse すると GrobnerFan を得ることが可能. +% BUG: 多項式の属する環のデータの保存はまだしてない. +%> +/saveGrobnerFan { + [/rr] pushVariables + [ + (cone.withGblist=) messagen cone.withGblist message + [ +% ユーザの設定するパラメータ. cone.gb, cone.parametrizeWeightSpace 等の関数もあり. + (cone.comment) + (cone.type) (cone.local) (cone.h0) + (cone.vlist) (cone.vv) + (cone.input) + +% プログラム中で利用する, 大事な大域変数. weight vector の射影行列が重要. + (cone.n) (cone.m) (cone.d) + (cone.W) (cone.Wpos) (cone.Wt) + (cone.L) (cone.Lp) (cone.Lpt) + (cone.weightBorder) + (cone.w_ineq) + (cone.w_ineq_projectedWt) + (cone.epsilon) + +% 結果の要約. + (cone.fan) + cone.withGblist { (cone.gblist) } { } ifelse + (cone.incidence) + + ] { inputForm.value nl } map /rr set + rr cat /rr set +% ring を save してないので当座の対処. + [ ([) cone.vv inputForm ( ring_of_differential_operators 0 ] define_ring ) + nl nl rr] cat /arg1 set + ] pop + popVariables + arg1 +} def + +/printGrobnerFan.1 { + /arg1 set + [/key /rr] pushVariables + [ + /key arg1 def + key boundp { + [(parse) key] extension pop /rr set + rr isArray { + key messagen ( = ) message rr pmat + } { + key messagen ( = ) messagen rr message + } ifelse + }{ + key messagen ( = ) message + } ifelse + ] pop + popVariables +} def + +/printGrobnerFan { + [/i] pushVariables + [ + (========== Grobner Fan ====================) message + [ + (cone.comment) + (cone.vlist) (cone.vv) + (cone.input) + (cone.type) (cone.local) (cone.h0) + (cone.n) (cone.m) (cone.d) + (cone.W) (cone.Wpos) (cone.Wt) + (cone.L) (cone.Lp) (cone.Lpt) + (cone.weightBorder) + (cone.incidence) + ] { printGrobnerFan.1 } map + ( ) message + 0 1 cone.fan length 1 sub { + /ii set + ii messagen ( : ) messagen + cone.fan ii get printTree + } for + cone.withGblist { + 0 1 cone.gblist length 1 sub { + /ii set + ii messagen ( : ) messagen + cone.gblist ii get printTree + } for + } { } ifelse + + + (=========================================) message + (cone.withGblist = ) messagen cone.withGblist message + ( ) message + ] pop + popVariables +} def + +%< +% Usages: m uniq +% Remove duplicated lines. +%> +/uniq { + /arg1 set + [/mm /prev /i /rr] pushVariables + [ + /mm arg1 def + { + mm length 0 eq { [ ] /rr set exit } { } ifelse + /prev mm 0 get def + [ + prev + 1 1 mm length 1 sub { + /i set + mm i get prev sub isZero { } + { /prev mm i get def prev } ifelse + } for + ] /rr set + exit + } loop + rr /arg1 set + ] pop + popVariables + arg1 +} def + +%< +% Usages: [vlist vw_vector] getGrRing [vlist vGlobal sublist] +% example: [(x,y,z) [(x) -1 (Dx) 1 (y) 1 (Dy) 2]] getGrRing +% [(x,y,z,y') [(x)] [[(Dy) (y')]]] +% h[0,1](D_0) 専用の getGrRing. +% u_i + v_i > 0 なら Dx_i ==> x_i' (可換な変数). sublist へ. +% u_i < 0 なら x_i は vGlobal へ. +% ii [vlist vGlobal sublist] toGrRing /ii set +% [ii jj vlist [(partialEcartGlobalVarX) vGlobal]] ecart.isSameIdeal と使う. +%> +/getGrRing { + /arg1 set + [/vlist /vw_vector /ans /vGlobal /sublist /newvlist + /dlist /tt /i /u /v /k + ] pushVariables + [ + /vlist arg1 0 get def + /vw_vector arg1 1 get def + + vlist isString { [vlist to_records pop] /vlist set } { } ifelse + vlist { toString } map /vlist set +% dlist は [(Dx) (Dy) (Dz)] のリスト. + vlist { /tt set [@@@.Dsymbol tt] cat } map /dlist set + + /newvlist [ ] def /sublist [ ] def /vGlobal [ ] def +% 可換な新しい変数を newvlist へ. 置換表を sublist へ. + 0 1 vlist length 1 sub { + /i set +% (u,v) は (x_i, Dx_i) に対する weight vector + /u vlist i get , vw_vector getGrRing.find def + u -1 gt { + vw_vector , u 1 add , get /u set + } { /u 0 def } ifelse + + /v dlist i get , vw_vector getGrRing.find def + v -1 gt { + vw_vector , v 1 add , get /v set + } { /v 0 def } ifelse + u to_int32 /u set , v to_int32 /v set + + u v add , 0 gt { + newvlist [vlist i get] join /newvlist set + } { } ifelse + u 0 lt { + vGlobal [vlist i get] join /vGlobal set + } { } ifelse + } for + + newvlist { /tt set [ [@@@.Dsymbol tt] cat [tt (')] cat ] } map + /sublist set + + /ans [ vlist , newvlist { /tt set [tt (')] cat } map , join from_records + vGlobal sublist] def + /arg1 ans def + ] pop + popVariables + arg1 +} def + +%< +% Usages: a uset getGrRing.find index +%> +/getGrRing.find { + /arg2 set /arg1 set + [/a /uset /ans /i] pushVariables + [ + /a arg1 def /uset arg2 def + /ans -1 def + { /ans -1 def + 0 1 , uset length 1 sub { + /i set + a tag , uset i get tag eq { + a , uset i get eq { + /ans i def exit + } { } ifelse + } { } ifelse + } for + exit + } loop + /arg1 ans def + ] pop + popVariables + arg1 +} def + +%< +% Usages: g1 g2 isSameGrRing bool +% g1, g2 は getGrRing の戻り値. +%> +/isSameGrRing { + /arg2 set /arg1 set + [/g1 /g2 /ans] pushVariables + [ + /g1 arg1 def /g2 arg2 def + { + /ans 1 def + g1 0 get , g2 0 get eq { } { /ans 0 def exit } ifelse + exit + g1 1 get , g2 1 get eq { } { /ans 0 def exit } ifelse + } loop + /arg1 ans def + ] pop + popVariables + arg1 +} def + +%< +% Usages: [[ii i_vw_vector] [jj j_vw_vector] vlist] isSameInGrRing_h +% It computes gb. +%> +/isSameInGrRing_h { + /arg1 set + [/ii /i_vw_vector /jj /j_vw_vector /vlist + /i_gr /j_gr /rrule /ans] pushVariables + [ + /ii arg1 [0 0] get def + /i_vw_vector arg1 [0 1] get def + /jj arg1 [1 0] get def + /j_vw_vector arg1 [1 1] get def + /vlist arg1 2 get def + { + [vlist i_vw_vector] getGrRing /i_gr set + [vlist j_vw_vector] getGrRing /j_gr set + i_gr j_gr isSameGrRing { } { /ans [0 [i_gr j_gr]] def exit} ifelse + +% bug: in case of module + [i_gr 0 get , ring_of_differential_operators 0] define_ring + +% H を 1 に. + /rrule [ [@@@.Hsymbol . (1).] ] def + + i_gr 2 get length 0 eq { + } { + rrule i_gr 2 get { { . } map } map join /rrule set + } ifelse + ii { toString . rrule replace toString } map /ii set + jj { toString . rrule replace toString } map /jj set + + [ii jj i_gr 0 get , i_gr 1 get] ecartd.isSameIdeal_h /ans set + [ans [i_gr] rrule ecartd.isSameIdeal_h.failed] /ans set + + exit + } loop + /arg1 ans def + ] pop + popVariables + arg1 +} def + +/test1.isSameInGrRing_h { + [(parse) (data/test8-data.sm1) pushfile] extension + + cone.gblist 0 get (initial) getNode 2 get /ii set + cone.gblist 0 get (weight) getNode [2 0 2] get /iiw set + + cone.gblist 1 get (initial) getNode 2 get /jj set + cone.gblist 1 get (weight) getNode [2 0 2] get /jjw set + + (Doing [ [ii iiw] [jj jjw] cone.vv ] isSameInGrRing_h /ff set) message + [ [ii iiw] [jj jjw] cone.vv ] isSameInGrRing_h /ff set + + ff pmat + +} def + + +%< +% Usages: i j isSameCone_h.0 [bool, ...] +% テスト方法. (data/test8.sm1) run (data/test8-data.sm1) run 0 1 isSameCone_h.0 +% gb を再度計算する stand alone 版. gr(Local ring) で比較. +%> +/isSameCone_h.0 { + /arg2 set /arg1 set + [/i /j /ans /ii /iiw /jj /jjw] pushVariables + [ + /i arg1 def /j arg2 def + i to_int32 /i set , j to_int32 /j set + cone.debug { (Comparing ) messagen [i j] message } { } ifelse + + cone.gblist i get (initial) getNode 2 get /ii set + cone.gblist i get (weight) getNode [2 0 2] get /iiw set + + cone.gblist j get (initial) getNode 2 get /jj set + cone.gblist j get (weight) getNode [2 0 2] get /jjw set + + [ [ii iiw] [jj jjw] cone.vv ] isSameInGrRing_h /ans set + + ans /arg1 set + ] pop + popVariables + arg1 +} def + +%< +% Usages: [ii vv i_vw_vector] getGbInGrRing_h [ii_gr i_gr] +% Get Grobner Basis of ii in the graded ring. +% The graded ring is obtained automatically from vv and i_vw_vector. +% ii_gr is the Grobner basis. i_gr is the output of getGrRing. +% cf. isSameInGrRing_h, ecart.isSameIdeal_h with [(noRecomputation) 1] +%> +/getGbInGrRing_h { + /arg1 set + [/ii /i_vw_vector /vlist /rng /vv /vvGlobal /wv /iigg + /i_gr /rrule /ans] pushVariables + [ + /ii arg1 0 get def + /vlist arg1 1 get def + /i_vw_vector arg1 2 get def + [vlist i_vw_vector] getGrRing /i_gr set + +% bug: in case of module + [i_gr 0 get , ring_of_differential_operators 0] define_ring + +% H を 1 に. + /rrule [ [@@@.Hsymbol . (1).] ] def + + i_gr 2 get length 0 eq { + } { + rrule i_gr 2 get { { . } map } map join /rrule set + } ifelse + /vvGlobal i_gr 1 get def + /vv i_gr 0 get def + + ii { toString . rrule replace toString } map /ii set + + [vv vvGlobal] ecart.stdBlockOrder /wv set + vvGlobal length 0 eq { + /rng [vv wv ] def + }{ + /rng [vv wv [(partialEcartGlobalVarX) vvGlobal]] def + } ifelse + /save-cone.autoHomogenize ecart.autoHomogenize def + /ecart.autoHomogenize 0 def + [ii] rng join ecartd.gb /iigg set + save-cone.autoHomogenize /ecart.autoHomogenize set + /ans [iigg 0 get i_gr] def + /arg1 ans def + ] pop + popVariables + arg1 +} def + +/test1.getGbInGrRing_h { + [(parse) (data/test8-data.sm1) pushfile] extension + + cone.gblist 0 get (initial) getNode 2 get /ii set + cone.gblist 0 get (weight) getNode [2 0 2] get /iiw set + [ii cone.vv iiw] getGbInGrRing_h /ff1 set + + cone.gblist 1 get (initial) getNode 2 get /jj set + cone.gblist 1 get (weight) getNode [2 0 2] get /jjw set + [jj cone.vv jjw] getGbInGrRing_h /ff2 set + + (ff1 and ff2) message + +} def + + +%< +% setGrGblist +% cone.grGblist を設定する. +%> +/setGrGblist { + [/ii /ww /gg] pushVariables + [ + cone.gblist { + /gg set + gg (initial) getNode 2 get /ii set + gg (weight) getNode [2 0 2] get /ww set + [ii cone.vv ww] getGbInGrRing_h + } map /cone.grGblist set + ] pop + popVariables +} def + +%< +% Usages: i j isSameCone_h.2 [bool, ...] +% gb を再度計算しない. +%> +/isSameCone_h.2 { + /arg2 set /arg1 set + [/i /j /ans /ii /iiw /jj /jjw] pushVariables + [ + /i arg1 def /j arg2 def + i to_int32 /i set , j to_int32 /j set + (cone.grGblist) boundp { } { setGrGblist } ifelse + cone.debug { (Comparing ) messagen [i j] message } { } ifelse + + cone.grGblist i get /ii set + cone.grGblist j get /jj set + + ii 1 get , jj 1 get isSameGrRing { } + { /ans [0 [ii 1 get jj 1 get]] def exit} ifelse + + [ii 0 get , jj 0 get cone.vv [[(noRecomputation) 1]] ] + ecartd.isSameIdeal_h /ans set + [ans [ii 1 get] ii 1 get , ecartd.isSameIdeal_h.failed] /ans set + + ans /arg1 set + ] pop + popVariables + arg1 +} def + +%< +% test1.isSameCone_h.2 は cone.grGblist に initial の gb を graded ring +% でまず計算し, それから ideal の比較をおこなう. isSameCone_h.1 に比べて +% gb の再度の計算がないので経済的. +%> +/test1.isSameCone_h.2 { + /cone.loaded boundp { } + { + [(parse) (cohom.sm1) pushfile] extension + [(parse) (dhecart.sm1) pushfile] extension + /cone.loaded 1 def + } ifelse + %[(parse) (cone.sm1) pushfile] extension + [(parse) (data/test8-data.sm1) pushfile] extension + setGrGblist + (cone.grGblist is set.) message + 0 1 isSameCone_h.2 pmat +} def + +%< +% dhcone は DeHomogenized Cone の略. H->1 として cone を merge していく関数 +% や大域変数に使う. +% cone.gblist, cone.fan が正しく設定されていること. +% (setGrGblist を実行済であること. 自動実行されるが... ) +% +%> + +/isSameCone_h { isSameCone_h.2 } def + +%< +% Usages: genDhcone.init +% dhcone.checked (dehomogenized 済の cone番号), dhcone.unchecked の初期化. +%> +/genDhcone.init { + /dhcone.checked [ ] def + /dhcone.unchecked [ + 0 1 cone.fan length 1 sub { + to_univNum + } for + ] def +} def + +%< +% Usages: k genDhcone dhcone +% cone.fan[k] を出発点として cone を dehomogenize する (merge する). +% +% テスト1. (data/test14.sm1) run (data/test14-data.sm1) run +% genDhcone.init +% 0 genDhcone /ff set +%> + +/genDhcone { + /arg1 set + [/k /facets /merged /nextcid /nextfid /coneid + /newfacets /newmerged /newnextcid /newnextfid /newconeid /vv + /i /j /p /q /rr /cones /differentC + ] pushVariables + [ + /k arg1 def + /facets [ ] def /merged [ ] def /nextcid [ ] def + /nextfid [ ] def /coneid [ ] def + /cones [ ] def + /differentC [ ] def + + k to_univNum /k set + + { +% Step1. cone.fan[k] を 加える. new... へ初期データを書き込む. + cone.debug {(Step 1. Adding ) messagen k messagen (-th cone.) message} { } ifelse + cones [k to_univNum] join /cones set + cone.fan k get , (facets) getNode 2 get /vv set + /newfacets [ ] vv join def + + cone.fan k get , (nextcid) getNode 2 get /vv set + /newnextcid [ ] vv join def + + cone.fan k get , (nextfid) getNode 2 get /vv set + /newnextfid [ ] vv join def + +% newmerged はまず 0 でうめる. 0 : まだ調べてない. +% 1 : merged で消えた. 2 : boundary. 3 : となりは異なる. +% [ ] join をやって ベクトルの clone を作る. + cone.fan k get , (flipped) getNode 2 get /vv set + /newmerged [ ] vv join def + 0 1 , newmerged length 1 sub { + /i set + newmerged i get , (2).. eq { } + { newmerged i (0).. put } ifelse + } for +% newconeid は k でうめる. + /newconeid newfacets length newVector { pop k to_univNum } map def + +% merged と newmerged を cone の隣接関係のみで更新する. +% 同じ init を持つことはわかっているので facet vector のみの check で十分. +% merged の i 番目 と newmerged の j 番目で比較. + 0 1 , merged length 1 sub { + /i set + 0 1 , newmerged length 1 sub { + /j set + merged i get , (0).. eq , + newmerged j get , (0).. eq , and + nextcid i get , k to_univNum eq , and + { + facets i get , newfacets j get , add isZero { +% merged[i], newmerged[j] に 1 を入れて消す. +% 上の判定は nextfid, newnextfid を用いてもよいのでは? + merged i (1).. put + newmerged j (1).. put + } { } ifelse + } { } ifelse + } for + } for + +% Step2. 結合してから, まだ調べてない facet を探す. + cone.debug { (Step 2. Joining *** and new***) message } { } ifelse + /facets facets newfacets join def + /merged merged newmerged join def + /nextcid nextcid newnextcid join def + /nextfid nextfid newnextfid join + /coneid coneid newconeid join def + + cone.debug{ ( Checking facets.) message } { } ifelse + /k null def + 0 1 , merged length 1 sub { + /i set + % i message + merged i get (0).. eq { +% i 番目をまだ調べていない. + coneid i get , /p set + nextcid i get , /q set + cone.debug { [p q] message } { } ifelse + q (0).. ge { +% cone.fan [p] と cone.fan [q] の initial を比較する. +% 同じなら k を設定. exit for. 違えば merged[i] = 3 (違う) を代入. +% differentC はすでに 現在の dhcone と違うと check された cone 番号. +% dhcone.checked は dhcone がすでに生成されている cone 番号のリスト. +% これにはいっていても違う. + q differentC memberQ , q dhcone.checked memberQ , or + { /rr [0 ] def } + { p q isSameCone_h /rr set } ifelse + + rr 0 get 1 eq { + cone.debug { (Found next cone. ) message } { } ifelse + /k q to_univNum def exit + } { + cone.debug { ( It is a different cone. ) message } { } ifelse + differentC [ q ] join /differentC set + merged i (3).. put + } ifelse + } { } ifelse + } { } ifelse + } for + + k tag 0 eq { exit } { } ifelse + } loop + + [(-1)..] cones join shell rest /cones set +% dhcone.checked, dhcone.unchecked を更新. + dhcone.checked cones join /dhcone.checked set + dhcone.unchecked cones setMinus /dhcone.unchecked set + + [(dhcone) [ ] + [ + [(cones) [ ] cones] arrayToTree + [(facets) [ ] facets] arrayToTree + [(merged) [ ] merged] arrayToTree + [(nextcid) [ ] nextcid] arrayToTree + [(nextfid) [ ] nextfid] arrayToTree + [(coneid) [ ] coneid] arrayToTree + ] + ] arrayToTree /arg1 set + ] pop + popVariables + arg1 +} def + + +%< +% Usages: dhCones_h +% cone.fan は doubly homogenized (local) で生成された Grobner fan. +% cone.fan を dehomogenize (H->1) して init を比べて dhcone.fan を生成する. +% +% テスト1. (data/test14.sm1) run (data/test14-data.sm1) run +% dhCones_h +% test22 +%> +/dhCones_h { + (cone.grGblist) boundp { } {setGrGblist} ifelse + genDhcone.init + /dhcone.fan [ ] def + { + (-----------------------------------------) message + (#dhcone.unchecked = ) messagen dhcone.unchecked length message + dhcone.unchecked length 0 eq { exit } { } ifelse + dhcone.fan + [ dhcone.unchecked 0 get , genDhcone ] join /dhcone.fan set + (#dhcone.fan = ) messagen dhcone.fan length message + } loop + dhcone.fan +} def + +%< +% Usages: dhcone.rtable +% dhcone の番号と cone の番号の 置換表を生成し dhcone2.fan (merge した cone の情報) +% を dhcone.fan から作る. dhcone2.gblist も作る補助関数. +% dhCones_h してから dhcone.rable する. +%> +/dhcone.rtable { + [/i /j /vv /cones /facets /facets2 /merged /nextcid /nextcid2 /ii /ww] pushVariables + [ +% 置換表 dhcone.h2dh を作る. + /dhcone.h2dh cone.fan length newVector.with-1 def + 0 1 , dhcone.fan length 1 sub { + /i set + dhcone.fan i get , (cones) getNode 2 get /vv set + 0 1 vv length 1 sub { + /j set + dhcone.h2dh , vv j get , i to_univNum , put + } for + } for +% merge した dhcone を整理したもの, dhcone2.fan を作る. + /dhcone2.fan dhcone.fan length newVector def + 0 1 , dhcone.fan length 1 sub { + /i set + dhcone.fan i get (facets) getNode 2 get /facets set + dhcone.fan i get (merged) getNode 2 get /merged set + dhcone.fan i get (nextcid) getNode 2 get /nextcid set + dhcone.fan i get (cones) getNode 2 get /cones set + /facets2 [ ] def + /nextcid2 [ ] def + 0 1 , facets length 1 sub { + /j set + merged j get , (3).. eq { + facets2 [ facets j get ] join /facets2 set +% となりの cone があるとき 変換表にしたがい, cone 番号を変換 + nextcid2 [ dhcone.h2dh , nextcid j get , get ] join /nextcid2 set + } { } ifelse + merged j get , (2).. eq { + facets2 [ facets j get ] join /facets2 set +% 境界のとき -2 を入れる. + nextcid2 [ (-2).. ] join /nextcid2 set + } { } ifelse + } for + + dhcone2.fan i , + [(dhcone) [ ] + [ + [(facets) [ ] facets2] arrayToTree + [(nextcid) [ ] nextcid2] arrayToTree + [(cones) [ ] cones] arrayToTree + ] + ] arrayToTree , put + + } for + +% 最後に dhcone2.gblist を作る. + /dhcone2.gblist , dhcone2.fan length newVector , def + 0 1 , dhcone2.fan length 1 sub { + /i set + dhcone2.fan i get (cones) getNode 2 get /cones set + cone.grGblist , cones 0 get , get , /ii set % GB of initial (H->1). + cone.gblist i get , (weight) getNode , [ 2 0 2 ] get /ww set + + dhcone2.gblist i, + [(gbasis) [ ] + [ + [(initial) [ ] ii] arrayToTree + [(weight) [ ] ww] arrayToTree + ] + ] arrayToTree , put + + } for + (dhcone2.fan, dhcone2.gblist, dhcone.h2dh are set.) message + + ] pop + popVariables +} def + +%< +% 表の見方の解説を印刷する関数. +% Usages: dhcone.explain +%> +/dhcone.explain { + [ + ( ) nl + (Data format in << dhcone2.fan >>, which is a dehomogenized Grobner fan.) nl nl + (<< cone.vlist >> is the list of the variables.) nl + @@@.Hsymbol ( is the homogenization variable to be dehomogenized.) nl nl + (<< cone.input >> is generators of a given ideal.) nl nl + (<< cone.d >> is the dimension of parametrization space of the weights P_w) nl + ( P_w is a cone in R^m where the number m is stored in << cone.m >>) nl + ( P_w --- W ---> R^n [weight space]. ) nl + ( W is stored in << cone.W >> ) nl + ( << u cone.W mul >> gives the weight vector standing for u) nl nl + (All cones in the data lie in the weight parametrization space P_w.) nl + ( "facets" are the inner normal vector of the cone. ) nl + ( "nextcid" is a list of the cone id's of the adjacent cones.) nl + ( -2 in "nextcid" means that this facet lies on the border of the weight space.) nl + ( "cones" is a list of the cone id's of the NON-dehomonized Grobner fan) nl + ( stored in << cone.fan >>) nl + ] cat +} def + +%< +% dhcone.printGrobnerFan +% dhcone の印刷関数 +%> +/dhcone.printGrobnerFan { + [/i] pushVariables + [ + (========== Grobner Fan (for dehomogenized cones) ============) message + [ + (cone.comment) + (cone.vlist) (cone.vv) + (cone.input) + (cone.type) (cone.local) (cone.h0) + (cone.n) (cone.m) (cone.d) + (cone.W) (cone.Wpos) (cone.Wt) + (cone.L) (cone.Lp) (cone.Lpt) + (cone.weightBorder) + (cone.incidence) + ] { printGrobnerFan.1 } map + ( ) message + (The number of cones = ) messagen dhcone.fan length message + ( ) message + 0 1 dhcone2.fan length 1 sub { + /ii set + ii messagen ( : ) messagen + dhcone2.fan ii get printTree + } for + 1 { + 0 1 dhcone2.gblist length 1 sub { + /ii set + ii messagen ( : ) messagen + dhcone2.gblist ii get printTree + } for + } { } ifelse + + + (=========================================) message + %(cone.withGblist = ) messagen cone.withGblist message + dhcone.explain message + ( ) message + ] pop + popVariables +} def + +% +% 試し方 test14, 22, 25 +% +% (data/test14.sm1) run (data/test14-data.sm1) run +% printGrobnerFan ; % H 付きで印刷. +% dhCones_h ; % dehomogenize Cones. +% dhcone.rtable ; % dhcone2.fan 等を生成. +% dhcone.printGrobnerFan ; % 印刷. +% 印刷したものは test*-print.txt へ格納してある. +% + +% Todo: save functions. + +%< +% Collart, Kalkbrener, Mall のアルゴリズムによる gb の flip. +% See also Sturmfels' book, p.22, 23. +% Usages: [reducedGb, vlist, oldWeight, facetWeight, newWeight] ckmFlip rGb +% If it fails, then it returns null, else it returns the reducedGb for the +% newWeight. +% gb の check をやるので, それに失敗したら null を戻す. +% weight はすべて vw 形式で. vw 形式 = variable weight の繰り返しの形式 +% reducedGb は文字列のリストではなく多項式の形式のこと. +% 理由は reducedGb より ring の構造を読むため. +%> +/ckmFlip { + /arg1 set + [/arg_ckmFlip /gOld /vlist /oldWeight /facetWeight /newWeight + /gNew + /ww /ww1 /ww2 % 本の中の w1, w, w2 (古い, facet, 新しい) + /ch1 /ch2 % 本の中の {\cal H}_1, {\cal H}_2 + /grData /rTable + /rTable2 % rTable の反対の変換. + /facetWeight_gr /vlist_gr % graded ring 用. + /oldWeight_gr + /ccf % reduction した係数. + /rwork /ccf2 /gNew + ] pushVariables + [ + arg1 /arg_ckmFlip set + arg_ckmFlip 0 get /gOld set + arg_ckmFlip 1 get /vlist set + arg_ckmFlip 2 get /oldWeight set + arg_ckmFlip 3 get /facetWeight set + arg_ckmFlip 4 get /newWeight set + +% facet weight vector ww についての initial を取り出す. ch1 へいれる. + gOld getRing ring_def + facetWeight weightv /ww set + gOld { ww init } map /ch1 set % facetWeight による initial の取り出し. + + +% 例: [(x,y) [(x) -1 (Dx) 1 (y) -1 (Dy) 2]] getGrRing +% [$x,y,y',$ , [ $x$ , $y$ ] , [ [ $Dy$ , $y'$ ] ] ] +% 変数リスト 置換表 +% ch1 を gr_ww の元に変換. + [vlist facetWeight] getGrRing /grData set + [grData 0 get ring_of_differential_operators 0] define_ring /rwork set + grData 2 get { { . } map } map /rTable set + rTable { reverse } map /rTable2 set + grData 0 get /vlist_gr set + ch1 { toString . rTable replace toString } map /ch1 set + + oldWeight { dup isString { . rTable replace toString } + { } ifelse } map /oldWeight_gr set + +% facetWeight も 新しい環 gr_ww の weight に変換. +% 例. [(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 + %% ttt cone.gb_gr /ch1 set %再度の計算は不要. + [[(1)] vlist_gr oldWeight_gr] cone.gb_gr getRing ring_def % Set Ring. + ch1 {toString .} map /ch1 set +%% ここまででとりあえずテストをしよう. +%% ch1 /arg1 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 で計算. +% どちらをとるかは cone.reduction_gr で区別するしかなし + ch1 getRing ring_def ; + ch2 {toString .} map {ch1 cone.reduction} map /ccf set + %ccf pmat + % とりあえずテスト. + % [ch1 ch2] /arg1 set + %% ccf[i][0] は 0 でないと矛盾. check まだしてない. + + %% ccf[i][2] (syzygy) を gr から もとの ring へ戻し, + %% 新しい reduced gbasis を ccf[i][2] * gOld で作る. + rwork ring_def + ccf { 2 get {toString . rTable2 replace toString} map } map /ccf2 set + %% ccf2 は gr でない ring の元. + gOld getRing ring_def + cone.DhH { cone.begin_DhH } { } ifelse % Hh か h^2 か. + ccf2 { {.} map gOld mul } map /gNew set + gNew { toString } map /gNew set + cone.DhH { cone.end_DhH } { } ifelse % Hh か h^2 か. + % gNew /arg1 set + %gNew が newWeight での GB か check. Yes なら reduced basis へ. + %No なら null を戻す. +%%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 + ] pop + popVariables + arg1 +} 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.begin_DhH + ff ggbasis reduction /ans set + 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 + /cone.comment [ + (BS for y and y-(x-1)^2, t1, t2 space, in doubly homogenized Weyl algebra.) nl + (The Grobner cones are dehomogenized to get local Grobner fan.) nl + ] cat def + /cone.vlist [(t1) (t2) (x) (y) (Dt1) (Dt2) (Dx) (Dy) (h) (H)] def + /cone.vv (t1,t2,x,y) def + /cone.type 1 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 + /cone.input + [ + (t1-y) (t2 - (y-(x-1)^2)) + ((-2 x + 2)*Dt2+Dx) + (Dt1+Dt2+Dy) + ] + def + % homogenize + [cone.vv ring_of_differential_operators + [[(t1) -1 (t2) -1 (Dt1) 1 (Dt2) 1]] ecart.weight_vector + 0] define_ring + dh.begin + cone.input { . homogenize toString } map /cone.input set + dh.end + + +% テストを開始する. +% getStartingCone /cone.ncone set +% cone.ncone updateFan +% cone.gblist 0 get message +% cone.ncone /cone.ccone set +% getNextFlip /cone.nextflip set +% cone.nextflip message + + /wOld [(t1) , -29 , (t2) , -38 , (Dt1) , 29 , (Dt2) , 38 ] def + /wFacet [(t1) , -1 , (t2) , -1 , (Dt1) , 1 , (Dt2) , 1 ] def + /wNew [(t1) , -39 , (t2) , -38 , (Dt1) , 39 , (Dt2) , 38 ] def + cone.input wOld cone.gb /ff set + [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 + +%%change +/cone_ir_input { + /arg1 set + [/msg ] pushVariables + [ + /msg arg1 def + (---------------) message + msg message + ( ) message + (Please also refer to the value of the variables cone.getConeInfo.rr0) message + ( cone.getConeInfo.rr1 cone.Lp cone.cinit) message + $ cone.cinit (FACETS) getNode :: $ message + (We are sorry that we cannot accept this input.) error + ] pop + popVariables +} def