[BACK]Return to bfunc.sm1 CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Doc / Old

Annotation of OpenXM/src/kan96xx/Doc/Old/bfunc.sm1, Revision 1.1.1.1

1.1       maekawa     1: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      2: %      The b-function b_f(s),                                      %
                      3: %      generators of the annihilators of f^s,                      %
                      4: %      and the 1st algebraic local cohomology group                %
                      5: %      for f = f(x,y,z)                                            %
                      6: %                                  18 Dec. 1995  by T. Oaku        %
                      7: %                        modified  26 Feb. 1996                    %
                      8: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      9: 
                     10: (Type in <<  f bf3 >> for the b-function for f(x,y,z).) message
                     11: (Type in bf for the b-function for x^3-y^2 z^2.) message
                     12: (  ) message
                     13: /bf {(x^3 - y^2*z^2 ) bf3} def
                     14: %%%%%%%%%%%%% Template to compute b-function for f(x,y,z) %%%%%%%%%%%%
                     15: /bf3 {
                     16:   /f set
                     17:   %%% s is used both for F-homogenization and for tDt %%%%%%%
                     18:    [(s,t,x,y,z) ring_of_differential_operators
                     19:     [[(s) 1] ]
                     20:     weight_vector 0 ] define_ring
                     21:   %%% GIVE THE POLYNOMIAL f(x,y,z) HERE %%%%%%%%%%%%%%%%%%%%%
                     22:     f .  /f set
                     23:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     24:     f fw_delta /ff set
                     25:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     26:    (Computing the b-function of) message 
                     27:    f ::
                     28:    (The generators are) message
                     29:    ff ::
                     30:    ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                     31:    (Computing FW-groebner basis in Q[t,x,y,z]<Dt,Dx,Dy,Dz> ) message
                     32:    [ff] groebner 0 get /ansG set 
                     33:    (  ) message
                     34:   %%%%% ansG is an FW-Groebner basis in Q[t,x,y,z]<Dt,Dx,Dy,Dz> %%%%%%
                     35:    ansG {fw_symbol} map /ansG0 set 
                     36:    ansG0 {fw_psi} map /ansH set
                     37:   %%%%  ansH generates an ideal in Q[s,x,y,z]<Dx,Dy,Dz> %%%%%
                     38:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     39:    [(s,x,y,z) ring_of_differential_operators
                     40:     [[(Dx) 1 (Dy) 1 (Dz) 1] ] weight_vector  
                     41:    0
                     42:    ] define_ring
                     43:    ansH {mymap} map /ansH set
                     44:    ansH {[[(h). (1).]] replace} map {homogenize} map /ansH set
                     45:    (Eliminating Dx, Dy, Dz  ) message
                     46:    [ansH] groebner 0 get /ansH set 
                     47:    (  ) message
                     48:    ansH (Dx) eliminate0
                     49:         (Dy) eliminate0
                     50:         (Dz) eliminate0
                     51:    /ansH0 set
                     52:    ansH0 {[[(h). (1).]] replace} map /ansH01 set
                     53:    ansH0 {[[(s). (-s-1).]] replace} map /ansH0 set
                     54:    ansH0 minimal /ansH0 set
                     55: %%%% ansH0 generates an ideal in Q[s,x,y,z] %%%%
                     56: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     57:    [(s,x,y,z) ring_of_polynomials 
                     58:     (x,y,z) elimination_order  0] define_ring 
                     59:    ansH0 {mymap} map /ansH0 set
                     60:    ansH0 {[[(x). (0).] [(y). (0).] [(z). (0).]] replace} map /ansH00 set
                     61:  %%% ansH00 is the restriction of ansH0 to x=y=z=0 %%%
                     62:    ansH0 {[[(h). (1).]] replace} map {homogenize} map /ansH0 set 
                     63:    (Eliminating x,y,z ) message
                     64:    [ansH0] groebner 0 get /ansH0 set 
                     65:    ansH0 (x) eliminate0
                     66:          (y) eliminate0
                     67:          (z) eliminate0
                     68:    /ansbff set
                     69:    ansbff minimal /ansbff set
                     70:    ansbff 0 get /ansbf set
                     71:    (the global b-function b_f(s) [ansbf] is ) message 
                     72:    ansbf ::
                     73:  %%%%%% restriction to x=y=z=0 %%%%%%%%%%%%%%%%% 
                     74:    ansH00 remove0 /ansH00 set
                     75:    [(s) ring_of_polynomials 
                     76:     ( ) elimination_order  0] define_ring 
                     77:    ansH00 {mymap} map /ansH00 set
                     78:    ansH00 {[[(h). (1).]] replace} map {homogenize} map /ansH00 set 
                     79:    [ansH00] groebner 0 get /ansbff0 set
                     80:    ansbff0 minimal /ansbff0 set
                     81:    ansbff0  0 get /ansbf0 set
                     82:    (a divisor of the local b-function b_f(s) [ansbf0] is ) message 
                     83:    (  ) message
                     84:    ansbf0 ::
                     85: } def
                     86: 
                     87: %%%%%%%%%%%%% finding a P s.t. Pf^{s+1} = b_f(s)f^s %%%%%%%%%%%%%%%%%%%%%%%
                     88: 
                     89: /bf0 {
                     90:   %%% s is used both for F-homogenization and for tDt %%%%%%%
                     91:    [(s,t,x,y,z) ring_of_differential_operators
                     92:     [[(s) 1] [(Dx) 1 (Dy) 1 (Dz) 1 (x) 1 (y) 1 (z) 1]]
                     93:     weight_vector 0 ] define_ring
                     94:   %%% Give the polynomial f(x,y,z) here %%%%%%%%%%%%%%%%%%%%%
                     95:     ( x^5-y^2*z^2 ). /f set
                     96:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     97:     f fw_delta /ff set
                     98:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     99:    (Computing the b-function of) message 
                    100:    f ::
                    101:    (The generators are) message
                    102:    ff ::
                    103:    ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                    104:    (Computing FW-groebner basis in Q[t,x,y,z]<Dt,Dx,Dy,Dz> ) message
                    105:    [ff] groebner 0 get /ansG set 
                    106:    ansG {fw_order} map /ansGford set
                    107:    ansG {[[(h). (1).]] replace} map /ansG set
                    108:    ansG (Dx) eliminatepsi0
                    109:         (Dy) eliminatepsi0
                    110:         (Dz) eliminatepsi0
                    111:         (x) eliminatepsi0
                    112:         (y) eliminatepsi0
                    113:         (z) eliminatepsi0
                    114:    /ansbft set
                    115:    ansbft 0 get fw_rhorest /ansP set
                    116:    (Completed: P is [ansP]) ansP :: 
                    117:    (F) f toa
                    118:    (P) ansP toa
                    119:     ansbft
                    120:  } def
                    121: 
                    122: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    123: % 2nd algorithm for b-function
                    124: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    125: (Type in bf2 for the b-function via saturation.) message
                    126: (  ) message
                    127: /bf2 {
                    128:   %%% s is used both for F-homogenization and for t*Dt.
                    129:   %%% u is used for the computation of saturation.
                    130:    [(s,t,u,x,y,z) ring_of_differential_operators
                    131:     [[(s) 1 (u) 1]] weight_vector 
                    132:    0 ] define_ring
                    133:   %%% Write f(x) here.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    134:    ( y*(x^5- y^2*z^2) ). /f set
                    135:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    136:    (Computing the b-function of ) message 
                    137:    f ::  
                    138:    f  fw_deltasat /ff set
                    139:    ff print
                    140:    (  are generators.) message (   ) message
                    141:    ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                    142:    (Computing the saturation...) message
                    143:    [ff] groebner 0 get /ansS set 
                    144:    (     ) message
                    145:    ansS {[[(h). (1). ]] replace} map /ansS set
                    146:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    147:   ansS (s) eliminate0 
                    148:        (u) eliminate0 
                    149:   /ansS0 set
                    150:   ansS0 {fw_psi} map /ansS1 set
                    151:   ansS1 {[[(s). (-s-1).]] replace} map /ansS1 set
                    152:   ansS1 [f] concat /ansS1 set
                    153:   %%%%  ansS1 generates an ideal in Q[s,x,y,z]<Dx,Dy,Dz> %%%%%
                    154:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    155:    [(s,x,y,z) ring_of_differential_operators
                    156:     [[(Dx) 1 (Dy) 1 (Dz) 1] [(x) 1 (y) 1 (z) 1]] weight_vector  
                    157:    0
                    158:    ] define_ring
                    159:    ansS1 {mymap} map /ansS1 set
                    160:    ansS1 {[[(h). (1).]] replace} map {homogenize} map /ansS1 set
                    161:    (Eliminating Dx, Dy, Dz  ) message
                    162:    [ansS1] groebner 0 get /ansS1 set 
                    163:    (  ) message
                    164:    ansS1 (Dx) eliminate0
                    165:          (Dy) eliminate0
                    166:          (Dz) eliminate0
                    167:    /ansJ set
                    168: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    169:    [(s,x,y,z) ring_of_polynomials 
                    170:     (x,y,z) elimination_order  0] define_ring 
                    171:    ansJ {mymap} map /ansJ set
                    172:    ansJ {[[(x). (0).] [(y). (0).] [(z). (0).]] replace} map /ansJ0 set
                    173:  %%% ansJ0 is the restriction of ansJ to x=y=z=0 %%%
                    174:    ansJ {[[(h). (1).]] replace} map {homogenize} map /ansJ set 
                    175:    (Eliminating x,y,z ) message
                    176:    [ansJ] groebner 0 get /ansJ set 
                    177:    ansJ  (x) eliminate0
                    178:          (y) eliminate0
                    179:          (z) eliminate0
                    180:    /ansbffS set
                    181:    ansbffS minimal /ansbffS set
                    182:    ansbffS 0 get /ansbfS set
                    183:    (the global b-function b_f(s) [ansbfS] is ) message 
                    184:    ansbfS ::
                    185:  %%%%%% restriction to x=y=z=0 %%%%%%%%%%%%%%%%% 
                    186:    ansJ0 remove0 /ansJ0 set
                    187:    [(s) ring_of_polynomials 
                    188:     ( ) elimination_order  0] define_ring 
                    189:    ansJ0 {mymap} map /ansJ0 set
                    190:    ansJ0 {[[(h). (1).]] replace} map {homogenize} map /ansJ0 set 
                    191:    [ansJ0] groebner 0 get /ansbffS0 set
                    192:    ansbffS0 minimal /ansbffS0 set
                    193:    ansbffS0  0 get /ansbfS0 set
                    194:    (a divisor of the local b-function b_f(s) [ansbfS0] is ) message 
                    195:    (  ) message
                    196:    ansbfS0 ::
                    197: } def
                    198: 
                    199: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    200: (Type in fs for the annihilators of f^s.)  message
                    201: (  ) message
                    202: %%%%%%%%%%%%%% Computing the annihilators of f^s %%%%%%%%%%%%%%%%%%%%
                    203: /fs {
                    204:   %%% s is used both for F-homogenization and for t*Dt.
                    205:   %%% u is used for the computation of saturation.
                    206:    [(s,t,u,x,y,z) ring_of_differential_operators
                    207:     [[(s) 1 (u) 1]] weight_vector 
                    208:    0 ] define_ring
                    209:   %%% Write f(x) here.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    210:    ( y*(x^5-y^2*z^2) ). /f set
                    211:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    212:    (Computing involutory generators for f^s with f = ) message 
                    213:    f ::  
                    214:    f  fw_deltasat /ff set
                    215:    ff print
                    216:    (  are generators.) message (   ) message
                    217:    ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                    218:    (Computing groebner basis) message
                    219:    [ff] groebner 0 get /ans set 
                    220:    (     ) message
                    221:    ans {[[(h). (1).]] replace} map /ans0 set
                    222:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    223:   ans0 (s) eliminate0 
                    224:        (u) eliminate0 
                    225:   /ans1 set
                    226:   ans1 {fw_psi} map /ans1 set
                    227:   ans1 {[[(s). (-s-1).]] replace} map /ans1 set
                    228:   ans1 involutory /ans2 set
                    229:   ans2 minimal /ansfs set
                    230:   (The answer [ansfs] is ) message
                    231:   ansfs print (  ) message
                    232:   (F) f toa
                    233:   (FS) ansfs toa_l
                    234: } def
                    235: 
                    236: %% Computing the D-module for f^s as D-module (not as D[s]-module)
                    237: /fs0 {
                    238:   %%% s is used both for F-homogenization and for t*Dt.
                    239:   %%% u is used for the computation of saturation.
                    240:    [(s,t,u,x,y,z) ring_of_differential_operators
                    241:     [[(s) 1 (u) 1]] weight_vector 
                    242:    0 ] define_ring
                    243:   %%% Write f(x) here.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    244:    ( y*(x^5-y^2*z^2) ). /f set
                    245:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    246:    (Computing involutory generators for f^s with f = ) message 
                    247:    f ::  
                    248:    f  fw_deltasat /ff set
                    249:    ff print
                    250:    (  are generators.) message (   ) message
                    251:    ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                    252:    (Computing groebner basis) message
                    253:    [ff] groebner 0 get /ans set 
                    254:    (     ) message
                    255:    ans {[[(h). (1).]] replace} map /ans0 set
                    256:   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    257:   ans0 (s) eliminate0 
                    258:        (u) eliminate0 
                    259:   /ans1 set
                    260:   ans1 {fw_psi} map /ans1 set
                    261:   ans1 {[[(s). (-s-1).]] replace} map /ans1 set
                    262:   
                    263:   ans1 {[[(h). (1).]] replace} map /ans1 set
                    264:   ans1 {homogenize} map /ans1 set 
                    265:   [ans1] groebner 0 get /ans1 set
                    266:   ans1 {[[(h). (1).]] replace} map /ans1 set
                    267:   ans1 (s) eliminate0 /ans1 set
                    268:   ans1 involutory /ans2 set
                    269:   ans2 minimal /ansfs set
                    270:   (The answer [ansfs] is ) message
                    271:   ansfs print (  ) message
                    272:   (F) f toa
                    273:   (FS) ansfs toa_l
                    274: } def
                    275: 
                    276: %%%% algebraic local cohomology %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    277: (Type in alc for the 1st algebraic local cohomology group.) message
                    278: ((Make sure for alc that b_f(s) has no integer roots other than -1.)) message
                    279: (  ) message
                    280: /alc {
                    281:   %%% s is used for FW-filtration.
                    282:    [(s,t,x,y,z) ring_of_differential_operators
                    283:     [[(s) 1]] weight_vector  0 ] define_ring
                    284: %%% give the polynomial f(x,y,z) here %%%%%%%%%%%%%%%%%%%%%%%%
                    285:   ( x^3 + y^3 + z^3 ). /f set
                    286: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    287:   f fw_delta /ff set
                    288: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    289:   (Computing the algebraic local cohomology for) message 
                    290:   f ::
                    291:   ff {[[(h). (1).]] replace} map {homogenize} map /ff set
                    292:   (Computing an FW-groebner basis) message
                    293:   [ff] groebner 0 get /ansfw set 
                    294:   (     ) message
                    295: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    296: % selecting the elements of F-order 0
                    297: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    298:   /gb ansfw def
                    299:   gb {fw_order} map /gbford set
                    300:   
                    301:   /ansind0 [
                    302:     0 1 << gb length 1 sub >> {
                    303:       /n set 
                    304:       gb n get /ff set
                    305:       ff fw_order (integer) data_conversion /m set 
                    306:       << m 2 lt >> 
                    307:        { << m 1 1 >> {pop ff (Dt). ff mul /ff set } for 
                    308:        }
                    309:        {    } ifelse
                    310:     } for
                    311:   ] def
                    312:   ansind0 {[[(h). (1).]] replace} map /ansind0 set
                    313:   ansind0 {[[(s). (1).]] replace} map /ansind0 set
                    314:   ansind0 {[[(t). (0).]] replace} map /ansind1 set
                    315:   ansind1 remove0 /ansind1 set
                    316:   ansind1 involutory /ansind2 set
                    317:   ansind2 minimal /ansind set
                    318:   (The answer [ansind] is ) message
                    319:   ansind ::
                    320:   (F) f toa
                    321:   (ALC) ansind toa_l
                    322: } def
                    323: 
                    324: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    325: %%%% involutory base in K(s)[x,y,z]<Dx,Dy,Dz> 
                    326: /involutory  {
                    327:    /ansff0 set
                    328:    [/ansff1 /ansff2 /ansff3 ] pushVariables
                    329:    [
                    330:      [(s,t,x,y,z) ring_of_differential_operators
                    331:       [[(Dx) 1 (Dy) 1 (Dz) 1] [ (x) 1 (y) 1 (z) 1 ]] weight_vector
                    332:      0
                    333:      ] define_ring
                    334:      ansff0 {mymap} map /ansff1 set
                    335:      ansff1 {[[(h). (1).]] replace} map {homogenize} map /ansff2 set
                    336:      (Computing an involutory base ) message
                    337:      [ansff2] groebner 0 get /ansff3 set
                    338:      (  ) message
                    339:      ansff3 {[[(h). (1).]] replace} map /ansff3 set
                    340:      ansff3 /ansinv set
                    341:    ] pop
                    342:    popVariables
                    343:    ansinv
                    344: } def
                    345: 
                    346: %%%%%% for FW-filtration, etc. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    347: % F-order
                    348: /fw_order {
                    349:   fw_symbol /fws set
                    350:   fws [[(s). (1).]] replace /fws set
                    351:   fws (Dt). coefficients 0 get 0 get /m set
                    352:   fws (Dt). coefficients 1 get 0 get /fwsDt set
                    353:   fwsDt (t). coefficients 0 get 0 get /k set
                    354:   m k sub
                    355: } def
                    356: 
                    357: % remove 0 elements from a list
                    358: /remove0 {
                    359:   /arg1 set
                    360:   [/gb /ans /n] pushVariables
                    361:   [
                    362:   /gb arg1 def
                    363:   /ans [
                    364:     0 1 << gb length 1 sub >> {
                    365:       /n set 
                    366:       gb n get /ff set
                    367:       ff (0). eq
                    368:         {  }
                    369:         { ff } ifelse
                    370:     } for
                    371:   ] def
                    372:   /arg1 ans def
                    373:   ] pop
                    374:   popVariables
                    375:   arg1
                    376: } def
                    377: 
                    378: % dehomogenize and obtain a minimal base
                    379: % in variables s,t,x,y,z,Dt,Dx,Dy,Dz 
                    380: (Note that the current ring changes once you get a minimal base.) message
                    381: /minimal {
                    382:   /arg1 set
                    383:   [/gb /inits /ans /n /syzlist /cc /nin /aa /j /cj] pushVariables
                    384:   [
                    385:   /gb arg1 def
                    386:   /inits gb {init} map def
                    387:   gb {[[(h). (1).]] replace} map /gb set
                    388:   [(Dx,Dy,Dz,Dt,t,x,y,z,s) ring_of_polynomials  
                    389:    (  ) elimination_order  0] define_ring
                    390:   inits {mymap} map /inits set
                    391:   gb    {mymap} map /gb    set
                    392:   inits {[[(h). (1).]] replace} map /inits set
                    393:   inits length /nin set
                    394:   [inits [(needBack)]] groebner 1 get /syzlist set
                    395:   syzlist length ::
                    396:   /ans [
                    397:     0 1 << syzlist length 1 sub >> {
                    398:       /n set
                    399:       syzlist n get /cc set
                    400:       (0). /gg set
                    401:       0 1 << nin 1 sub >> {
                    402:         /j set
                    403:         gb j get /aa set
                    404:         cc j get /cj set
                    405:         << cj aa mul >> gg add /gg set
                    406:       } for
                    407:       gg
                    408:     } for
                    409:   ] def
                    410:   /ansmin ans def
                    411:   ] pop
                    412:   popVariables
                    413:   ansmin 
                    414: } def
                    415: 
                    416: %%%%% The formal symbol %%%%%%%%%%%%%%%%%%%%%%
                    417: /fw_symbol {
                    418:   [[(h). (1).]] replace (s). coefficients 1 get 0 get
                    419: } def
                    420: 
                    421: %%%%% psi(P)(s) %%%%%%
                    422: /fw_psi {
                    423:   fw_symbol /P set
                    424:   P fw_order (integer) data_conversion /k set
                    425:     << 1 1 k >> 
                    426:     {(t). P mul /P set pop}
                    427:     for
                    428:     << -1 -1 k >>
                    429:     {(Dt). P mul /P set pop}
                    430:     for 
                    431:   (0). /Q set
                    432:   P (Dt). coefficients 0 get length /m set
                    433:   0 /i set
                    434:   1 1 m  
                    435:   {
                    436:     P (Dt). coefficients 0 get i get /kk set 
                    437:     P (Dt). coefficients 1 get i get /PPt set
                    438:     PPt (t). coefficients 1 get 0 get /PPC set
                    439:     kk (integer) data_conversion /kk set
                    440:     (s). /Ss set
                    441:     0 1 << kk 1 sub >> {
                    442:       PPC Ss mul /PPC set
                    443:       Ss (1). sub /Ss set
                    444:       pop
                    445:     } for
                    446:     Q PPC add /Q set
                    447:     i 1 add /i set
                    448:     pop
                    449:   } for
                    450:   Q  
                    451: } def
                    452: 
                    453: /fw_psi0 {
                    454:   fw_symbol /P set
                    455:   P fw_order (integer) data_conversion /k set
                    456:     << 1 1 k >> 
                    457:     {(t). P mul /P set pop}
                    458:     for
                    459:   (0). /Q set
                    460:   P (Dt). coefficients 0 get length /m set
                    461:   0 /i set
                    462:   1 1 m  
                    463:   {
                    464:     P (Dt). coefficients 0 get i get /kk set 
                    465:     P (Dt). coefficients 1 get i get /PPt set
                    466:     PPt (t). coefficients 1 get 0 get /PPC set
                    467:     kk (integer) data_conversion /kk set
                    468:     (s). /Ss set
                    469:     0 1 << kk 1 sub >> {
                    470:       PPC Ss mul /PPC set
                    471:       Ss (1). sub /Ss set
                    472:       pop
                    473:     } for
                    474:     Q PPC add /Q set
                    475:     i 1 add /i set
                    476:     pop
                    477:   } for
                    478:   Q  
                    479: } def
                    480: 
                    481: %%%%% rho(P)(s) %%%%%%
                    482: /fw_rho {
                    483:   /P0 set
                    484:   P0 fw_order (integer) data_conversion /k set
                    485:   << 1 1 k >> 
                    486:     {(t). P0 mul /P0 set}
                    487:   for
                    488:   << -1 -1 k >>
                    489:     {(Dt). P0 mul /P0 set}
                    490:   for 
                    491:   P0 (s). coefficients 0 get /sdegs set
                    492:   sdegs length /n set
                    493:   sdegs 0 get (integer) data_conversion /k0 set
                    494:   (0). /PP set 
                    495:   0 /jj set
                    496:   0 1 << n 1 sub >>
                    497:   {
                    498:     sdegs jj get (integer) data_conversion /kkp set
                    499:     P0 (s). coefficients 1 get jj get /Pj set
                    500:     Pj fw_psi0 /Pj set
                    501:     << 1 1 << k0 kkp sub >> >>    
                    502:       {Pj f mul /Pj set pop} for
                    503:     k0 kkp sub  /l set
                    504:     Pj [[(s). << (-s-1). l . sub >> ]] replace /Pj set
                    505:     PP Pj add /PP set 
                    506:     jj 1 add /jj set
                    507:     pop
                    508:   } for 
                    509:   PP [[(h). (1).]] replace /PP set 
                    510:   pop 
                    511:   PP
                    512: } def 
                    513: 
                    514: /fw_rhorest {
                    515:   /P0 set
                    516:   P0 fw_order (integer) data_conversion /k set
                    517:   << 1 1 k >> 
                    518:     {(t). P0 mul /P0 set}
                    519:   for
                    520:   << -1 -1 k >>
                    521:     {(Dt). P0 mul /P0 set}
                    522:   for 
                    523:   P0 (s). coefficients 0 get /sdegs set
                    524:   sdegs length /n set
                    525:   sdegs 0 get (integer) data_conversion /k0 set
                    526:   (0). /PP set 
                    527:   1 /jj set
                    528:   1 1 << n 1 sub >>
                    529:   { 
                    530:     sdegs jj get (integer) data_conversion /kkp set
                    531:     P0 (s). coefficients 1 get jj get /Pj set
                    532:     Pj fw_psi0 /Pj set
                    533:      << 2 1 << k0 kkp sub >> >>    
                    534:       {Pj f mul /Pj set } for
                    535:     k0 kkp sub  /l set
                    536:     Pj [[(s). << (-s-1). l . sub >> ]] replace /Pj set
                    537:     PP Pj add /PP set
                    538:     jj 1 add /jj set
                    539:    } for 
                    540:   PP [[(h). (1).]] replace /PP set
                    541:   (-1). PP mul /PP set 
                    542:   pop 
                    543:   PP
                    544: } def 
                    545: 
                    546: /fw_rhotest {
                    547:    [(s,t,x,y,z) ring_of_differential_operators
                    548:     [[(s) 1] ]
                    549:     weight_vector 0 ] define_ring
                    550:     (x^2-y^3). /f set 
                    551:     (t*Dt^2*s + t*Dt). /Pex  set
                    552:     Pex fw_rho /PPex set
                    553:     PPex ::
                    554: } def 
                    555: 
                    556:  %%%%%%%%%%%% [t - s*f, Dx + f_xDt, ...] %%%%%%%%%%%%%%%
                    557: /fw_delta {
                    558:   /F set
                    559:   << (Dx). F mul >> << F (Dx). mul >> sub [[(h). (1).]] replace /Fx set
                    560:   << (Dy). F mul >> << F (Dy). mul >> sub [[(h). (1).]] replace /Fy set
                    561:   << (Dz). F mul >> << F (Dz). mul >> sub [[(h). (1).]] replace /Fz set
                    562:   (t). << (s). F mul >>  sub /F0 set
                    563:   (Dx). << (s*Dt). Fx mul >> add /Fx set 
                    564:   (Dy). << (s*Dt). Fy mul >> add /Fy set
                    565:   (Dz). << (s*Dt). Fz mul >> add /Fz set
                    566:   [ F0 Fx Fy Fz ]
                    567: } def
                    568: 
                    569:  %%%%%%%%%%%% [1-s*u,t - s*f, Dx + f_xDt, ...] %%%%%%%%%%%%%%%
                    570: /fw_deltasat {
                    571:   /F set
                    572:   << (Dx). F mul >> << F (Dx). mul >> sub [[(h). (1).]] replace /Fx set
                    573:   << (Dy). F mul >> << F (Dy). mul >> sub [[(h). (1).]] replace /Fy set
                    574:   << (Dz). F mul >> << F (Dz). mul >> sub [[(h). (1).]] replace /Fz set
                    575:   (t). << (s). F mul >>  sub /F0 set
                    576:   (Dx). << (s*Dt). Fx mul >> add /Fx set 
                    577:   (Dy). << (s*Dt). Fy mul >> add /Fy set
                    578:   (Dz). << (s*Dt). Fz mul >> add /Fz set
                    579:   [ F0 Fx Fy Fz (1-s*u). ]
                    580: } def
                    581: 
                    582: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    583: %   convert to a Risa/Asir input file         %
                    584: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
                    585: % (Varriable name) expression  toa
                    586: /toa {
                    587:   /expr set /Varname set
                    588:   expr (string) data_conversion /expr set
                    589:   (toa.t) (a) file /fd set
                    590:   fd (Dx = dx$ Dy = dy$ Dz = dz$ Dt = dt$) writestring
                    591:   fd Varname writestring
                    592:   fd ( = ) writestring
                    593:   fd expr  writestring
                    594:   fd ($) writestring
                    595:   fd ( ) writestring
                    596:   fd closefile
                    597: } def
                    598: 
                    599: % (Varriable name) expression(list)  toa_l
                    600: /toa_l {
                    601:   /expr set /Varname set
                    602:   (toa.t) (a) file /fd set
                    603:   fd (Dx = dx$ Dy = dy$ Dz = dz$ Dt = dt$) writestring
                    604:   fd 10 (string) data_conversion writestring
                    605:   fd Varname writestring
                    606:   fd ( = [ ) writestring
                    607:   fd 10 (string) data_conversion writestring
                    608:   expr length /n set
                    609:   0 1 << n 1 sub >> {
                    610:     /k set
                    611:     expr k get /expr1 set
                    612:     expr1 (string) data_conversion /expr1 set
                    613:     fd expr1  writestring
                    614:     k << n 1 sub >> eq 
                    615:       {fd (]$ )  writestring }
                    616:       {fd ( , )  writestring } 
                    617:     ifelse 
                    618:     fd 10 (string) data_conversion writestring
                    619:   } for  
                    620:   fd closefile
                    621: } def
                    622: 
                    623: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    624: 
                    625: /mymap {
                    626:   (string) data_conversion .
                    627: } def
                    628: 
                    629: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    630: %% [   ] {outputans1} map ;
                    631: /outputans1 {
                    632:  (t.t) (a) file /fd set
                    633:  (string) data_conversion /tmp0 set
                    634:  fd tmp0 writestring
                    635:  fd (  ,) writestring
                    636:  fd 10 (string) data_conversion writestring
                    637:  fd closefile
                    638: } def
                    639: 
                    640: %%%%%%%% Do not touch the below. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    641: [(position)
                    642:   [(set element position number)
                    643:    (Example: [(cat) (dog) (hot chocolate)] (cat) position ===> 0.)
                    644:   ]
                    645: ] putUsages
                    646: /position {
                    647:   /arg2 set /arg1 set
                    648:   [/univ /elem /num /flag] pushVariables
                    649:   [
                    650:      /univ arg1 def 
                    651:      /elem arg2 def
                    652:      /num -1 def /flag -1 def
                    653:      0 1 << univ length 1 sub >> 
                    654:      {
                    655:         /num set
                    656:         univ num get  elem  eq
                    657:         { /flag 0 def exit }
                    658:         {    }
                    659:         ifelse   
                    660:      }  for
                    661:      flag -1 eq
                    662:      {/num -1 def}
                    663:      {  }
                    664:      ifelse
                    665:   ] pop
                    666:   /arg1 num def
                    667:   popVariables
                    668:   arg1
                    669: } def
                    670: 
                    671: 
                    672: [(evecw)
                    673:   [(size position weight evecw  [0 0 ... 0 weight 0 ... 0] )
                    674:    (Example: 3 0 113 evecw ===> [113  0  0])
                    675:   ]
                    676: ] putUsages
                    677: /evecw {
                    678:  /arg3 set /arg2 set /arg1 set
                    679:  [/size /iii /www] pushVariables
                    680:  /size arg1 def  /iii arg2 def /www arg3 def
                    681:  [
                    682:    0 1 << size 1 sub >>
                    683:    {
                    684:       iii eq
                    685:       {  www }
                    686:       {  0 }
                    687:       ifelse
                    688:    } for
                    689:   ] /arg1 set
                    690:   popVariables
                    691:   arg1
                    692: } def
                    693: 
                    694: [(weight_vector)
                    695:  [ ([x-list d-list params] [[(name) weight ...] [...] ...] weight_vector)
                    696:    ([x-list d-list params order])
                    697:    (Example:)
                    698:    (   [(x,y,z) ring_of_polynomials [[(x) 100 (y) 10]] weight_vector 0] )
                    699:    (   define_ring )
                    700:   ]
                    701: ] putUsages
                    702: /weight_vector {
                    703:   /arg2 set  /arg1 set
                    704:   [/vars /univ /w-vectors /www /k /order1 /order2] pushVariables
                    705:   /vars arg1 def /w-vectors arg2 def
                    706:   [
                    707:     /univ vars 0 get reverse
                    708:           vars 1 get reverse join
                    709:     def
                    710:     [
                    711:     0 1 << w-vectors length 1 sub >> 
                    712:     {
                    713:       /k set
                    714:       univ w-vectors k get w_to_vec
                    715:     } for
                    716:     ] /order1 set
                    717:     %% order1 ::
                    718:     
                    719:     vars ( ) elimination_order 3 get /order2 set
                    720:     vars [ << order1 order2 join >> ] join /arg1 set
                    721:   ] pop
                    722:   popVariables
                    723:   arg1
                    724: } def
                    725: 
                    726: %% [(e) (x) (y) (h)] [(x) 100 (y) 10] w_to_vec [0 100 10 0]
                    727: %%  univ              www
                    728: /w_to_vec {
                    729:   /arg2 set  /arg1 set
                    730:   [/univ /www /k /vname /vweight /ans] pushVariables
                    731:   /univ arg1 def /www arg2 def
                    732:   [ 
                    733:     /ans << univ length >> -1 0 evecw def
                    734:     0  2  << www length 2 sub >>
                    735:     {
                    736:       %% ans ::
                    737:       /k set
                    738:       www k get /vname set
                    739:       www << k 1 add >> get /vweight set
                    740:       << univ length >> 
                    741:       << univ vname position >>
                    742:       vweight evecw
                    743:       ans add /ans set
                    744:     } for
                    745:     /arg1 ans def
                    746:   ] pop
                    747:   popVariables
                    748:   arg1
                    749: } def
                    750: 
                    751: 
                    752: /fw_principal {
                    753:    {[[(h). (1).]] replace} map {(s). coefficients 1 get 0 get} map
                    754: } def
                    755: 
                    756: 
                    757: %%%%%%%%%%%%%%%%%%%%%
                    758: % [g1 g2 g3 ...] var eliminate0
                    759: /eliminate0 {
                    760:   /arg2 set /arg1 set
                    761:   [/gb /degs /ans /n /var] pushVariables
                    762:   [
                    763:   /gb arg1 def
                    764:   /var arg2 def
                    765:   /degs gb {var . degree} map def
                    766:   /ans [
                    767:     0 1 << gb length 1 sub >> {
                    768:       /n set
                    769:       << degs n get  >>  0 eq
                    770:       { gb n get /ff set
                    771:         ff (0). eq
                    772:         {  }
                    773:         { ff } ifelse
                    774:       }
                    775:       {   } ifelse
                    776:     } for
                    777:   ] def
                    778:   /arg1 ans def
                    779:   ] pop
                    780:   popVariables
                    781:   arg1
                    782: } def
                    783: 
                    784: % [g1 g2 g3 ...] var eliminate0
                    785: /eliminatepsi0 {
                    786:   /arg2 set /arg1 set
                    787:   [/gb /degs /ans /n /var] pushVariables
                    788:   [
                    789:   /gb arg1 def
                    790:   /var arg2 def
                    791:   /degs gb {fw_symbol} map {var . degree} map def
                    792:   /ans [
                    793:     0 1 << gb length 1 sub >> {
                    794:       /n set
                    795:       << degs n get  >>  0 eq
                    796:       { gb n get /ff set
                    797:         ff (0). eq
                    798:         {  }
                    799:         { ff } ifelse
                    800:       }
                    801:       {   } ifelse
                    802:     } for
                    803:   ] def
                    804:   /arg1 ans def
                    805:   ] pop
                    806:   popVariables
                    807:   arg1
                    808: } def
                    809: 
                    810: %%%%%%%%% concatenate two lists %%%%%%% 
                    811: /concat {
                    812:   /listB set /listA set
                    813:   listA length /NA set
                    814:   listB length /NB set
                    815:   /listAB [
                    816:     0 1 << NA 1 sub >> {
                    817:       /n set
                    818:       listA n get
                    819:     } for
                    820:    0 1 << NB 1 sub >> {
                    821:       /n set
                    822:       listB n get
                    823:     } for
                    824:   ] def
                    825:   listAB
                    826: } def 

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>