Annotation of OpenXM/src/kan96xx/Doc/hol.sm1, Revision 1.1.1.1
1.1 maekawa 1: %% hol.sm1, 1998, 11/8, 11/10, 11/14, 11/25, 1999, 5/18, 6/5.
2: %% rank, rrank, characteristic
3: %% This file is error clean.
4: /hol.version (2.990515) def
5: hol.version [(Version)] system_variable gt
6: { [(This package hol.sm1 requires the latest version of kan/sm1) nl
7: (Please get it from http://www.math.kobe-u.ac.jp/KAN)
8: ] cat
9: error
10: } { } ifelse
11:
12: $hol.sm1, basic package for holonomic systems (C) N.Takayama, 1999, 6/05 $
13: message-quiet
14:
15: /rank.v [(x) (y) (z)] def %% default value of v (variables).
16: /rank.ch [ ] def %% characteristic variety.
17: /rank.verbose 0 def
18: /rank {
19: /arg1 set
20: [/in-rank /aa /typev /setarg /f /v /vsss /vddd
21: /gg /wv /vd /vdweight /chv
22: /one
23: ] pushVariables
24: [(CurrentRingp) (KanGBmessage)] pushEnv
25: [
26:
27: /aa arg1 def
28: aa isArray { } { ( << array >> rank) error } ifelse
29: /setarg 0 def
30: aa { tag } map /typev set
31: typev [ ArrayP ] eq
32: { /f aa 0 get def
33: /v rank.v def
34: /setarg 1 def
35: } { } ifelse
36: typev [ArrayP StringP] eq
37: { /f aa 0 get def
38: /v [ aa 1 get to_records pop ] def
39: /setarg 1 def
40: } { } ifelse
41: typev [ArrayP ArrayP] eq
42: { /f aa 0 get def
43: /v aa 1 get def
44: /setarg 1 def
45: } { } ifelse
46: setarg { } { (rank : Argument mismatch) error } ifelse
47:
48: [(KanGBmessage) rank.verbose ] system_variable
49:
50: f { toString } map /f set
51: v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map
52: /vddd set %% vddd = [(Dx) 1 (Dy) 1 (Dz) 1]
53: v { @@@.Dsymbol 2 1 roll 2 cat_n } map
54: /vd set %% vd = [(Dx) (Dy) (Dz)]
55: /vdweight
56: vd { [ 2 1 roll -1 ] } map %% vdweight=[[(Dx) -1] [(Dy) -1] [(Dz) -1]]
57: def
58:
59: [v from_records
60: ring_of_differential_operators [vddd] weight_vector 0] define_ring
61: f { . dehomogenize } map /f set
62: [f] groebner_sugar 0 get /gg set
63:
64: /wv vddd weightv def
65: gg { wv init } map /chv set %%obtained the characteristic variety.
66: /rank.ch chv def
67: chv { toString } map /chv set
68:
69: [ v vd join from_records
70: ring_of_polynomials
71: [vddd] vdweight join weight_vector
72: 0
73: ] define_ring
74: [chv {.} map] groebner_sugar 0 get { init } map /chii set
75:
76: /rank.chii chii def
77: rank.verbose { chii message } { } ifelse
78: v {[ 2 1 roll . (1).]} map /one set
79: %% [[(x). (1).] [(y). (1).] [(z). (1).]]
80: %% chii { one replace } map %% buggy code.
81: %% Arg of hilb should be a reduced GB.
82: [chii { one replace } map] groebner 0 get
83: vd hilb /arg1 set
84: ] pop
85: popEnv
86: popVariables
87: arg1
88: } def
89:
90:
91: [(rank)
92: [( a rank b)
93: ( array a; number b)
94: (Example 1 : )
95: $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] rank :: $
96: (Example 2 : )
97: $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] rank :: $
98: ]
99: ] putUsages
100: (rank ) messagen-quiet
101:
102: /characteristic.verbose 0 def
103: /characteristic.v [(x) (y) (z)] def
104: /characteristic.ch [ ] def
105: /ch { characteristic } def
106: /characteristic {
107: /arg1 set
108: [/in-rank /aa /typev /setarg /f /v /vsss /vddd
109: /gg /wv /vd /chv
110: /one
111: ] pushVariables
112: [(CurrentRingp) (KanGBmessage)] pushEnv
113: [
114:
115: /aa arg1 def
116: aa isArray { } { ( << array >> characteristic) error } ifelse
117: /setarg 0 def
118: aa { tag } map /typev set
119: typev [ ArrayP ] eq
120: { /f aa 0 get def
121: /v characteristic.v def
122: /setarg 1 def
123: } { } ifelse
124: typev [ArrayP StringP] eq
125: { /f aa 0 get def
126: /v [ aa 1 get to_records pop ] def
127: /setarg 1 def
128: } { } ifelse
129: typev [ArrayP ArrayP] eq
130: { /f aa 0 get def
131: /v aa 1 get def
132: /setarg 1 def
133: } { } ifelse
134: setarg { } { (rank : Argument mismatch) error } ifelse
135:
136: [(KanGBmessage) characteristic.verbose ] system_variable
137:
138: f { toString } map /f set
139: v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map
140: /vddd set %% vddd = [(Dx) 1 (Dy) 1 (Dz) 1]
141: v { @@@.Dsymbol 2 1 roll 2 cat_n } map
142: /vd set %% vd = [(Dx) (Dy) (Dz)]
143:
144: [v from_records
145: ring_of_differential_operators [vddd] weight_vector 0] define_ring
146: f { . dehomogenize } map /f set
147: [f] groebner_sugar 0 get /gg set
148:
149: /wv vddd weightv def
150: gg { wv init } map /chv set
151: /characteristic.ch [chv] def
152: %% gg { wv init toString} map /chv set %%obtained the characteristic variety.
153: %% /characteristic.ch chv def
154:
155: %% [ v vd join from_records
156: %% ring_of_polynomials
157: %% [vddd] weight_vector
158: %% 0
159: %% ] define_ring
160: %% [chv {.} map] groebner_sugar 0 get /characteristic.ch set
161:
162: characteristic.ch /arg1 set
163: ] pop
164: popEnv
165: popVariables
166: arg1
167: } def
168:
169: [(characteristic)
170: [( a characteristic b)
171: ( array a; number b)
172: (b is the generator of the characteristic variety of a.)
173: (For the algorithm, see Japan J. of Industrial and Applied Math., 1994, 485--497.)
174: (Example 1 : )
175: $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] characteristic :: $
176: (Example 2 : )
177: $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] characteristic :: $
178: ]
179: ] putUsages
180: (characteristic ) messagen-quiet
181: [(ch)
182: [(ch is the abbreviation of characteristic.)
183: ( a ch b)
184: ( array a; number b)
185: (b is the generator of the characteristic variety of a.)
186: (For the algorithm, see, Japan J. of Industrial and Applied Math., 1994, 485--497.)
187: (Example 1 : )
188: $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] ch :: $
189: (Example 2 : )
190: $[ [( (x^3-y^2) Dx + 3 x^2) ( (x^3-y^2) Dy - 2 y)] (x,y)] ch :: $
191: ]
192: ] putUsages
193: (ch ) messagen-quiet
194:
195: %%%% developing rrank.sm1
196: /rrank.v [(x) (y) (z)] def %% default value of v (variables).
197: /rrank.init [ ] def %% initial ideal.
198: /rrank.verbose 0 def
199: /rrank {
200: /arg1 set
201: [/in-rrank /aa /typev /setarg /f /v /vsss /vddd
202: /gg /wv /vd /vdweight
203: /one /i /chv
204: ] pushVariables
205: [(CurrentRingp) (KanGBmessage)] pushEnv
206: [
207:
208: /aa arg1 def
209: aa isArray { } { ( << array >> rrank) error } ifelse
210: /setarg 0 def
211: aa { tag } map /typev set
212: typev [ ArrayP ] eq
213: { /f aa 0 get def
214: /v rrank.v def
215: /setarg 1 def
216: } { } ifelse
217: typev [ArrayP StringP] eq
218: { /f aa 0 get def
219: /v [ aa 1 get to_records pop ] def
220: /setarg 1 def
221: } { } ifelse
222: typev [ArrayP ArrayP] eq
223: { /f aa 0 get def
224: /v aa 1 get def
225: /setarg 1 def
226: } { } ifelse
227: setarg { } { (rrank : Argument mismatch) error } ifelse
228:
229: [(KanGBmessage) rrank.verbose ] system_variable
230:
231: f { toString } map /f set
232: v { @@@.Dsymbol 2 1 roll 2 cat_n 1 } map
233:
234: v { @@@.Dsymbol 2 1 roll 2 cat_n } map
235: /vd set %% vd = [(Dx) (Dy) (Dz)] , v = [(x) (y) (z)]
236: /vdweight
237: [ 0 1 v length 1 sub { /i set v i get << 0 i sub >>
238: vd i get << i >> } for ]
239: def
240: rrank.verbose { vdweight message } { } ifelse
241:
242: [v from_records
243: ring_of_differential_operators [vdweight] weight_vector 0] define_ring
244: f { . dehomogenize homogenize } map /f set
245: [f] groebner 0 get {dehomogenize} map /gg set
246:
247: /wv vdweight weightv def
248: gg { wv init } map /rrank.init set %%obtained the initial ideal
249: rrank.init {toString} map /chv set
250: /arg1 [chv v] rank def
251: ] pop
252: popEnv
253: popVariables
254: arg1
255: } def
256:
257:
258: [(rrank)
259: [( a rrank b)
260: ( array a; number b)
261: (It computes the holonomic rank for regular holonomic system.)
262: (For the algorithm, see Grobner deformations of hypergeometric differential equations, 1999, Springer.)
263: (Chapter 2.)
264: (Example 1 : )
265: $ [ [( (x Dx)^2 + ( y Dy)^2) ( x Dx y Dy -1)] (x,y)] rrank :: $
266: ]
267: ] putUsages
268: (rrank ) messagen-quiet
269:
270: /gb.v 1 def
271: /gb.verbose 0 def
272: /gb {
273: /arg1 set
274: [/in-gb /aa /typev /setarg /f /v
275: /gg /wv /termorder /vec /ans /rr /mm
276: ] pushVariables
277: [(CurrentRingp) (KanGBmessage)] pushEnv
278: [
279:
280: /aa arg1 def
281: aa isArray { } { ( << array >> gb) error } ifelse
282: /setarg 0 def
283: /wv 0 def
284: aa { tag } map /typev set
285: typev [ ArrayP ] eq
286: { /f aa 0 get def
287: /v gb.v def
288: /setarg 1 def
289: } { } ifelse
290: typev [ArrayP StringP] eq
291: { /f aa 0 get def
292: /v aa 1 get def
293: /setarg 1 def
294: } { } ifelse
295: typev [ArrayP ArrayP] eq
296: { /f aa 0 get def
297: /v aa 1 get from_records def
298: /setarg 1 def
299: } { } ifelse
300: typev [ArrayP StringP ArrayP] eq
301: { /f aa 0 get def
302: /v aa 1 get def
303: /wv aa 2 get def
304: /setarg 1 def
305: } { } ifelse
306: typev [ArrayP ArrayP ArrayP] eq
307: { /f aa 0 get def
308: /v aa 1 get from_records def
309: /wv aa 2 get def
310: /setarg 1 def
311: } { } ifelse
312:
313: setarg { } { (gb : Argument mismatch) error } ifelse
314:
315: [(KanGBmessage) gb.verbose ] system_variable
316:
317: %%% Start of the preprocess
318: f getRing /rr set
319: %% To the normal form : matrix expression.
320: f gb.toMatrixOfString /f set
321: /mm gb.itWasMatrix def
322:
323: rr tag 0 eq {
324: %% Define our own ring
325: v isInteger {
326: (Error in gb: Specify variables) error
327: } { } ifelse
328: wv isInteger {
329: [v ring_of_differential_operators
330: 0] define_ring
331: /termorder 1 def
332: }{
333: [v ring_of_differential_operators
334: wv weight_vector
335: 0] define_ring
336: wv gb.isTermOrder /termorder set
337: } ifelse
338: } {
339: %% Use the ring structre given by the input.
340: v isInteger not {
341: (Warning : the given ring definition is not used.) message
342: } { } ifelse
343: rr ring_def
344: /wv rr gb.getWeight def
345: wv gb.isTermOrder /termorder set
346: } ifelse
347: %%% Enf of the preprocess
348:
349:
350: termorder {
351: f { {. dehomogenize} map } map /f set
352: [f] groebner_sugar 0 get /gg set
353: }{
354: f { {. dehomogenize} map} map /f set
355: f fromVectors { homogenize } map /f set
356: [f] groebner 0 get /gg set
357: }ifelse
358: wv isInteger {
359: /ans [gg gg {init} map] def
360: }{
361: /ans [gg gg {wv 0 get weightv init} map] def
362: }ifelse
363:
364: %% Postprocess : recover the matrix expression.
365: mm {
366: ans { /tmp set [mm tmp] toVectors } map
367: /ans set
368: }{ }
369: ifelse
370: %%
371:
372: /arg1 ans def
373: ] pop
374: popEnv
375: popVariables
376: arg1
377: } def
378: (gb ) messagen-quiet
379:
380: /pgb {
381: /arg1 set
382: [/in-pgb /aa /typev /setarg /f /v
383: /gg /wv /termorder /vec /ans /rr /mm
384: ] pushVariables
385: [(CurrentRingp) (KanGBmessage) (UseCriterion1)] pushEnv
386: [
387:
388: /aa arg1 def
389: aa isArray { } { (<< array >> pgb) error } ifelse
390: /setarg 0 def
391: /wv 0 def
392: aa { tag } map /typev set
393: typev [ ArrayP ] eq
394: { /f aa 0 get def
395: /v gb.v def
396: /setarg 1 def
397: } { } ifelse
398: typev [ArrayP StringP] eq
399: { /f aa 0 get def
400: /v aa 1 get def
401: /setarg 1 def
402: } { } ifelse
403: typev [ArrayP ArrayP] eq
404: { /f aa 0 get def
405: /v aa 1 get from_records def
406: /setarg 1 def
407: } { } ifelse
408: typev [ArrayP StringP ArrayP] eq
409: { /f aa 0 get def
410: /v aa 1 get def
411: /wv aa 2 get def
412: /setarg 1 def
413: } { } ifelse
414: typev [ArrayP ArrayP ArrayP] eq
415: { /f aa 0 get def
416: /v aa 1 get from_records def
417: /wv aa 2 get def
418: /setarg 1 def
419: } { } ifelse
420:
421: setarg { } { (pgb : Argument mismatch) error } ifelse
422:
423: [(KanGBmessage) gb.verbose ] system_variable
424:
425: %%% Start of the preprocess
426: f getRing /rr set
427: %% To the normal form : matrix expression.
428: f gb.toMatrixOfString /f set
429: /mm gb.itWasMatrix def
430:
431: rr tag 0 eq {
432: %% Define our own ring
433: v isInteger {
434: (Error in pgb: Specify variables) error
435: } { } ifelse
436: wv isInteger {
437: [v ring_of_polynomials
438: 0] define_ring
439: /termorder 1 def
440: }{
441: [v ring_of_polynomials
442: wv weight_vector
443: 0] define_ring
444: wv gb.isTermOrder /termorder set
445: } ifelse
446: } {
447: %% Use the ring structre given by the input.
448: v isInteger not {
449: (Warning : the given ring definition is not used.) message
450: } { } ifelse
451: rr ring_def
452: /wv rr gb.getWeight def
453: wv gb.isTermOrder /termorder set
454: } ifelse
455: %%% Enf of the preprocess
456:
457:
458: termorder {
459: f { {. dehomogenize} map } map /f set
460: [(UseCriterion1) 1] system_variable
461: [f] groebner_sugar 0 get /gg set
462: [(UseCriterion1) 0] system_variable
463: }{
464: f { {. dehomogenize} map} map /f set
465: f fromVectors { homogenize } map /f set
466: [(UseCriterion1) 1] system_variable
467: [f] groebner 0 get /gg set
468: [(UseCriterion1) 0] system_variable
469: }ifelse
470: wv isInteger {
471: /ans [gg gg {init} map] def
472: }{
473: /ans [gg gg {wv 0 get weightv init} map] def
474: }ifelse
475:
476: %% Postprocess : recover the matrix expression.
477: mm {
478: ans { /tmp set [mm tmp] toVectors } map
479: /ans set
480: }{ }
481: ifelse
482: %%
483:
484: /arg1 ans def
485: ] pop
486: popEnv
487: popVariables
488: arg1
489: } def
490:
491: /pgb.old {
492: /arg1 set
493: [/in-pgb /aa /typev /setarg /f /v
494: /gg /wv /termorder /vec /ans
495: ] pushVariables
496: [(CurrentRingp) (KanGBmessage) (UseCriterion1)] pushEnv
497: [
498:
499: /aa arg1 def
500: aa isArray { } { (array pgb) message (pgb) usage error } ifelse
501: /setarg 0 def
502: /wv 0 def
503: aa { tag } map /typev set
504: typev [ ArrayP ] eq
505: { /f aa 0 get def
506: /v gb.v def
507: /setarg 1 def
508: } { } ifelse
509: typev [ArrayP StringP] eq
510: { /f aa 0 get def
511: /v aa 1 get def
512: /setarg 1 def
513: } { } ifelse
514: typev [ArrayP ArrayP] eq
515: { /f aa 0 get def
516: /v aa 1 get from_records def
517: /setarg 1 def
518: } { } ifelse
519: typev [ArrayP StringP ArrayP] eq
520: { /f aa 0 get def
521: /v aa 1 get def
522: /wv aa 2 get def
523: /setarg 1 def
524: } { } ifelse
525: typev [ArrayP ArrayP ArrayP] eq
526: { /f aa 0 get def
527: /v aa 1 get from_records def
528: /wv aa 2 get def
529: /setarg 1 def
530: } { } ifelse
531:
532: setarg { } { (pgb : Argument mismatch) message error } ifelse
533:
534: [(KanGBmessage) gb.verbose ] system_variable
535:
536: %% Input must not be vectors.
537: f { toString } map /f set
538:
539: wv isInteger {
540: [v ring_of_polynomials
541: 0] define_ring
542: /termorder 1 def
543: }{
544: [v ring_of_polynomials
545: wv weight_vector
546: 0] define_ring
547: wv gb.isTermOrder /termorder set
548: } ifelse
549: termorder {
550: f { . dehomogenize } map /f set
551: [(UseCriterion1) 1] system_variable
552: [f] groebner_sugar 0 get /gg set
553: [(UseCriterion1) 0] system_variable
554: }{
555: f { . dehomogenize homogenize} map /f set
556: [(UseCriterion1) 1] system_variable
557: [f] groebner 0 get /gg set
558: [(UseCriterion1) 0] system_variable
559: }ifelse
560: wv isInteger {
561: /ans [gg gg {init} map] def
562: }{
563: /ans [gg gg {wv 0 get weightv init} map] def
564: }ifelse
565: /arg1 ans def
566: ] pop
567: popEnv
568: popVariables
569: arg1
570: } def
571: (pgb ) messagen-quiet
572:
573: /gb.toMatrixOfString {
574: /arg1 set
575: [/in-gb.toMatrixOfString /ff /aa /ans] pushVariables
576: [
577: /aa arg1 def
578: aa length 0 eq { /ans [ ] def /gb.toMatrixOfString.LLL goto }{ } ifelse
579: aa 0 get isArray {
580: /gb.itWasMatrix aa 0 get length def
581: }{
582: /gb.itWasMatrix 0 def
583: } ifelse
584: aa {
585: /ff set
586: ff isArray {
587: ff {toString} map /ff set
588: }{
589: [ff toString] /ff set
590: } ifelse
591: ff
592: } map /ans set
593: /gb.toMatrixOfString.LLL
594: /arg1 ans def
595: ] pop
596: popVariables
597: arg1
598: } def
599: [(gb.toMatrixOfString)
600: [(It translates given input into a matrix form which is a data structure)
601: (for computations of kernel, image, cokernel, etc.)
602: (gb.itWasMatrix is set to the length of the input vector.)
603: $Example 1: $
604: $ [ (x). (y).] gb.toMatrixOfString ==> [[(x)] [(y)]] $
605: $ gb.itWasMatrix is 0.$
606: $Example 2: $
607: $ [ [(x). (1).] [(y). (0).]] gb.toMatrixOfString ==> [ [(x) (1)] [(y) (0)]] $
608: $ gb.itWasMatrix is 2.$
609: ]] putUsages
610:
611: /gb.toMatrixOfPoly {
612: /arg1 set
613: [/in-gb.toMatrixOfPoly /ff /aa /ans] pushVariables
614: [
615: /aa arg1 def
616: aa length 0 eq { /ans [ ] def /gb.toMatrixOfPoly.LLL goto }{ } ifelse
617: aa 0 get isArray {
618: /gb.itWasMatrix aa 0 get length def
619: }{
620: /gb.itWasMatrix 0 def
621: } ifelse
622: aa {
623: /ff set
624: ff isArray {
625: }{
626: [ff] /ff set
627: } ifelse
628: ff
629: } map /ans set
630: /gb.toMatrixOfPoly.LLL
631: /arg1 ans def
632: ] pop
633: popVariables
634: arg1
635: } def
636: [(gb.toMatrixOfPoly)
637: [(It translates given input into a matrix form which is a data structure)
638: (for computations of kernel, image, cokernel, etc.)
639: (gb.itWasMatrix is set to the length of the input vector.)
640: $Example 1: $
641: $ [ (x). (y).] gb.toMatrixOfPoly ==> [[(x)] [(y)]] $
642: $ gb.itWasMatrix is 0.$
643: $Example 2: $
644: $ [ [(x). (1).] [(y). (0).]] gb.toMatrixOfPoly ==> [ [(x) (1)] [(y) (0)]] $
645: $ gb.itWasMatrix is 2.$
646: ]] putUsages
647:
648: /gb.getWeight {
649: /arg1 set
650: [/in-gb.getWeight /rr /ww /vv /ans /nn /ii] pushVariables
651: [(CurrentRingp)] pushEnv
652: [
653: /rr arg1 def
654: rr ring_def
655: getVariableNames /vv set
656: [(orderMatrix)] system_variable 0 get /ww set
657: /nn vv length 1 sub def
658: [0 1 nn {
659: /ii set
660: ww ii get 0 eq {
661: } {
662: vv ii get
663: ww ii get
664: } ifelse
665: } for
666: ] /ans set
667: /arg1 [ans] def
668: ] pop
669: popEnv
670: popVariables
671: arg1
672: } def
673: [(gb.getWeight)
674: [(ring gb.getWeight wv)
675: (It gets the weight vector field of the ring ring.)
676: ]] putUsages
677:
678:
679: /gb.isTermOrder {
680: /arg1 set
681: [/in-gb.isTermOrder /vv /ww /yes /i /j] pushVariables
682: [
683: /vv arg1 def
684: /yes 1 def
685: 0 1 vv length 1 sub {
686: /i set
687: /ww vv i get def
688: 0 1 ww length 1 sub {
689: /j set
690: ww j get isInteger {
691: ww j get 0 lt { /yes 0 def } { } ifelse
692: }{ } ifelse
693: }for
694: }for
695: /arg1 yes def
696: ] pop
697: popVariables
698: arg1
699: } def
700: [(gb)
701: [(a gb b)
702: (array a; array b;)
703: (b : [g ii]; array g; array in; g is a Grobner basis of f)
704: ( in the ring of differential operators.)
705: $ ii is the initial ideal in case of w is given or <<a>> belongs$
706: $ to a ring. In the other cases, it returns the initial monominal.$
707: (a : [f ]; array f; f is a set of generators of an ideal in a ring.)
708: (a : [f v]; array f; string v; v is the variables. )
709: (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
710: ( )
711: $Example 1: [ [( (x Dx)^2 + (y Dy)^2 -1) ( x y Dx Dy -1)] (x,y) $
712: $ [ [ (Dx) 1 ] ] ] gb pmat ; $
713: (Example 2: )
714: (To put h=1, type in, e.g., )
715: $ [ [(2 x Dx + 3 y Dy+6) (2 y Dx + 3 x^2 Dy)] (x,y) $
716: $ [[(x) -1 (Dx) 1 (y) -1 (Dy) 1]]] gb /gg set gg dehomogenize pmat ;$
717: ( )
718: $Example 3: [ [( (x Dx)^2 + (y Dy)^2 -1) ( x y Dx Dy -1)] (x,y) $
719: $ [ [ (Dx) 1 (Dy) 1] ] ] gb pmat ; $
720: ( )
721: $Example 4: [[ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] (x,y) $
722: $ [ [ (x) -1 (y) -1] ] ] gb pmat ; $
723: ( )
724: (cf. gb, groebner, groebner_sugar, syz. )
725: ]] putUsages
726:
727: [(pgb)
728: [(a pgb b)
729: (array a; array b;)
730: (b : [g ii]; array g; array in; g is a Grobner basis of f)
731: ( in the ring of polynomials.)
732: $ ii is the initial ideal in case of w is given or <<a>>belongs$
733: $ to a ring. In the other cases, it returns the initial monominal.$
734: (a : [f ]; array f; f is a set of generators of an ideal in a ring.)
735: (a : [f v]; array f; string v; v is the variables.)
736: (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
737: $Example 1: [(x,y) ring_of_polynomials 0] define_ring $
738: $ [ [(x^2+y^2-4). (x y -1).] ] pgb :: $
739: $Example 2: [ [(x^2+y^2) (x y)] (x,y) [ [(x) -1 (y) -1] ] ] pgb :: $
740: (cf. gb, groebner, groebner_sugar, syz. )
741: ]] putUsages
742:
743:
744: %/syz.v 1 def
745: /syz.v 1 def
746: /syz.verbose 0 def
747: /syz {
748: /arg1 set
749: [/in-syz /aa /typev /setarg /f /v
750: /gg /wv /termorder /vec /ans /ggall /vectorInput /vsize /gtmp /gtmp2
751: /rr /mm
752: ] pushVariables
753: [(CurrentRingp) (KanGBmessage)] pushEnv
754: [
755:
756: /aa arg1 def
757: aa isArray { } { (<< array >> syz) error } ifelse
758: /setarg 0 def
759: /wv 0 def
760: aa { tag } map /typev set
761: typev [ ArrayP ] eq
762: { /f aa 0 get def
763: /v syz.v def
764: /setarg 1 def
765: } { } ifelse
766: typev [ArrayP StringP] eq
767: { /f aa 0 get def
768: /v aa 1 get def
769: /setarg 1 def
770: } { } ifelse
771: typev [ArrayP ArrayP] eq
772: { /f aa 0 get def
773: /v aa 1 get from_records def
774: /setarg 1 def
775: } { } ifelse
776: typev [ArrayP StringP ArrayP] eq
777: { /f aa 0 get def
778: /v aa 1 get def
779: /wv aa 2 get def
780: /setarg 1 def
781: } { } ifelse
782: typev [ArrayP ArrayP ArrayP] eq
783: { /f aa 0 get def
784: /v aa 1 get from_records def
785: /wv aa 2 get def
786: /setarg 1 def
787: } { } ifelse
788:
789: setarg { } { (syz : Argument mismatch) error } ifelse
790:
791: [(KanGBmessage) syz.verbose ] system_variable
792:
793:
794:
795: %%% Start of the preprocess
796: f getRing /rr set
797: %% To the normal form : matrix expression.
798: f gb.toMatrixOfString /f set
799: /mm gb.itWasMatrix def
800: mm 0 gt {
801: /vectorInput 1 def
802: }{
803: /vectorInput 1 def
804: } ifelse
805:
806: rr tag 0 eq {
807: %% Define our own ring
808: v isInteger {
809: (Error in syz: Specify variables) error
810: } { } ifelse
811: wv isInteger {
812: [v ring_of_differential_operators
813: 0] define_ring
814: /termorder 1 def
815: }{
816: [v ring_of_differential_operators
817: wv weight_vector
818: 0] define_ring
819: wv gb.isTermOrder /termorder set
820: } ifelse
821: }{
822: %% Use the ring structre given by the input.
823: v isInteger not {
824: (Warning : the given ring definition is not used.) message
825: } { } ifelse
826: rr ring_def
827: /wv rr gb.getWeight def
828: wv gb.isTermOrder /termorder set
829: } ifelse
830: %%% Enf of the preprocess
831:
832: termorder {
833: f { {. dehomogenize} map } map /f set
834: [f [(needBack) (needSyz)]] groebner_sugar /ggall set
835: ggall 2 get /gg set
836: }{
837: f { {. dehomogenize } map homogenize } map /f set
838: [f [(needBack) (needSyz)]] groebner /ggall set
839: ggall 2 get /gg set
840: }ifelse
841: vectorInput {
842: /vsize f 0 get length def %% input vector size.
843: /gtmp ggall 0 get def
844: [vsize gtmp] toVectors /gtmp set
845: ggall 0 gtmp put
846: }{ } ifelse
847: /arg1 [gg dehomogenize ggall] def
848: ] pop
849: popEnv
850: popVariables
851: arg1
852: } def
853: (syz ) messagen-quiet
854:
855: [(syz)
856: [(a syz [b c])
857: (array a; array b; array c)
858: (b is a set of generators of the syzygies of f.)
859: (c = [gb, backward transformation, syzygy without dehomogenization].)
860: (See groebner.)
861: (a : [f ]; array f; f is a set of generators of an ideal in a ring.)
862: (a : [f v]; array f; string v; v is the variables.)
863: (a : [f v w]; array f; string v; array of array w; w is the weight matirx.)
864: $Example 1: [(x,y) ring_of_polynomials 0] define_ring $
865: $ [ [(x^2+y^2-4). (x y -1).] ] syz :: $
866: $Example 2: [ [(x^2+y^2) (x y)] (x,y) [ [(x) -1 (y) -1] ] ] syz :: $
867: $Example 3: [ [( (x Dx)^2 + (y Dy)^2 -1) ( x y Dx Dy -1)] (x,y) $
868: $ [ [ (Dx) 1 ] ] ] syz pmat ; $
869: $Example 4: [ [(2 x Dx + 3 y Dy+6) (2 y Dx + 3 x^2 Dy)] (x,y) $
870: $ [[(x) -1 (Dx) 1 (y) -1 (Dy) 1]]] syz pmat ;$
871: $Example 5: [ [ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] $
872: $ (x,y) ] syz pmat ;$
873: $Example 6: [ [ [(x^2) (y+x)] [(x+y) (y^3)] [(2 x^2+x y) (y+x+x y^3)]] $
874: $ (x,y) [[(x) -1 (y) -2]] ] syz pmat ;$
875: $Example 7: [ [ [(0) (0)] [(0) (0)] [(x) (y)]] $
876: $ [(x) (y)]] syz pmat ;$
877: ]] putUsages
878:
879:
880: %%%%%%%%%%%%%%%%%% package fs %%%%%%%%%%%%%%%%%%%%%%%
881: [(genericAnn)
882: [ (f [s v1 v2 ... vn] genericAnn [L1 ... Lm])
883: (L1, ..., Lm are annihilating ideal for f^s.)
884: (f is a polynomial of v1, ..., vn)
885: (<string> | <poly> f, s, v1, ..., vn ; <poly> L1, ..., Lm )
886: $Example: (x^3+y^3+z^3) [(s) (x) (y) (z)] genericAnn$
887: ]
888: ] putUsages ( genericAnn ) messagen-quiet
889: /fs.verbose 0 def
890: /genericAnn {
891: /arg2 set /arg1 set
892: [/in-genericAnn /f /vlist /s /vvv /nnn /rrr
893: /v1 /ops /ggg /ggg0
894: ] pushVariables
895: [(CurrentRingp) (KanGBmessage)] pushEnv
896: [
897: /f arg1 def /vlist arg2 def
898: f toString /f set
899: vlist { toString } map /vlist set
900: [(KanGBmessage) fs.verbose] system_variable
901: /s vlist 0 get def
902: /vvv (_u,_v,_t,) vlist rest { (,) 2 cat_n } map aload length /nnn set
903: s nnn 2 add cat_n def
904: fs.verbose { vvv message } { }ifelse
905: [vvv ring_of_differential_operators
906: [[(_u) 1 (_v) 1]] weight_vector 0] define_ring /rrr set
907:
908: [ (_u*_t). f . sub (_u*_v-1). ]
909: vlist rest { /v1 set
910: %%D-clean f . (D) v1 2 cat_n . 1 diff0 (_v*D_t). mul
911: f . @@@.Dsymbol v1 2 cat_n . 1 diff0 [(_v*) @@@.Dsymbol (_t)] cat . mul
912: @@@.Dsymbol v1 2 cat_n . add } map
913: join
914: /ops set
915: ops {[[(h). (1).]] replace } map /ops set
916: fs.verbose { ops message } { }ifelse
917: [ops] groebner_sugar 0 get /ggg0 set
918: fs.verbose { ggg0 message } { } ifelse
919: ggg0 [(_u) (_v)] eliminatev
920: %%D-clean { [(_t).] [ (D_t).] [s .] distraction
921: { [(_t).] [ [@@@.Dsymbol (_t)] cat .] [s .] distraction
922: [[s . << (0). s . sub (1). sub >>]] replace
923: } map /arg1 set
924: ] pop
925: popEnv
926: popVariables
927: arg1
928: } def
929:
930: %% Find differential equations for f^(m), r0 the minimal integral root.
931: [(annfs)
932: [( [ f v m r0] annfs g )
933: (It returns the annihilating ideal of f^m where r0 must be smaller)
934: (or equal to the minimal integral root of the b-function.)
935: (Or, it returns the annihilating ideal of f^r0, r0 and the b-function)
936: (where r0 is the minial integral root of b.)
937: (For the algorithm, see J. Pure and Applied Algebra 117&118(1997), 495--518.)
938: (Example 1: [(x^2+y^2+z^2+t^2) (x,y,z,t) -1 -2] annfs :: )
939: $ It returns the annihilating ideal of (x^2+y^2+z^2+t^2)^(-1).$
940: (Example 2: [(x^2+y^2+z^2+t^2) (x,y,z,t)] annfs :: )
941: $ It returns the annihilating ideal of f^r0 and [r0, b-function]$
942: $ where r0 is the minimal integral root of the b-function.$
943: (Example 3: [(x^2+y^2+z^2) (x,y,z) -1 -1] annfs :: )
944: (Example 4: [(x^3+y^3+z^3) (x,y,z)] annfs :: )
945: (Example 5: [((x1+x2+x3)(x1 x2 + x2 x3 + x1 x3) - t x1 x2 x3 ) )
946: ( (t,x1,x2,x3) -1 -2] annfs :: )
947: ( Note that the example 4 uses huge memory space.)
948: ]] putUsages
949: ( annfs ) messagen-quiet
950: /annfs.verbose fs.verbose def
951: /annfs.v [(x) (y) (z)] def
952: /annfs.s (_s) def
953: %% The first variable must be s.
954: /annfs {
955: /arg1 set
956: [/in-annfs /aa /typev /setarg /v /m /r0 /gg /ss /fs /gg2
957: /ans /vtmp /w2 /velim /bbb /rrr /r0
958: ] pushVariables
959: [(CurrentRingp) (KanGBmessage)] pushEnv
960: [
961:
962: /aa arg1 def
963: aa isArray { } { ( << array >> annfs) error } ifelse
964: /setarg 0 def
965: aa { tag } map /typev set
966: /r0 [ ] def
967: /m [ ] def
968: /v annfs.v def
969: aa 0 << aa 0 get toString >> put
970: typev [ StringP ] eq
971: { /f aa 0 get def
972: /setarg 1 def
973: } { } ifelse
974: typev [StringP StringP] eq
975: { /f aa 0 get def
976: /v [ aa 1 get to_records pop ] def
977: /setarg 1 def
978: } { } ifelse
979: typev [StringP ArrayP] eq
980: { /f aa 0 get def
981: /v aa 1 get def
982: /setarg 1 def
983: } { } ifelse
984: typev [StringP ArrayP IntegerP IntegerP] eq
985: { /f aa 0 get def
986: /v aa 1 get def
987: /m aa 2 get def
988: /r0 aa 3 get def
989: /setarg 1 def
990: } { } ifelse
991: typev [StringP StringP IntegerP IntegerP] eq
992: { /f aa 0 get def
993: /v [ aa 1 get to_records pop ] def
994: /m aa 2 get def
995: /r0 aa 3 get def
996: /setarg 1 def
997: } { } ifelse
998: setarg 1 eq { } { (annfs : wrong argument) error } ifelse
999:
1000: [annfs.s] v join /v set
1001:
1002: /ss v 0 get def
1003: annfs.verbose {
1004: (f, v, s, f^{m}, m+r0 = ) messagen
1005: [ f (, ) v (, ) ss (, )
1006: (f^) m (,) m (+) r0 ] {messagen} map ( ) message
1007: } { } ifelse
1008:
1009: f v genericAnn /fs set
1010:
1011: annfs.verbose {
1012: (genericAnn is ) messagen fs message
1013: } { } ifelse
1014: [(KanGBmessage) annfs.verbose] system_variable
1015:
1016: m isArray {
1017: %% Now, let us find the b-function. /vtmp /w2 /velim /bbb /rrr /r0
1018: v rest { /vtmp set vtmp @@@.Dsymbol vtmp 2 cat_n } map /velim set
1019: velim { 1 } map /w2 set
1020: annfs.verbose { w2 message } { } ifelse
1021: [v from_records ring_of_differential_operators
1022: [w2] weight_vector 0] define_ring
1023: [ fs { toString . } map [ f toString . ] join ]
1024: groebner_sugar 0 get velim eliminatev 0 get /bbb set
1025: [[(s) annfs.s] from_records ring_of_polynomials 0] define_ring
1026: bbb toString . [[annfs.s . (s).]] replace /bbb set
1027: annfs.verbose { bbb message } { } ifelse
1028: bbb findIntegralRoots /rrr set
1029: rrr 0 get /r0 set %% minimal integral root.
1030: annfs.verbose { rrr message } { } ifelse
1031: fs 0 get (ring) dc ring_def
1032: fs { [[annfs.s . r0 toString .]] replace } map /ans set
1033: /ans [ans [r0 bbb]] def
1034: /annfs.label1 goto
1035: } { } ifelse
1036: m 0 ge {
1037: (annfs works only for getting annihilating ideal for f^(negative))
1038: error
1039: } { } ifelse
1040: r0 isArray {
1041: [(Need to compute the minimal root of b-function) nl
1042: (It has not been implemented.) ] cat
1043: error
1044: } { } ifelse
1045:
1046: [v from_records ring_of_differential_operators 0] define_ring
1047: fs {toString . dehomogenize [[ss . r0 (poly) dc]] replace}
1048: map /gg set
1049: annfs.verbose { gg message } { } ifelse
1050:
1051: [ [f . << m r0 sub >> npower ] gg join
1052: [(needBack) (needSyz)]] groebner_sugar 2 get /gg2 set
1053:
1054: gg2 { 0 get } map /ans set
1055: /ans ans { dup (0). eq {pop} { } ifelse } map def
1056:
1057: /annfs.label1
1058: /arg1 ans def
1059: ] pop
1060: popEnv
1061: popVariables
1062: arg1
1063: } def
1064:
1065: /genericAnnWithL.s (s) def
1066: /annfs.verify 0 def
1067: /genericAnnWithL {
1068: /arg1 set
1069: [/in-genericAnnWithL /aa /typev /setarg /v /m /r0 /gg /ss /fs /gg2
1070: /ans /vtmp /w2 /velim /bbb /rrr /r0 /myL /mygb /jj
1071: ] pushVariables
1072: [(CurrentRingp) (KanGBmessage) (Homogenize)] pushEnv
1073: [
1074:
1075: /aa arg1 def
1076: aa isArray { } { ( << array >> annfs) error } ifelse
1077: /setarg 0 def
1078: aa { tag } map /typev set
1079: /r0 [ ] def
1080: /m [ ] def
1081: /v annfs.v def
1082: aa 0 << aa 0 get toString >> put
1083: typev [ StringP ] eq
1084: { /f aa 0 get def
1085: /setarg 1 def
1086: } { } ifelse
1087: typev [StringP StringP] eq
1088: { /f aa 0 get def
1089: /v [ aa 1 get to_records pop ] def
1090: /setarg 1 def
1091: } { } ifelse
1092: typev [StringP ArrayP] eq
1093: { /f aa 0 get def
1094: /v aa 1 get def
1095: /setarg 1 def
1096: } { } ifelse
1097: setarg 1 eq { } { (genericAnnWithL : wrong argument) error } ifelse
1098:
1099: [genericAnnWithL.s] v join /v set
1100:
1101: /ss v 0 get def
1102: annfs.verbose {
1103: (f, v, s, f^{m}, m+r0 = ) messagen
1104: [ f (, ) v (, ) ss (, )
1105: (f^) m (,) m (+) r0 ] {messagen} map ( ) message
1106: } { } ifelse
1107:
1108: f v genericAnn /fs set
1109:
1110: annfs.verbose {
1111: (genericAnn is ) messagen fs message
1112: } { } ifelse
1113: [(KanGBmessage) annfs.verbose] system_variable
1114:
1115: m isArray {
1116: %% Now, let us find the b-function. /vtmp /w2 /velim /bbb /rrr /r0
1117: v rest { /vtmp set vtmp @@@.Dsymbol vtmp 2 cat_n } map /velim set
1118: velim { 1 } map /w2 set
1119: annfs.verbose { w2 message } { } ifelse
1120: [v from_records ring_of_differential_operators
1121: [w2] weight_vector 0] define_ring
1122:
1123: [ [ f toString . ] fs { toString . } map join [(needBack)]]
1124: groebner_sugar /mygb set
1125: mygb 0 get velim eliminatev 0 get /bbb set
1126: mygb 0 get bbb position /jj set
1127: mygb 1 get jj get 0 get /myL set
1128:
1129: annfs.verbose { bbb message } { } ifelse
1130:
1131: annfs.verify {
1132: (Verifying L f - b belongs to genericAnn(f)) message
1133: [(Homogenize) 0] system_variable
1134: << myL f . mul bbb sub >>
1135: [fs { toString . } map] groebner_sugar 0 get
1136: reduction 0 get message
1137: (Is it zero? Then it's fine.) message
1138: } { } ifelse
1139:
1140: /ans [bbb [myL fs] ] def
1141: /annfs.label1 goto
1142: } { } ifelse
1143:
1144: /annfs.label1
1145: /arg1 ans def
1146: ] pop
1147: popEnv
1148: popVariables
1149: arg1
1150: } def
1151:
1152:
1153: [(genericAnnWithL)
1154: [$[f v] genericAnnWithL [b [L I]]$
1155: $String f,v; poly b,L; array of poly I;$
1156: $f is a polynomial given by a string. v is the variables.$
1157: $ v must not contain names s, e.$
1158: $b is the b-function (Bernstein-Sato polynomial) for f and$
1159: $ L is the operator satisfying L f^{s+1} = b(s) f^s $
1160: $ I is the annihilating ideal of f^s.$
1161: $cf. bfunction, annfs, genericAnn.$
1162: $Example 1: [(x^2+y^2) (x,y)] genericAnnWithL ::$
1163: $Example 2: [(x^2+y^2+z^2) (x,y,z)] genericAnnWithL ::$
1164: $Example 3: [(x^3-y^2 z^2) (x,y,z)] genericAnnWithL ::$
1165: ]] putUsages
1166:
1167:
1168: ( ) message-quiet ;
1169:
1170:
1171:
1172:
1173:
1174:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>