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