% $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.16 2009/09/04 02:59:55 takayama Exp $ % cp cone.sm1 $OpenXM_HOME/src/kan96xx/Doc/gfan.sm1 % $Id: cone.sm1,v 1.81 2005/07/07 07:53:27 taka Exp $ % iso-2022-jp %%Ref: @s/2004/08/21-note.pdf %% gfan.sm1 works only for polymake 2.0 Use webservice of 2.0. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 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 /@@@polymake.web 1 def %% Choose it automatically. [(which) (polymake)] oxshell tag 0 eq @@@polymake.web 1 eq or { (Polymake is not installed in this system or @@@polymake.web is set.) 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 { } { [(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 も見るか? % 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 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 % cone.h0: 1 なら h の weight 0 での Grobner fan を計算する. /cone.h0 1 def % Global: cone.n (number of variables in GB) % cone.m (freedom of the weight space. cf. cone.W) % 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 の満たす % 不等式制約を求める. %> /coneEq1 { /arg1 set [/g /eqs /gsize /i /j /n /f /exps /m % Do not use "eq" as a variable /expsTop ] pushVariables [ /g arg1 def % Reduced Grobner basis /eqs [ ] def % 不等式系の係数 /gsize g length def 0 1 gsize 1 sub { /i set g i get /f set % f は i 番目の reduced Grobner basis の元 [(exponents) f cone.type] gbext /exps set % exps は f の exponent vector exps length /m set m 1 eq not { /expsTop exps 0 get def % expsTop は f の先頭の exponent vector. 1 1 exps length 1 sub { /j set eqs [expsTop exps j get sub] join /eqs set % exps[0]-exps[j] を eqs へ格納していくだけ. % Cone の closure をだすので >= で OK. } for } { } ifelse } for /arg1 eqs def ] pop popVariables arg1 } def %< % Usage: ww g coneEq % ww は [v1 w1 v2 w2 ... ] 形式. (v-w 形式) w1, w2 は univNumber でもいい. % g は reduced Grobner basis % in(f) が monomial でない場合も扱う. % in_w(f) = in_ww(f) となる weight w の満たす % 不等式制約を求める. % ord_w, init (weightv) を用いる. %> /coneEq { /arg2 set /arg1 set [/g /eqs /gsize /i /j /n /f /exps /m /expsTop /ww /ww2 /iterms ] pushVariables [ /g arg2 def % Reduced Grobner basis /ww arg1 def % weight vector. v-w 形式 ww to_int32 /ww set % univNum があれば int32 に直しておく. /ww2 ww weightv def % v-w 形式を 数字のベクトルに. (init 用) /eqs null def % 不等式系の係数 /gsize g length def 0 1 gsize 1 sub { /i set g i get /f set % f は i 番目の reduced Grobner basis の元 [(exponents) f cone.type] gbext /exps set % exps は f の exponent vector exps length /m set m 1 eq not { /expsTop exps 0 get def % expsTop は f の先頭の exponent vector. /iterms f ww2 init length def % f の initial term の項の数. % in_ww(f) > f_j となる項の処理. iterms 1 exps length 1 sub { /j 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 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 arg1 } def %< % Usage: wv g coneEq genPo % polymake 形式の INEQUALITIES を生成する. coneEq -> genPo と利用 %> /genPo { /arg1 set [/outConeEq /rr /nn /ii /mm /jj /ee] pushVariables [ /outConeEq arg1 def /rr [(INEQUALITIES) nl] cat def % 文字列 rr に足していく. outConeEq length /nn set 0 1 nn 1 sub { /ii set outConeEq ii get /ee set [ rr (0 ) % 非せいじ用の 0 を加える. 0 1 ee length 1 sub { /jj set ee jj get toString ( ) } for nl ] cat /rr set } for /arg1 rr def ] pop popVariables arg1 } def %< % Usage: wv g coneEq genPo2 % doPolyamke 形式の INEQUALITIES を生成する. coneEq -> genPo2 と利用 % tfb 形式文字列. %> /genPo2 { /arg1 set [/outConeEq /rr /nn /ii /mm /jj /ee] pushVariables [ /outConeEq arg1 def /rr $polymake.data(polymake.INEQUALITIES([$ def % 文字列 rr に足していく. outConeEq length /nn set 0 1 nn 1 sub { /ii set outConeEq ii get /ee set [ rr ([0,) % 非せいじ用の 0 を加える. 0 1 ee length 1 sub { /jj set ee jj get toString jj ee length 1 sub eq { } { (,) } ifelse } for (]) ii nn 1 sub eq { } { (,) } ifelse ] cat /rr set } for [rr $]))$ ] cat /rr set /arg1 rr def ] pop popVariables arg1 } def /test1 { [(x,y) ring_of_differential_operators 0] define_ring [ (x + y + Dx + Dy). (x ^2 Dx^2 + y^2 Dy^2). (x). ] /gg set gg coneEq1 /ggc set gg message ggc pmat ggc genPo message } def /test2 { [(parse) (dhecart.sm1) pushfile] extension dh.test.p1 /ff set ff 0 get coneEq1 /ggc set ggc message ggc genPo /ss set ss message (Data is in ss) message } def /test3 { % [(parse) (cohom.sm1) pushfile] extension /ww [(Dx) 1 (Dy) 1] def [(x,y) ring_of_differential_operators [ww] weight_vector 0] define_ring [ (x Dx + y Dy -1). (y^2 Dy^2 + 2 + y Dy ). ] /gg set gg {homogenize} map /gg set [gg] groebner 0 get /gg set ww message ww gg coneEq /ggc set gg message ggc pmat ggc genPo message } def %< % Usage: test3b % Grobner cone を決定して, polymake 用のデータを生成するテスト. % weight (0,0,1,1) だと max dim cone でない. %> /test3b { % [(parse) (cohom.sm1) pushfile] extension /ww [(Dx) 1 (Dy) 2] def [(x,y) ring_of_differential_operators [ww] weight_vector 0] define_ring [ (x Dx + y Dy -1). (y^2 Dy^2 + 2 + y Dy ). ] /gg set gg {homogenize} map /gg set [gg] groebner 0 get /gg set ww message ww gg coneEq /ggc set gg message ggc pmat % ggc genPo /ggs set % INEQ を文字列形式で % ggs message % ggs output % (mv sm1out.txt test3b.poly) system % (Type in polymake-pear.sh test3b.poly FACETS) message ggc genPo2 /ggs set % INEQ を文字列形式 for doPolymake ggs message } def % commit (dr.sm1): lcm, denominator, ngcd, to_univNum, numerator, reduce % 8/22, changelog-ja まだ. % to do : nnormalize_vec, sort_vec --> shell で OK. % 8/27, getNode /test4 { $polymake.data(polymake.INEQUALITIES([[0,1,0,0],[0,0,1,0]]))$ /ff set [(FACETS) ff] doPolymake /rr set rr 1 get /rr1 set rr1 getLinearitySubspace pmat } def %< % Usage: vv ineq isInLinearSpace % vv が ineq[i] > 0 で定義される半空間のどれかにはいっているなら 0 % vv が 全ての i について ineq[i] = 0 にはいっていたら 1. %> /isInLinearSpace { /arg2 set /arg1 set [/vv /ineq /ii /rr] pushVariables [ /vv arg1 def /ineq arg2 def /rr 1 def { 0 1 ineq length 1 sub { /ii set % vv . ineq[ii] != 0 なら vv は linearity space の元でない. vv ineq ii get mul to_univNum isZero { } { /rr 0 def exit} ifelse } for exit } loop /arg1 rr def ] pop popVariables arg1 } def %< % Usages: doPolymakeObj getLinearitySubspace % INEQUALITIES と VERTICES から maximal linearity subspace % の生成ベクトルを求める. % 例: VERTICES [[0,1,0,0],[0,0,1,0],[0,0,0,-1],[0,0,0,1]]] % 例: INEQUALITIES [[0,1,0,0],[0,0,1,0]] % 入力は polymake の tree (doPolymake の 1 get) %> /getLinearitySubspace { /arg1 set [/pdata /vv /ineq /rr /ii] pushVariables [ /pdata arg1 def { /rr [ ] def % POINTED なら max lin subspace は 0. pdata (POINTED) getNode tag 0 eq { } { exit} ifelse pdata (INEQUALITIES) getNode 2 get 0 get /ineq set pdata (VERTICES) getNode 2 get 0 get /vv set 0 1 vv length 1 sub { /ii set % -vv[ii] が ineq を満たすか調べる. vv ii get ineq isInLinearSpace { rr [vv ii get] join /rr set } { } ifelse } for exit } loop /arg1 rr def ] pop popVariables arg1 } def %< % Usages: mm asir_matrix_image % 生成元より線形空間の基底を得る. %> /asir_matrix_image { /arg1 set [/mm /rr] pushVariables [(CurrentRingp)] pushEnv [ /mm arg1 def mm to_univNum /mm set oxasir.ccc [ ] eq { (Starting ox_asir server.) message ox_asirConnectMethod } { } ifelse { oxasir.ccc [(matrix_image) mm] asir /rr set rr null_to_zero /rr set exit (asir_matrix_image: not implemented) error exit } loop rr numerator /rr set /arg1 rr def ] pop popEnv popVariables arg1 } def [(asir_matrix_image) [(Calling the function matrix_image of asir. It gets a reduced basis of a given matrix.) (Example: [[1 2 3] [2 4 6]] asir_matrix_image) ]] putUsages %< % Usages: mm asir_matrix_kernel % 直交する空間の基底. %> /asir_matrix_kernel { /arg1 set [/mm /rr] pushVariables [(CurrentRingp)] pushEnv [ /mm arg1 def mm to_univNum /mm set oxasir.ccc [ ] eq { (Starting ox_asir server.) message ox_asirConnectMethod } { } ifelse { oxasir.ccc [(matrix_kernel) mm] asir /rr set rr null_to_zero /rr set exit (asir_matrix_image: not implemented) error exit } loop rr 1 get numerator /rr set /arg1 rr def ] pop popEnv popVariables arg1 } def [(asir_matrix_kernel) [(Calling the function matrix_kernel of asir.) (It gets a reduced basis of the kernel of a given matrix.) (Example: [[1 2 3] [2 4 6]] asir_matrix_kernel) ]] putUsages %< % Usages: v null_to_zero %> /null_to_zero { /arg1 set [/pp /rr] pushVariables [ /pp arg1 def { /rr pp def pp isArray { pp {null_to_zero} map /rr set exit }{ } ifelse pp tag 0 eq { /rr (0).. def exit }{ } ifelse exit } loop /arg1 rr def ] pop popVariables arg1 } def [(null_to_zero) [(obj null_to_zero rob) $It translates null to (0)..$ ]] putUsages %< % Usages: newVector.with-1 % (-1).. で埋めたベクトルを作る. %> /newVector.with-1 { newVector { pop (-1).. } map } def % [2 0] lcm は 0 をもどすがいいか? --> OK. %< % Usages: mm addZeroForPolymake % 以下の二つの関数は, toQuotientSpace にも利用. % Polymake INEQUALITIES 用に 0 を始めに足す. % 入力は リストのリスト % [[1,2], [3,4],[5,6]] --> [[0,1,2],[0,3,4],[0,5,6]] %> /addZeroForPolymake { /arg1 set [/mm /rr] pushVariables [ /mm arg1 def mm to_univNum /mm set mm { [(0)..] 2 1 roll join } map /mm set /arg1 mm def ] pop popVariables arg1 } def %< % Usages: mm cone.appendZero %> /cone.appendZero { /arg1 set [/mm /rr] pushVariables [ /mm arg1 def mm to_univNum /mm set mm { [(0)..] join } map /mm set /arg1 mm def ] pop popVariables arg1 } def %< % Usages: mm removeFirstFromPolymake % 始めの 0 を取り除く. % 入力は リストのリスト % [[0,1,2],[0,3,4],[0,5,6]] ---> [[1,2], [3,4],[5,6]] %> /removeFirstFromPolymake { /arg1 set [/mm /rr] pushVariables [ /mm arg1 def mm to_univNum /mm set mm {rest} map /mm set /arg1 mm def ] pop popVariables arg1 } def %< % Usages: mm genUnit % [1,0,0,...] を加えるために生成. % [[0,1,2], [0,3,4],[0,5,6]]--> [1,0,0] %> /genUnit { /arg1 set [/mm /rr /i] pushVariables [ /mm arg1 def mm 0 get length newVector /rr set rr null_to_zero /rr set rr 0 (1).. put /arg1 rr def ] pop popVariables arg1 } def %< % Usages: mm genUnitMatrix % [[0,1,2], [0,3,4],[0,5,6]]--> [[1,0,0],[0,1,0],[0,0,1]] %> /genUnitMatrix { /arg1 set [/mm /rr /nn /i] pushVariables [ /mm arg1 def mm 0 get length /nn set [ 0 1 nn 1 sub { /i set nn newVector null_to_zero /mm set mm i (1).. put mm } for ] /arg1 set ] pop popVariables arg1 } def %< %%note: 2004, 8/29 (sun) % toQuotientSpace : Linearity space で割る. % Usages: ineq mm toQuotientSpace % 入力は coneEq の出力 ineq % および doPolymake --> getLinearitySubspace ==> L % [L,[1,0,0,...]] asir_matrix_kernel removeFirstFromPolymake で得られた mm % 出力から 0 ベクトルは削除. % 出力も coneEq 形式. 特に polymake 用に 0 を加えるのが必要. % ref: getUnit, removeFirstFromPolymake, addZeroForPolymake, % asir_matrix_kernel, getLinearitySubspace %> /toQuotientSpace { /arg2 set /arg1 set [/ineq /mm /rr] pushVariables [ /ineq arg1 def /mm arg2 def ineq mm transpose mul /rr set /arg1 rr def ] pop popVariables arg1 } def /test5.data $polymake.data(polymake.INEQUALITIES([[0,1,-1,1,-1,0],[0,0,-1,0,-1,2],[0,0,-1,0,-1,2],[0,0,-2,0,-2,4],[0,-1,0,-1,0,2],[0,-2,0,-2,0,4]]),polymake.VERTICES([[0,0,-1,0,0,0],[0,-1,-1,0,0,0],[0,1,0,-1,0,0],[0,-1,0,1,0,0],[0,0,1,0,-1,0],[0,0,-1,0,1,0],[0,-2,-2,0,0,-1],[0,2,2,0,0,1]]),polymake.FACETS([[0,1,-1,1,-1,0],[0,-1,0,-1,0,2]]),polymake.AFFINE_HULL(),polymake.FEASIBLE(),polymake.NOT__POINTED(),polymake.FAR_FACE([polymake._set([0,1,2,3,4,5,6,7])]),polymake.VERTICES_IN_INEQUALITIES([polymake._set([1,2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7])]),polymake.DIM([[5]]),polymake.AMBIENT_DIM([[5]]))$ def %< % Usages: test5 %% getConeInfo を変更すれば polymake を呼ばずにテストできる. %> /test5 { % test3b より /ww [(Dx) 1 (Dy) 2] def % /ww [(x) 1 (y) -2 (Dx) 3 (Dy) 6] def [(x,y) ring_of_differential_operators [ww] weight_vector 0] define_ring [ (x Dx + y Dy -1). (y^2 Dy^2 + 2 + y Dy ). ] /gg set gg {homogenize} map /gg set [(AutoReduce) 1] system_variable [gg] groebner 0 get /gg set ww message ww gg coneEq getConeInfo /rr set (Type in rr 0 get :: ) message } def %[5, [[1,0,1,0,-2],[0,1,0,1,-2]], $NOT__POINTED$ ] % この場合は 2 次元まで落すと pointed cone になる. % coneEq mmc transpose をもとに FACETS を計算すればよい. %< % Usage: ceq getConeInfo % vw は [v1 w1 v2 w2 ... ] 形式. (v-w 形式) w1, w2 は univNumber でもいい. % g は reduced Grobner basis として vw g coneEq を計算. これを getConeInfo へ. % Grobner cone の 次元 cdim (DIM), 補空間 (linearity space ) への行列 mmc % linearity space 自体, pointed or not__pointed % つまり [cdim, L', L, PointedQ] % を計算して戻す. (polymake 形式の余分な部分なし) % polymake 必要. % ref: coneEq % Global: % cone.getConeInfo.rr0, cone.getConeInfo.rr1 に polymake よりの戻り値がはいる. %> /getConeInfo { /arg1 set [/ww /g /ceq /ceq2 /cdim /mmc /mmL /rr /ineq /ppt /rr0 /mm0 /mm1] pushVariables [ /ceq arg1 def ceq pruneZeroVector /ceq set ceq length 0 eq { (Monomial ideal is not accepted as an input.) cone_ir_input } { } ifelse /mm1 ( Use [(keep_tmp_files) 1] oxshell to check the input to polymake2tfb. See /tmp or $TMP ) def ceq genPo2 /ceq2 set % ceq2 は polymake.data(polymake.INEQUALITIES(...)) 形式 % polymake で ceq2 の次元の計算. /getConeInfo.ceq ceq def /getConeInfo.ceq2 ceq2 def cone.debug { (Calling polymake DIM.) message } { } ifelse [(DIM) ceq2] doPolymake /rr0 set % rr0 2 get message rr0 2 get 1 get 0 get /mm0 set mm0 length 0 eq { } { [mm0 mm1] cat error } ifelse rr0 1 get /rr set cone.debug {(Done.) message } { } ifelse % test5 には次のコメントとりさる. 上の行をコメントアウト. % test5.data tfbToTree /rr set /cone.getConeInfo.rr0 rr def rr (DIM) getNode /cdim set cdim 2 get 0 get 0 get 0 get to_univNum /cdim set % polymake の DIM は一つ小さいので 1 足す. cdim (1).. add /cdim set rr (FACETS) getNode tag 0 eq { % FACETS を持っていないなら再度計算する. % POINTED, NOT__POINTED も得られる cone.debug { (Calling polymake FACETS.) message } { } ifelse [(FACETS) ceq2] doPolymake /rr0 set % rr0 2 get message rr0 2 get 1 get 0 get /mm0 set mm0 length 0 eq { } { [mm0 mm1] cat error } ifelse rr0 1 get /rr set cone.debug { (Done.) message } { } ifelse } { } ifelse rr (VERTICES) getNode tag 0 eq { (internal error: VERTICES is not found.) error } { rr (VERTICES) getNode (UNDEF) getNode tag 0 eq { } { (internal error: VERTICES is UNDEF. See rr. Set /@@@polymake.web 1 def) error } ifelse } ifelse /cone.getConeInfo.rr1 rr def rr (NOT__POINTED) getNode tag 0 eq { % cone が pointed の時は mmc は単位行列. genUnitMatrix を使う. % VERTICES より一つ小さいサイズ. /mmc [ rr (VERTICES) getNode 2 get 0 get 0 get rest] genUnitMatrix def /mmL [ ] def /ppt (POINTED) def } { % pointed でない場合, % cone の線形部分空間を計算. rr getLinearitySubspace /mmL set [mmL genUnit] mmL join /mmc set % [1,0,0,...] を足す. mmc asir_matrix_kernel /mmc set % 補空間 mmc removeFirstFromPolymake /mmc set % ひとつ小さいサイズに. [mmL genUnit] mmL join asir_matrix_image removeFirstFromPolymake /mmL set mmL asir_matrix_image /mmL set % Linearity space を求める. rm 0vector /ppt (NOT__POINTED) def } ifelse /arg1 [[cdim mmc mmL ppt] rr] def ] pop popVariables arg1 } def /test.put { /dog [(dog) [[(legs) 4] ] [1 2 3 ]] [(class) (tree)] dc def /man [(man) [[(legs) 2] ] [1 2 3 ]] [(class) (tree)] dc def /ma [(mammal) [ ] [man dog]] [(class) (tree)] dc def /fan [ma 1 copy] def ma (dog) getNode /dd set dd 2 get /dd2 set dd2 1 0 put ma message fan message } def /test6.data $polymake.data(polymake.INEQUALITIES([[0,1,-1,1,-1,0],[0,0,-1,0,-1,2],[0,0,-1,0,-1,2],[0,0,-2,0,-2,4],[0,-1,0,-1,0,2],[0,-2,0,-2,0,4]]),polymake.VERTICES([[0,0,-1,0,0,0],[0,-1,-1,0,0,0],[0,1,0,-1,0,0],[0,-1,0,1,0,0],[0,0,1,0,-1,0],[0,0,-1,0,1,0],[0,-2,-2,0,0,-1],[0,2,2,0,0,1]]),polymake.FACETS([[0,1,-1,1,-1,0],[0,-1,0,-1,0,2]]),polymake.AFFINE_HULL(),polymake.FEASIBLE(),polymake.NOT__POINTED(),polymake.FAR_FACE([polymake._set([0,1,2,3,4,5,6,7])]),polymake.VERTICES_IN_INEQUALITIES([polymake._set([1,2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7])]))$ def % tfbToTree /arrayToTree { [(class) (tree)] dc } def %< % polymake より得られた TreeObject から TreeObject cone を生成する. % Usages: test6.data tfbToTree newCone で動作テスト %> /test6 { test6.data tfbToTree /rr set rr newCone /rr2 set } def %< % Usages: doPolymakeObj newCone %> /newCone { /arg1 set [/polydata /cone /facets /vertices /flipped /ineq /facetsv /rr] pushVariables [ /polydata arg1 def polydata (FACETS) getNode tag 0 eq { (newCone : no FACETS data.) error } { } ifelse % facets は有理数の場合正規化する. data/test11 で 有理数でる. 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 [[ ] ] vertices join shell rest removeFirstFromPolymake /vertices set % inequalities は有理数の場合正規化する. polydata (INEQUALITIES) getNode 2 get 0 get to_univNum { 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 ] ] arrayToTree /cone set /arg1 cone def ] pop popVariables arg1 } def %< % Usages: newCone_facetv % facet vertices newCone_facetv % facet にのっている vertices をすべて列挙. %> /newCone_facetv { /arg2 set /arg1 set [/facet /vertices] pushVariables [ /facet arg1 def /vertices arg2 def [ 0 1 vertices length 1 sub { /ii set facet vertices ii get mul isZero { vertices ii get } { } ifelse } for ] /arg1 set ] pop popVariables arg1 } def %< % Usages: newCone_facetsv % facets vertices newCone_facetv % facets にのっている vertices をすべて列挙. リストを作る. %> /newCone_facetsv { /arg2 set /arg1 set [/facets /vertices] pushVariables [ /facets arg1 def /vertices arg2 def facets { vertices newCone_facetv } map /arg1 set ] pop popVariables arg1 } 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 /cone_random { [(tdiv_qr) cone_random.start (1103515245).. mul (12345).. add (2147483646).. ] mpzext 1 get /cone_random.start set cone_random.start } def /cone_random.limit 40 def /cone_random_vec { /arg1 set [/nn /rr] pushVariables [ /nn arg1 def [ 0 1 nn 1 sub { pop [(tdiv_qr) cone_random cone_random.limit] mpzext 1 get } for ] /arg1 set ] pop popVariables arg1 } def %< % Usages: getNewRandomWeight %% max dim の cone を生成するために, random な weight を生成する. %% h, H の処理も必要. %% 制約条件 u+v >= 2t をみたす weight が必要. これをどのように作るのか? %> /getNewRandomWeight { /arg1 set [/vv /vvd /rr] pushVariables [ /vv arg1 def vv { (D) 2 1 roll 2 cat_n } map /vvd set ] pop popVariables arg1 } def % test7 : univNum の weight が正しく認識されるかのテスト % aux-cone.sm1 %< % Usages: n d coneEqForSmallFan.2 (cone.type 2 専用: x,y,Dx,Dy,h) % n 変数の数, d zero にしない変数の数. d は max dim cone の次元となる. % はじめから d 個の変数. % 4, 2 , s,t,x,y なら weight は s,t,Ds,Dt のみ. % u_i + v_i >= 0 , u_i = v_i = 0. % homog 変数の条件 u_i+v_i >= t, i.e, -t >= 0 も入れる. % coneEq の結果と coneEqForSmallFan.2 の結果を join して % getConeInfo or newCone % note-cone.sm1 2004.8.31 を見よ. w_ineq あたり. % cone.local が設定されていると u_i <= 0 も条件に入る. %> /coneEqForSmallFan.2 { /arg2 set /arg1 set [/n /d /nn /dd /ii /tt] pushVariables [ /n arg1 def /d arg2 def n to_int32 /n set d to_int32 /d set /dd n d add def /nn n n add def % 0 ~ d-1, n ~ dd-1 では u_i + v_i = 0 % d ~ n-1, dd ~ nn-1 では u_i=v+i = 0. % -t >= 0 [ % d ~ n-1, dd ~ nn-1 では u_i=v+i = 0. d 1 n 1 sub { /ii set % [ 0,0, ..., 0,1,0,... ; 0] を生成 nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt % [ 0,0, ..., 0,-1,0,... ; 0] を生成 nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for dd 1 nn 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for % 0 ~ d-1, n ~ dd-1 では u_i + v_i = 0 0 1 d 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt ii n add (1).. put tt nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt ii n add (-1).. put tt } for % -t >= 0 cone.h0 { % t = 0 nn 1 add newVector null_to_zero /tt set tt nn (1).. put tt nn 1 add newVector null_to_zero /tt set tt nn (-1).. put tt } { % -t >= 0 nn 1 add newVector null_to_zero /tt set tt nn (-1).. put tt } ifelse % cone.local が 1 の時 % 0 ~ d-1 では -u_i >= 0 cone.local { 0 1 d 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for } { } ifelse ] /rr set /arg1 rr to_univNum def ] pop popVariables arg1 } def %< % Usages: n d coneEqForSmallFan.1 (cone.type 1 専用: x,y,Dx,Dy,h,H) % cone.type 2 では x,y,Dx,Dy,h % coneEqForSmallFan.2 の結果を用いて生成. % H の条件を加える. %> /coneEqForSmallFan.1 { /arg2 set /arg1 set [/n /d /i /j /rr /tt /tt2] pushVariables [ /n arg1 def /d arg2 def n d coneEqForSmallFan.2 /rr set rr cone.appendZero /rr set % H 用の 0 を加える. % とりあえず t' = 0 できめうち. cone.h0 { } { (cone.h0 = 0 has not yet been implemented.) error } ifelse n 2 mul 2 add newVector null_to_zero /tt set tt n 2 mul 2 add 1 sub (-1).. put n 2 mul 2 add newVector null_to_zero /tt2 set tt2 n 2 mul 2 add 1 sub (1).. put rr [tt tt2] join /rr set /arg1 rr to_univNum def ] pop popVariables arg1 } def %< % Usages: vv ineq toQuotientCone % weight space の パラメータつけのために使う. % cone.V を求めたい. vv は doPolymakeObj (VERTICES) getNode 2 get 0 get で得る. % vertices の non-negative combination が cone. % vertice cone.w_ineq isInLinearSubspace なら取り除く. % つまり vertice*cone.w_ineq = 0 なら取り除く. % % これで正しい? 証明は? まだ途中. cone.W を求めるのに使う. (BUG) % cone.w_cone 1 get (VERTICES) getNode :: と比較せよ. % この関数を呼んで cone.W を作るのは不要かも. % % Example: cf. parametrizeSmallFan % 4 2 coneEqForSmallFan.2 /cone.w_ineq set cone.w_ineq getConeInfo /rr set % rr 1 get (VERTICES) getNode 2 get 0 get removeFirstFromPolymake /vv set % vv cone.w_ineq toQuotientCone pmat %> /toQuotientCone { /arg2 set /arg1 set [/vv /ineq /rr] pushVariables [ /vv arg1 def /ineq arg2 def vv { dup ineq isInLinearSpace 1 eq { pop } { } ifelse } map /arg1 set ] pop popVariables arg1 } def %< % Usages: n d parametrizeSmallFan % n : x 変数の数. % d : 0 にしない weight の数. % 次の大域変数も設定される. % cone.W : weight をパラメータづけするベクトルの組. % cone.Wpos : i が 0 ~ Wpos-1 の範囲のとき V[i] へは N の元を掛け算してよい, % i が Wpos ~ の範囲のとき V[i] へは Z の元を掛け算してよい. % cone.w_ineq : weight space の不等式制約. 以後の計算で常に付加する. % cone.w_cone : w_ineq を polymake で getConeInfo した結果. % Example: /cone.local 1 def ; 4 2 parametrizeSmallFan pmat % Example: /cone.local 0 def ; 4 2 parametrizeSmallFan pmat %> /parametrizeSmallFan { /arg2 set /arg1 set [/n /d /vv /coneray] pushVariables [ /n arg1 def /d arg2 def { cone.type 1 eq { n d coneEqForSmallFan.1 /cone.w_ineq set exit } { } ifelse cone.type 2 eq { n d coneEqForSmallFan.2 /cone.w_ineq set exit } { } ifelse (This cone.type has not yet been implemented.) error } loop cone.w_ineq getConeInfo /cone.w_cone set cone.w_cone 1 get (VERTICES) getNode 2 get 0 get removeFirstFromPolymake /vv set vv cone.w_ineq toQuotientCone /coneray set coneray length /cone.Wpos set coneray cone.w_cone 0 get 2 get join /cone.W set /arg1 cone.W def ] pop popVariables arg1 } def %< % Usages: n d coneEqForTotalFan.2 (cone.type 2 専用: x,y,Dx,Dy,h) % n 変数の数, % d 0 にしない変数. % u_i + v_i >= 0 , % homog 変数の条件 u_i+v_i >= 0, t = 0 も入れる. % coneEq の結果と coneEqForSmallFan.2 の結果を join して % getConeInfo or newCone % cone.local が設定されていると u_i <= 0 も条件に入る. %> /coneEqForTotalFan.2 { /arg2 set /arg1 set [/n /nn /dd /ii /tt] pushVariables [ /n arg1 def /d arg2 def n to_int32 /n set d to_int32 /d set /nn n n add def /dd n d add def % 0 ~ d-1, n ~ dd-1 では u_i + v_i >= 0 % d ~ n-1, dd ~ nn-1 では u_i=v+i = 0. % t = 0 [ % d ~ n-1, dd ~ nn-1 では u_i=v+i = 0. d 1 n 1 sub { /ii set % [ 0,0, ..., 0,1,0,... ; 0] を生成 nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt % [ 0,0, ..., 0,-1,0,... ; 0] を生成 nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for dd 1 nn 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for % 0 ~ d-1, n ~ dd-1 では u_i + v_i >= 0 0 1 d 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (1).. put tt ii n add (1).. put tt } for % t = 0 cone.h0 { % t = 0 nn 1 add newVector null_to_zero /tt set tt nn (1).. put tt nn 1 add newVector null_to_zero /tt set tt nn (-1).. put tt } { (coneForTotalFan.2. Not implemented.) error } ifelse % cone.local が 1 の時 % 0 ~ d-1 では -u_i >= 0 cone.local { 0 1 d 1 sub { /ii set nn 1 add newVector null_to_zero /tt set tt ii (-1).. put tt } for } { } ifelse ] /rr set /arg1 rr to_univNum def ] pop popVariables arg1 } def %< % Usages: n d parametrizeTotalFan % n : x 変数の数. % d : 0 にしない数. % 次の大域変数も設定される. % cone.W : weight をパラメータづけするベクトルの組. % cone.Wpos : i が 0 ~ Wpos-1 の範囲のとき V[i] へは N の元を掛け算してよい, % i が Wpos ~ の範囲のとき V[i] へは Z の元を掛け算してよい. % cone.w_ineq : weight space の不等式制約. 以後の計算で常に付加する. % cone.w_ineq を getConeInfo した結果は cone.w_cone % Example: /cone.local 1 def ; 3 parametrizeSmallFan pmat % Example: /cone.local 0 def ; 3 parametrizeSmallFan pmat % local が 1 だと u_i <= 0 になる. %> /parametrizeTotalFan { /arg2 set /arg1 set [/n /d /vv /coneray] pushVariables [ /n arg1 def /d arg2 def { cone.type 2 eq { n d coneEqForTotalFan.2 /cone.w_ineq set exit} { } ifelse (This cone.type has not yet been implemented.) error } loop cone.w_ineq getConeInfo /cone.w_cone set cone.w_cone 1 get (VERTICES) getNode 2 get 0 get removeFirstFromPolymake /vv set vv cone.w_ineq toQuotientCone /coneray set coneray length /cone.Wpos set coneray cone.w_cone 0 get 2 get join /cone.W set /arg1 cone.W def ] pop popVariables arg1 } def %< % Usages: vlist wlist cone_wtowv % [x y Dx Dy h] [-1 0 1 0 0] ==> [(x) -1 (Dx) 1] を作る. %> /cone_wtowv { /arg2 set /arg1 set [/vlist /wlist /ii] pushVariables [ /vlist arg1 def /wlist arg2 def wlist length vlist length eq { } { (cone_wtowv: length of the argument must be the same. Please check the values of cone.vlist cone.vv cone.type parametrizeWeightSpace) error} ifelse wlist to_int32 /wlist set [ 0 1 wlist length 1 sub { /ii set wlist ii get 0 eq { } { vlist ii get wlist ii get } ifelse } for ] /arg1 set ] pop popVariables arg1 } 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 [/mm /ii /jj /tt] pushVariables [ /mm arg1 def mm to_univNum /mm set [ [ ] ] mm join shell rest uniq /mm set [ 0 1 mm length 1 sub { /ii set mm ii get /tt set { 0 1 tt length 1 sub { /jj set tt jj get (0).. eq { } { tt exit } ifelse } for exit } loop } for ] /arg1 set ] pop arg1 } def %< % Usages: a projectIneq v , dim(a) = n, dim(v) = d % a*cone.Wt*cone.Lpt %> /projectIneq { cone.Wt mul cone.Lpt mul } def %< % Usages: v liftWeight [w vw], dim(v) = d, dim(w) = n, vw : vw 形式の weight % v*cone.Lp*cone.W cone.vlist w cone_wtowv %> /liftWeight { /arg1 set [/v /w /vw] pushVariables [ /v arg1 def v cone.Lp mul cone.W mul /w set [w cone.vlist w cone_wtowv] /arg1 set ] pop popVariables arg1 } def %< % Usage: m isZero % dr.sm1 へ移す. %> /isZero { /arg1 set [/mm /ans /ii] pushVariables [ /mm arg1 def /ans 1 def mm isArray { 0 1 mm length 1 sub { /ii set mm ii get isZero /ans set ans 0 eq { exit } { } ifelse } for } { { mm tag 1 eq {/ans mm 0 eq def exit} { } ifelse mm isPolynomial { /ans mm (0). eq def exit } { } ifelse mm isUniversalNumber { /ans mm (0).. eq def exit } { } ifelse /ans 0 def exit } loop } ifelse /arg1 ans def ] pop popVariables arg1 } def [(isZero) [(m isZero bool)]] putUsages %< % Usage: m isNonNegative % dr.sm1 へ移す. %> /isNonNegative { /arg1 set [/mm /ans /ii] pushVariables [ /mm arg1 def /ans 1 def mm isArray { 0 1 mm length 1 sub { /ii set mm ii get isNonNegative /ans set ans 0 eq { exit } { } ifelse } for } { { mm tag 1 eq {/ans mm 0 gt mm 0 eq or def exit} { } ifelse mm isUniversalNumber { /ans mm (0).. gt mm (0).. eq or def exit } { } ifelse mm isRational { mm (numerator) dc mm (denominator) dc mul /mm set /ans mm (0).. gt mm (0).. eq or def exit } { } ifelse /ans 0 def exit } loop } ifelse /arg1 ans def ] pop popVariables arg1 } def [(isNonNegative) [(m isNonNegative bool) (In case of matrix, m[i,j] >= 0 must hold for all i,j.) ]] putUsages % Global variable: cone.weightBorder % /cone.weightBorder null def 不要であろう. getStartingCone で設定される. %< % Usages: cone i isOnWeigthBorder % cone の i 番目の facet が weight 空間の境界にあるか? % 大域変数 cone.weightBorder が設定されてること. % この変数は cone の facet ベクトルのリスト. % この変数は setWeightBorder で設定 % cone.weightBorder[0] or cone.weightBorder[1] or ... % /ccone cone.startingCone def ccone 0 isOnWeightBorder % ccone 1 isOnWeightBorder %> /isOnWeightBorder { /arg2 set /arg1 set [/cone /facet_i /i /j /vv /co /ans] pushVariables [ /cone arg1 def /facet_i arg2 def facet_i to_int32 /facet_i set /ans 0 def cone (facetsv) getNode 2 get facet_i get /vv set % Facet を vertex 表現. { 0 1 cone.weightBorder length 1 sub { /i set cone.weightBorder i get /co set % co に制約条件 vv cone.Lp mul % vv を weight space へ lift. co mul isZero { /ans 1 def exit } { } ifelse } for exit } loop /arg1 ans def ] pop popVariables arg1 } def %< % Usages: cone i markFlipped % cone の i 番目の facet に flipped の印をつける. cone 自体が変更される. % cone は class-tree. Constructor は newCone %> /markFlipped { /arg2 set /arg1 set [/cone /facet_i /vv] pushVariables [ /cone arg1 def /facet_i arg2 def facet_i to_int32 /facet_i set cone (flipped) getNode 2 get /vv set vv facet_i (1).. put ] pop 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 を戻す. % それがないときは null %> /getNextFacet { /arg1 set [/cone /facet_i /vv /ii] pushVariables [ /cone arg1 def /facet_i null def cone (flipped) getNode 2 get /vv set 0 1 vv length 1 sub { /ii set vv ii get to_int32 0 eq { /facet_i ii def exit } { } ifelse } for /arg1 facet_i def ] pop popVariables arg1 } def %< % Usages: cone i epsilon flipWeight % cone の i 番目の facet にかんして flip する. % 新しい weight を求める. cf. liftWeight %> /flipWeight { /arg3 set /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 /ep arg3 def ep to_univNum (1).. div /ep set % 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 vp v ii get add /vp set } for vp ep f mul sub /vp set vp nnormalize_vec /vp set /arg1 vp def ] pop popVariables arg1 } def %< % Usages: cone1 cone2 isSameCone bool % cone1 cone2 が等しいか? facet で比べる. % cone1, cone2 は pointed cone でないといけない. %> /isSameCone { /arg2 set /arg1 set [/cone1 /cone2 /facets1 /facets2 /ans] pushVariables [ /cone1 arg1 def /cone2 arg2 def /facets1 cone1 (facets) getNode 2 get def /facets2 cone2 (facets) getNode 2 get def facets1 length facets2 length eq { facets1 facets2 sub isZero /ans set } { /ans 0 def } ifelse /arg1 ans def ] pop popVariables arg1 } def %< % Usages: cone1 cone2 getCommonFacet list % cone1 の中で cone2 に含まれる facet のリスト % cone2 の中で cone1 に含まれる facet のリストをもどす. % [1 [i] [j]] あるとき. [0 [ ] [ ]] ないとき. % cone1 の facetsv[i] が cone2 に含まれるか調べる. % cone2 の facetsv[i] が cone1 に含まれるか調べる. % cone1, cone2 は pointed cone でないといけない. %> /getCommonFacet { /arg2 set /arg1 set [/cone1 /cone2 /facets /ineq /ans1 /ans2 /i /tt] pushVariables [ /cone1 arg1 def /cone2 arg2 def /facets cone1 (facetsv) getNode 2 get def /ineq cone2 (inequalities) getNode 2 get def /ans1 [ 0 1 facets length 1 sub { /i set facets i get /tt set % facetsv[i] を tt へ. ineq tt transpose mul isNonNegative { i } { } ifelse } for ] def /facets cone2 (facetsv) getNode 2 get def /ineq cone1 (inequalities) getNode 2 get def /ans2 [ 0 1 facets length 1 sub { /i set facets i get /tt set % facetsv[i] を tt へ. ineq tt transpose mul isNonNegative { i } { } ifelse } for ] def ans1 length 1 gt ans2 length 1 gt or { (getCommonFacet found more than 1 common facets.) error } { } ifelse % 共通 facet があれば 1, なければ 0. ans1 length 1 eq ans2 length 1 eq and { /tt 1 def } { /tt 0 def } ifelse /arg1 [tt ans1 ans2] def ] pop popVariables arg1 } def % % ------------------------------------------------- % test8 は aux-cone.sm1 へ移動. % 以下いよいよ一般のプログラムの作成開始. % ------------------------------------------------- % %< % Usages: setWeightBorder % cone.weightBorder (weight cone の facet ベクトルの集合) を設定する. % あと副産物として cone.w_cone_projectedWt (doPolymakeObj) % cone.w_ineq_projectedWt % cone.m 次元のベクトル. % cone.W, cone.Wt, cone.w_ineq がすでに計算ずみでないといけない. %> /setWeightBorder { [ (Entering setWeightBorder ) message cone.w_ineq cone.Wt mul pruneZeroVector /cone.w_ineq_projectedWt set { cone.w_ineq_projectedWt length 0 eq { % weight の空間に border がない場合. /cone.weightBorder [ ] def exit } { } ifelse % weight の空間に border がある場合. cone.w_ineq_projectedWt getConeInfo /cone.w_cone_projectedWt set cone.w_cone_projectedWt 0 get 0 get to_int32 cone.m to_int32 eq { } { (setWeightBorder : internal error.) message } ifelse cone.w_cone_projectedWt 1 get (FACETS) getNode 2 get 0 get removeFirstFromPolymake /cone.weightBorder set exit } loop (cone.weightBorder=) message cone.weightBorder pmat ] pop } def % % ------------------------------------------------- % プログラムの流れ. % Global: cone.fan cone を配列として格納する. % % ncone (next cone) が新規に得られた cone であるとする. % このとき次の操作をする. % 0. ncone が cone.fan にすでにないか調べる. あれば, internal error. % 1. ncone markBorder ; ncone の中の border 上の facet を mark % 2. cone.fan の中の cone と共通 facet がないか調べ (getCommonFacet), % あればそれらを mark する. % global: cone.incidence に 共通facet を持つ組みの情報を加える. % 3. ncone を cone.fan の最後に加える. % 以上の操作をまとめたものが ncone updateFan % % getNextFlip は cone.fan の中から flip してない cone と facet の組を戻す. % なければ null を戻す. null が戻ればプログラム終了. % % getStargingCone は計算を出発すべき新規の cone を計算する. 大域変数 cone.Lt, cone.W % などもこの中で設定する. % 変数リスト, weight space を生成する関数, 入力多項式, weight の候補 等は大域変数 % として入力しておく. % % reduced gb は 関数 input weight cone.gb reduced_G で計算する. % % % [ccone i] getNextCone ncone : flip により次の cone を得る. % % 1. clearGlobals ; 入力大域変数の設定. % 2. getStartingCone /ncone set % 3. { ncone updateFan % 4. getNextFlip /cone.nextflip set % 6. cone.nextflip isNull { exit } { } ifelse % 7. cone.nextflip getNextCone /ncone set % 8. } loop % % % ------------------------------------------------- % %< % Usages: input weight cone.gb_Dh reduced_G % gb in h[1,1](D) %> /cone.gb_Dh { /arg2 set /arg1 set [/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 %(---) 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 popVariables arg1 } def %< % Usages: cone.boundp % /cone.boundp { dup boundp 2 1 roll tag 0 eq not and } def %< % Usages: clearGlobals % cf. cone.boundp % polymake を再度呼ぶために global 変数をクリアする. % まだ途中. %> /clearGlobals { /cone.W null def /cone.Wt null def /cone.cinit null def /cone.weightBorder null def } def %< % Usages: getStartingCone ncone % getStargingCone は計算を出発すべき新規の cone を計算する. % 設定すべき大域変数は以下を見よ. %> /getStartingCone.test { %------------------Globals---------------------------------------- % --------------- 入力データ用大域変数の設定 -------------------------- % % cone.input : 入力多項式系 /cone.input [(t1-x-y) (h*t2-x^2-y^2) (2*x*Dt2+h*Dt1+h*Dx) (2*y*Dt2+h*Dt1+h*Dy)] def % cone.vlist : 全変数のリスト /cone.vlist [(t1) (t2) (x) (y) (Dt1) (Dt2) (Dx) (Dy) (h)] def % cone.vv : define_ring 形式の変数リスト. % t1,t2, x,y : t-space の Grobner fan (local) を求める. /cone.vv (t1,t2,x,y) def % cone.parametrizeWeightSpace : weight 空間を parametrize する関数. % 大域変数 cone.W , cone.Wpos もきまる. /cone.parametrizeWeightSpace { 4 2 parametrizeSmallFan } def % cone.w_start : weight空間における weight の初期値. % この値で max dim cone が得られないと random weight による サーチが始まる. /cone.w_start [ 1 4 ] def % cone.gb : gb を計算する関数. /cone.gb { cone.gb_Dh } def % % ----------------- おわり --------------------------- % } def % end of getStartingCone.test /getStartingCone { [/wv_start /w_start /reduced_G] pushVariables [ % cone.n は自動的にきめられる. % cone.n は GB を計算する空間の次元. /cone.n cone.vlist length def %[1] cone.W, cone.Wpos を求める. cone.m は cone.W より自動的にきまる. % cone.m は weight 空間の自由度. cone.W で射影される先の次元. /cone.W cone.boundp { (Skip cone.parametrizeWeightSpace. cf. clearGlobals) message } { cone.parametrizeWeightSpace } ifelse (parametrizing weight space: cone.W = ) messagen cone.W message /cone.Wt cone.W transpose def /cone.m cone.W length def % WeightBorder の条件判定 facet を設定. /cone.weightBorder cone.boundp { (Skip setWeightBorder cf. clearGlobals) message } { setWeightBorder } ifelse %[2] weight vector wv_start を生成する. % wv_start を設定. cone.w_start tag 0 eq { % cone.w_start が null なら random に weight を設定. /cone.w_start cone.m cone_random_vec def } { cone.w_start length cone.m to_int32 eq { } { (Error: cone.w_start has wrong length.) error /cone.w_start cone.m cone_random_vec def } ifelse } ifelse /w_start cone.w_start cone.W mul def { cone.vlist w_start cone_wtowv /wv_start set (Trying a starting weight vector : ) messagen wv_start pmat %[3] reduced GB の計算. cone.input wv_start cone.gb /reduced_G set (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 cone.g_ineq cone.w_ineq join /cone.gw_ineq set cone.gw_ineq cone.Wt mul /cone.gw_ineq_projectedWt set % 射影 /cone.cinit cone.boundp { (Skipping cone.gw_ineq_projectedWt getConeInfo. cf. clearGlobals) message } { cone.gw_ineq_projectedWt getConeInfo /cone.cinit set } ifelse (cone.cinit is --- the first number is the dim of cone.) messagen cone.cinit 0 get pmat % Maximal dimensional cone かどうかの検査. 検査にパスすれば loop を exit % パスしない場合 w_start を cone_random_vec を用いて変更する. cone.cinit 0 get 0 get to_int32 cone.m eq { exit } { (Failed to get the max dim cone. Updating the weight ...) messagen 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 } loop (cone.m = ) messagen cone.m message (Suceeded to get the maximal dimensional startingCone.) message % Linearity subspace の orth complement への射影行列. % 大域変数 cone.Lp, cone.Lpt を設定 cone.cinit 0 get 1 get /cone.Lp set cone.Lp transpose /cone.Lpt set % Linearity subspace の行列を設定. % 大域変数 cone.L を設定 cone.cinit 0 get 2 get /cone.L set % cone.d は cone.W および Linearity space で割った後, cone を考えるときの次元. % 大域変数 cone.d の設定. /cone.d cone.Lp length def cone.m cone.d eq { (There is no linearity space) message } { (Dim of the linearity space is ) messagen cone.m cone.d sub message (cone.Lp = ) messagen cone.Lp pmat } ifelse %[5] cone.g_ineq * cone.Wt * cone.Lpt % cone.w_ineq * cone.Wt * cone.Lpt % で制約を d 次元ベクトルに変換. % W (R^m) 空間の不等式制約を L' (R^d) 空間へ射影 % cone.gw_ineq_projectedWtLpt % = cone.g_ineq*cone.Wt*cone.Lpt \/ cone.w_ineq*coneWt*cone.Lpt /cone.gw_ineq_projectedWtLpt cone.gw_ineq_projectedWt cone.Lpt mul def cone.m cone.d eq { /cone.cinit.d cone.cinit def } { % cone.m > cone.d ならば, 再度 cone の計算が必要. % R^d の cone は cone.cinit.d へ入れる. cone.gw_ineq_projectedWtLpt getConeInfo /cone.cinit.d set } ifelse cone.cinit.d 1 get newCone /cone.startingCone set (cone.startingCone is ) message cone.startingCone message ] pop popVariables cone.startingCone } def % % data/test9.sm1 の test9 1-simplex X 2-simplex % % data/test10.sm1 1-simplex X 3-simplex % data/test11.sm1 SST, p.59 % % いよいよ, cone enumeration のプログラム書き開始 % %< % Usages: cone markBorder % cone->facets[i] が weight space の border にあるとき % cone->flipped[i] = 2 とする. % これを cone のすべての facet に対して計算. %> /markBorder { /arg1 set [/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 ] pop popVariables } def %< % Usages: ncone updateFan % グローバル変数 cone.fan を更新する. %> % % updateFan の debug は data/test8 でとりあえずやる. % test8 /ncone set を実行してから ncone updateFan % global: cone.fan /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 [/ncone /kk /cfacet /ii /jj /tcone /flipped_t] pushVariables [ /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 ncone cone.fan kk get isSameCone { (Internal error updateFan: ncone is already in cone.fan) error } { } ifelse } for % 1. ncone の中の border 上の facet をすべて mark. ncone markBorder % 2. ncone /\ cone.fan[kk] があるか調べる. あれば Mark する. incidence graph に加える 0 1 cone.fan.n 1 sub { /kk set ncone cone.fan kk get getCommonFacet /cfacet set cfacet 0 get { % 共通 facet がある場合. [[cone番号 face番号] [cone番号 face番号]] の形式で格納. /ii cfacet 1 get 0 get def /jj cfacet 2 get 0 get def cone.incidence [ [[cone.fan.n ii] [kk jj]] ] join /cone.incidence set % flipped を mark する. 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 を加える. cone.fan [ncone] join /cone.fan set ] pop popVariables } def %< % usages: getNextFlip [cone, k, cid] % cone.fan を検索して まだ flip してない cone と facet の組を戻す. % もうないときには null を戻す. % cid は cone が cone.fan の 何番目であるかの index. cone.gblist の検索等に % 用いる. %> /getNextFlip { [/tcone /ans /ii /cid] pushVariables [ /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 cid] def } ifelse ] pop popVariables arg1 } def % global variable : cone.epsilon , cone.epsilon.limit % 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.orig { /arg1 set [/ncone /ccone /kk /w /next_weight_w_wv] pushVariables [ /ccone arg1 def /ncone null def /kk ccone 1 get def 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.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 (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 (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 %< % Usages: set globals and getGrobnerFan % cf. clearGlobals % getStartingCone すると weightSpace とかの計算ができる. isOnWeightBorder が % 決められる. %> % とりあえず (data/test8.sm1) run してから getGrobnerFan /getGrobnerFan { getStartingCone /cone.ncone set { cone.ncone updateFan ( ) message (----------------------------------------------------------) message (getGrobnerFan #cone.fan=) messagen cone.fan length message cone.ncone /cone.ccone set getNextFlip /cone.nextflip set cone.nextflip tag 0 eq { exit } { } ifelse cone.nextflip getNextCone /cone.ncone set } loop (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