1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: center.lib,v 1.13 2005-02-23 18:10:44 levandov Exp $"; |
---|
3 | category="Noncommutative"; |
---|
4 | info=" |
---|
5 | LIBRARY: center.lib computation of central elements of G-algebras and Q-algebras. |
---|
6 | AUTHOR: Oleksandr Motsak, motsak@mathematik.uni-kl.de. |
---|
7 | OVERVIEW: |
---|
8 | This is a library for computing the central elements and centralizers of elements in various noncommutative algebras. |
---|
9 | Implementation is based on algorithms, written in the frame of the diploma thesis by O. Motsak (advisor: Prof. S.A. Ovsienko, support: V. Levandovskyy), at Kyiv Taras Shevchenko University (Ukraine) with the title 'An algorithm for the computation of the center of noncommutative polynomial algebra'. |
---|
10 | |
---|
11 | SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' |
---|
12 | |
---|
13 | PROCEDURES: |
---|
14 | center(MaxDeg[,N]); computes the generators of the center of a basering, |
---|
15 | centralizer(f, MaxDeg[,N]); computes the generators of the centralizer of f in a basering, |
---|
16 | inCenter(l); checks the centrality of elements of list/ideal/poly l |
---|
17 | inCentralizer(l, f); checks the commutativity wrt polynomial f of polynomials of list/ideal/poly l |
---|
18 | |
---|
19 | KEYWORDS: inCenter; inCentralizer; center; centralizer |
---|
20 | "; |
---|
21 | |
---|
22 | |
---|
23 | /******************************************************/ |
---|
24 | // stuff |
---|
25 | |
---|
26 | /******************************************************/ |
---|
27 | static proc myValue ( def s, list # ) |
---|
28 | " |
---|
29 | return s or (typeof(s))(#) |
---|
30 | " |
---|
31 | { |
---|
32 | def @p = s; |
---|
33 | if ( size(#) > 0 ) |
---|
34 | { |
---|
35 | if ( typeof( #[1] ) == typeof(s) ) |
---|
36 | { |
---|
37 | @p = #[1]; |
---|
38 | }; |
---|
39 | }; |
---|
40 | return (@p); |
---|
41 | }; |
---|
42 | |
---|
43 | /******************************************************/ |
---|
44 | static proc myInt ( list # ) |
---|
45 | " |
---|
46 | return 0 or int(#) |
---|
47 | " |
---|
48 | { |
---|
49 | int @p = 0; |
---|
50 | return (myValue ( @p, # )); |
---|
51 | }; |
---|
52 | |
---|
53 | /******************************************************/ |
---|
54 | static proc myPoly ( list # ) |
---|
55 | " |
---|
56 | return 0 or poly(#) |
---|
57 | " |
---|
58 | { |
---|
59 | poly @p = 0; |
---|
60 | return (myValue( @p, # )); |
---|
61 | }; |
---|
62 | |
---|
63 | /******************************************************/ |
---|
64 | static proc myRing ( list # ) |
---|
65 | " |
---|
66 | return basring or ring(#) |
---|
67 | " |
---|
68 | { |
---|
69 | def @p = basering; |
---|
70 | return (myValue ( @p, # )); |
---|
71 | }; |
---|
72 | |
---|
73 | /******************************************************/ |
---|
74 | static proc myString ( list # ) |
---|
75 | " |
---|
76 | return basring or ring(#) |
---|
77 | " |
---|
78 | { |
---|
79 | string @p = ""; |
---|
80 | return (myValue ( @p, # )); |
---|
81 | }; |
---|
82 | |
---|
83 | /******************************************************/ |
---|
84 | static proc myIdeal ( list # ) |
---|
85 | " |
---|
86 | return 0 or ideal(#) |
---|
87 | " |
---|
88 | { |
---|
89 | ideal @p; |
---|
90 | return (myValue ( @p, # )); |
---|
91 | }; |
---|
92 | |
---|
93 | /******************************************************/ |
---|
94 | // for debug |
---|
95 | |
---|
96 | /******************************************************/ |
---|
97 | static proc toprint() |
---|
98 | { |
---|
99 | return (0); |
---|
100 | }; |
---|
101 | |
---|
102 | /******************************************************/ |
---|
103 | static proc Print( list # ) |
---|
104 | { |
---|
105 | if ( toprint() ) |
---|
106 | { |
---|
107 | print (#); |
---|
108 | }; |
---|
109 | }; |
---|
110 | |
---|
111 | /******************************************************/ |
---|
112 | static proc BCall(string Name, list #) |
---|
113 | { |
---|
114 | if( toprint() ) |
---|
115 | { |
---|
116 | "CALL: ", Name, "( ", string(#), " )"; |
---|
117 | }; |
---|
118 | }; |
---|
119 | |
---|
120 | /******************************************************/ |
---|
121 | static proc ECall(string Name, list #) |
---|
122 | { |
---|
123 | if( toprint() ) |
---|
124 | { |
---|
125 | "RET : ", Name, "(...)={ ", string(#), " }"; |
---|
126 | }; |
---|
127 | }; |
---|
128 | |
---|
129 | |
---|
130 | /******************************************************/ |
---|
131 | // "my" degres functions and lie procedures... |
---|
132 | |
---|
133 | /******************************************************/ |
---|
134 | static proc maxDeg( def z ) |
---|
135 | " |
---|
136 | returns: int : max of givendeg( z_i ), among all z_i \in z |
---|
137 | returns -1 if z is empty or contain only 0 polynomial |
---|
138 | " |
---|
139 | { |
---|
140 | int max = -1; |
---|
141 | int d; |
---|
142 | |
---|
143 | for(int i=size(z); i>0; i-- ) |
---|
144 | { |
---|
145 | d = deg(z[i]); |
---|
146 | if( d > max ) |
---|
147 | { |
---|
148 | max = d; |
---|
149 | }; |
---|
150 | }; |
---|
151 | |
---|
152 | return(max); |
---|
153 | }; |
---|
154 | |
---|
155 | /******************************************************/ |
---|
156 | // some business for my_var |
---|
157 | |
---|
158 | /******************************************************/ |
---|
159 | static proc myCoeff ( poly f, poly m ) |
---|
160 | " |
---|
161 | input: poly f, |
---|
162 | return: coeff at m |
---|
163 | " |
---|
164 | { |
---|
165 | m = leadmonom(m); |
---|
166 | |
---|
167 | int i = size(f); |
---|
168 | |
---|
169 | while ( (i>0) and (leadmonom(f[i])<m) ) |
---|
170 | { |
---|
171 | i--; |
---|
172 | }; |
---|
173 | if( i == 0 ) |
---|
174 | { |
---|
175 | return ( 0 ); |
---|
176 | }; |
---|
177 | |
---|
178 | if(leadmonom(f[i]) == m) |
---|
179 | { |
---|
180 | return ( leadcoef(f[i]) ); |
---|
181 | }; |
---|
182 | |
---|
183 | return ( 0 ); |
---|
184 | }; |
---|
185 | |
---|
186 | /******************************************************/ |
---|
187 | static proc my_vars () |
---|
188 | { |
---|
189 | int N = nvars( basering ); |
---|
190 | list V = list(); |
---|
191 | |
---|
192 | for ( int k = 1; k <= N; k++ ) |
---|
193 | { |
---|
194 | V[k] = var(k); |
---|
195 | }; |
---|
196 | |
---|
197 | return (V); |
---|
198 | }; |
---|
199 | |
---|
200 | /******************************************************/ |
---|
201 | static proc my_commutators () |
---|
202 | { |
---|
203 | list V = my_vars(); |
---|
204 | int N = size( V ); |
---|
205 | |
---|
206 | matrix M[N][N] = 0; |
---|
207 | |
---|
208 | poly v, d; int i; |
---|
209 | |
---|
210 | for ( int k = 2; k <= N; k++ ) |
---|
211 | { |
---|
212 | v = V[k]; |
---|
213 | for ( i = 1; i < k; i++ ) |
---|
214 | { |
---|
215 | d = V[i]; |
---|
216 | M[k,i] = v*d - d*v; // [var(k),var(i)] |
---|
217 | M[i,k] = -M[k,i]; // [var(i),var(k)] == -[var(k),var(i)] |
---|
218 | }; |
---|
219 | }; |
---|
220 | |
---|
221 | return (M); |
---|
222 | }; |
---|
223 | |
---|
224 | /******************************************************/ |
---|
225 | static proc my_associated () |
---|
226 | " |
---|
227 | return list of records: |
---|
228 | [1] - poly: variable, |
---|
229 | [2] - bool: isCartan |
---|
230 | [3] - matrix: associated matrix (basering should be compatible) |
---|
231 | " |
---|
232 | { |
---|
233 | int i, j; poly a; poly d; |
---|
234 | |
---|
235 | int N = nvars( basering ); |
---|
236 | list V = my_vars(); |
---|
237 | matrix M = my_commutators(); |
---|
238 | |
---|
239 | matrix A[N][N]; |
---|
240 | |
---|
241 | int cartan; |
---|
242 | |
---|
243 | list RESULT = list(); |
---|
244 | int r_begin = 1; |
---|
245 | int r_end = N; |
---|
246 | |
---|
247 | |
---|
248 | for ( int k = 1; k <= N; k++ ) |
---|
249 | { |
---|
250 | A = 0; |
---|
251 | cartan = 1; |
---|
252 | for ( i = 1; i <= N; i++ ) |
---|
253 | { |
---|
254 | d = M[k,i]; |
---|
255 | for ( j = 1; j <= N; j++ ) |
---|
256 | { |
---|
257 | a = myCoeff( d, V[j] ); |
---|
258 | A[i,j] = a; |
---|
259 | |
---|
260 | if( (i!=j) and (a!=0) ) |
---|
261 | { |
---|
262 | cartan = 0; |
---|
263 | }; |
---|
264 | }; |
---|
265 | }; |
---|
266 | |
---|
267 | if ( cartan ) |
---|
268 | { |
---|
269 | RESULT[r_begin] = list( V[k], cartan, A ); |
---|
270 | r_begin++; |
---|
271 | } else |
---|
272 | { |
---|
273 | RESULT[r_end] = list( V[k], cartan, A ); |
---|
274 | r_end--; |
---|
275 | }; |
---|
276 | |
---|
277 | }; |
---|
278 | |
---|
279 | return (RESULT); |
---|
280 | |
---|
281 | }; |
---|
282 | |
---|
283 | /******************************************************/ |
---|
284 | static proc my_var_init() |
---|
285 | " |
---|
286 | must be called once before using my_var!!! |
---|
287 | " |
---|
288 | { |
---|
289 | if( defined( @@@SORTEDVARARRAY ) ) |
---|
290 | { |
---|
291 | kill ( @@@SORTEDVARARRAY ); |
---|
292 | }; |
---|
293 | |
---|
294 | list @@@SORTEDVARARRAY; |
---|
295 | list V = my_associated(); |
---|
296 | |
---|
297 | // cartans - first then others... |
---|
298 | |
---|
299 | for ( int i = 1; i<= size(V); i++ ) |
---|
300 | { |
---|
301 | @@@SORTEDVARARRAY[i] = V[i][1]; |
---|
302 | }; |
---|
303 | |
---|
304 | Print( "@@@SORTEDVARARRAY: " + string(@@@SORTEDVARARRAY) ); |
---|
305 | |
---|
306 | export(@@@SORTEDVARARRAY); |
---|
307 | }; |
---|
308 | |
---|
309 | /******************************************************/ |
---|
310 | static proc my_var(int @number ) |
---|
311 | " |
---|
312 | 'my' var (with my order of variables), |
---|
313 | before use call my_var_init() |
---|
314 | " |
---|
315 | { |
---|
316 | if( defined( @@@SORTEDVARARRAY ) ) |
---|
317 | { |
---|
318 | return( @@@SORTEDVARARRAY[@number] ); |
---|
319 | }; |
---|
320 | |
---|
321 | Print( "Error: my_var_init() was not called before this..." ); |
---|
322 | return( var(@number) ); |
---|
323 | |
---|
324 | }; |
---|
325 | |
---|
326 | /******************************************************/ |
---|
327 | static proc my_var_done() |
---|
328 | " |
---|
329 | this procedure should be called after finishing using my_var in order to clean up. |
---|
330 | " |
---|
331 | { |
---|
332 | if( defined( @@@SORTEDVARARRAY ) ) |
---|
333 | { |
---|
334 | kill ( @@@SORTEDVARARRAY ); |
---|
335 | }; |
---|
336 | }; |
---|
337 | |
---|
338 | |
---|
339 | /******************************************************/ |
---|
340 | // QNF stuff |
---|
341 | |
---|
342 | /******************************************************/ |
---|
343 | static proc myBitParam( int param, int n ) |
---|
344 | " |
---|
345 | return: string: param in bin within no less then n digits (adding zeroes) |
---|
346 | " |
---|
347 | { |
---|
348 | string s = ""; |
---|
349 | |
---|
350 | while ( param != 0 ) |
---|
351 | { |
---|
352 | s = string ( param % 2 ) + s; |
---|
353 | param = param / 2; |
---|
354 | n --; |
---|
355 | }; |
---|
356 | while ( n > 0 ) |
---|
357 | { |
---|
358 | s = "0" + s; |
---|
359 | n --; |
---|
360 | }; |
---|
361 | return ( s ); |
---|
362 | }; |
---|
363 | |
---|
364 | /******************************************************/ |
---|
365 | static proc QNF_init(list #) |
---|
366 | { |
---|
367 | if( defined(@@@MY_QIDEAL) ) |
---|
368 | { |
---|
369 | kill (@@@MY_QIDEAL); |
---|
370 | }; |
---|
371 | |
---|
372 | if( typeof(basering) == "qring" ) |
---|
373 | { |
---|
374 | ideal @@@MY_QIDEAL = twostd(myIdeal(#)); |
---|
375 | // setup redSB, redTail options |
---|
376 | export(@@@MY_QIDEAL); |
---|
377 | |
---|
378 | option(redSB); |
---|
379 | option(redTail); |
---|
380 | |
---|
381 | }; |
---|
382 | }; |
---|
383 | |
---|
384 | /******************************************************/ |
---|
385 | static proc QNF_done() |
---|
386 | { |
---|
387 | if( defined(@@@MY_QIDEAL) ) |
---|
388 | { |
---|
389 | kill (@@@MY_QIDEAL); |
---|
390 | }; |
---|
391 | }; |
---|
392 | |
---|
393 | /******************************************************/ |
---|
394 | static proc QNF_poly ( poly p, list # ) |
---|
395 | { |
---|
396 | if( defined(@@@MY_QIDEAL) ) |
---|
397 | { |
---|
398 | /* |
---|
399 | int param = 1 ; // by def: reduce both 1st and 2nd entries |
---|
400 | param = myValue(param, #); |
---|
401 | string bits = myBitParam( param, 8 ); |
---|
402 | */ |
---|
403 | BCall( "QNF_poly", p ); |
---|
404 | |
---|
405 | p = NF( p, @@@MY_QIDEAL ); |
---|
406 | |
---|
407 | ECall( "QNF_poly", p ); |
---|
408 | }; // QRing |
---|
409 | |
---|
410 | return(p); |
---|
411 | }; |
---|
412 | |
---|
413 | |
---|
414 | /******************************************************/ |
---|
415 | static proc QNF_list ( list l, list # ) |
---|
416 | " |
---|
417 | expects list of records with two polynomial entries |
---|
418 | and in case of Qring allpies NF( * , std(0) ) to entries (depends on parameter) |
---|
419 | |
---|
420 | 1234 5678 |
---|
421 | |
---|
422 | +1 8 => reduce every 1st entry |
---|
423 | +2 7 => reduce every 2nd entry |
---|
424 | |
---|
425 | --- (for PBW in qrings) --- |
---|
426 | |
---|
427 | +4 6 => kill out pbw entries where pbw monom (1) was affected |
---|
428 | |
---|
429 | //?? wofuer? |
---|
430 | +8 5 => reduce every 3rd entry |
---|
431 | |
---|
432 | " |
---|
433 | { |
---|
434 | |
---|
435 | if( defined(@@@MY_QIDEAL) ) |
---|
436 | { |
---|
437 | |
---|
438 | int param = 1 + 2; // by def: reduce both 1st and 2nd entries |
---|
439 | param = myValue(param, #); |
---|
440 | |
---|
441 | string bits = myBitParam( param, 8 ); |
---|
442 | |
---|
443 | BCall( "QNF_list", "list, " + bits ); |
---|
444 | |
---|
445 | poly temp; |
---|
446 | |
---|
447 | for ( int i = size(l); i>0 ; i -- ) |
---|
448 | { |
---|
449 | |
---|
450 | if ( typeof( l[i] ) == "poly" ) |
---|
451 | { |
---|
452 | |
---|
453 | if ( (bits[8] == "1") or (bits[6] == "1") ) |
---|
454 | { |
---|
455 | temp = NF( l[i], @@@MY_QIDEAL ); |
---|
456 | |
---|
457 | if ( bits[6] == "1" ) |
---|
458 | {// for PBW in qrings: kill out pbw entries where pbw monom was reduced |
---|
459 | |
---|
460 | if( temp != l[i] ) |
---|
461 | { |
---|
462 | l = delete( l, i ); |
---|
463 | i --; |
---|
464 | continue; |
---|
465 | }; |
---|
466 | }; |
---|
467 | |
---|
468 | if ( bits[8] == "1" ) |
---|
469 | { |
---|
470 | l[i] = temp; |
---|
471 | }; |
---|
472 | }; |
---|
473 | }; |
---|
474 | |
---|
475 | if ( typeof( l[i] ) == "list" ) |
---|
476 | { |
---|
477 | // 1st |
---|
478 | |
---|
479 | if ( size(l[i])>0 ) |
---|
480 | { |
---|
481 | if ( (bits[8] == "1") or (bits[6] == "1") ) |
---|
482 | { |
---|
483 | if( typeof( l[i][1] ) == "poly" ) |
---|
484 | { |
---|
485 | |
---|
486 | temp = NF( l[i][1], @@@MY_QIDEAL ); |
---|
487 | |
---|
488 | if ( bits[6] == "1" ) |
---|
489 | {// for PBW in qrings: kill out pbw entries where pbw monom was reduced |
---|
490 | |
---|
491 | if( temp != l[i][1] ) |
---|
492 | { |
---|
493 | l = delete( l, i ); |
---|
494 | i --; |
---|
495 | continue; |
---|
496 | }; |
---|
497 | }; |
---|
498 | |
---|
499 | if ( bits[8] == "1" ) |
---|
500 | { |
---|
501 | l[i][1] = temp; |
---|
502 | }; |
---|
503 | |
---|
504 | }; |
---|
505 | }; |
---|
506 | }; |
---|
507 | |
---|
508 | // 2nd |
---|
509 | |
---|
510 | if ( size(l[i])>1 ) |
---|
511 | { |
---|
512 | if ( bits[7] == "1" ) |
---|
513 | { |
---|
514 | if( typeof( l[i][2] ) == "poly" ) |
---|
515 | { |
---|
516 | temp = NF( l[i][2], @@@MY_QIDEAL ); |
---|
517 | |
---|
518 | l[i][2] = temp; |
---|
519 | }; |
---|
520 | }; |
---|
521 | }; |
---|
522 | }; |
---|
523 | }; |
---|
524 | |
---|
525 | ECall( "QNF_list", "list" ); |
---|
526 | |
---|
527 | }; // Qring |
---|
528 | |
---|
529 | return ( l ); |
---|
530 | }; |
---|
531 | |
---|
532 | /******************************************************/ |
---|
533 | // things for zCommon |
---|
534 | |
---|
535 | /******************************************************/ |
---|
536 | static proc index( int i ) |
---|
537 | " |
---|
538 | i need indexing of arrays from 0 => ARRAY[ index(i) ] - everywhere... |
---|
539 | " |
---|
540 | { |
---|
541 | return (i + 1); |
---|
542 | }; |
---|
543 | |
---|
544 | /******************************************************/ |
---|
545 | static proc uni_poly( poly p ) |
---|
546 | " |
---|
547 | returns polynomial with the same monomials but without coefficients. |
---|
548 | " |
---|
549 | { |
---|
550 | poly @tt = poly(0); |
---|
551 | for ( int @k = size(p); @k > 0; @k -- ) |
---|
552 | { |
---|
553 | @tt = @tt + leadmonom(p[@k]); |
---|
554 | }; |
---|
555 | return (@tt); |
---|
556 | }; |
---|
557 | |
---|
558 | /******************************************************/ |
---|
559 | static proc get_max ( poly @t, number @def ) |
---|
560 | { |
---|
561 | int @n = size( @t ); |
---|
562 | |
---|
563 | if ( @n == 0 ) |
---|
564 | { |
---|
565 | return (@def); |
---|
566 | }; |
---|
567 | |
---|
568 | number @max = leadcoef(@t[1]); |
---|
569 | number @mm; |
---|
570 | |
---|
571 | if ( @n > 1) |
---|
572 | { |
---|
573 | for ( int @i = 2; @i <= @n ;@i ++ ) |
---|
574 | { |
---|
575 | @mm = leadcoef ( @t[@i] ); |
---|
576 | if ( @mm < 0 ) |
---|
577 | { |
---|
578 | @mm = -@mm; |
---|
579 | }; |
---|
580 | |
---|
581 | if( @mm > @max ) |
---|
582 | { |
---|
583 | @max = @mm; |
---|
584 | }; |
---|
585 | }; |
---|
586 | }; |
---|
587 | |
---|
588 | @max = @max + 1; |
---|
589 | if ( @max == 0 ) |
---|
590 | { |
---|
591 | @max = @max + 1; |
---|
592 | }; |
---|
593 | return ( 1/@max ); |
---|
594 | }; |
---|
595 | |
---|
596 | |
---|
597 | /******************************************************/ |
---|
598 | static proc get_uni_sum(list @given, int @num) |
---|
599 | { |
---|
600 | int @l, @k; |
---|
601 | |
---|
602 | poly @t, @tt; |
---|
603 | |
---|
604 | @l = size (@given); |
---|
605 | |
---|
606 | @t = poly(0); |
---|
607 | for ( @k = @l; @k > 0; @k -- ) |
---|
608 | { |
---|
609 | if (@num == 1) |
---|
610 | { |
---|
611 | @tt = @given[@k]; |
---|
612 | } else |
---|
613 | { |
---|
614 | @tt = @given[@k][2]; |
---|
615 | }; |
---|
616 | |
---|
617 | @t = @t + uni_poly( @tt ); |
---|
618 | }; |
---|
619 | |
---|
620 | return ( uni_poly(@t) ); |
---|
621 | }; |
---|
622 | |
---|
623 | |
---|
624 | |
---|
625 | /******************************************************/ |
---|
626 | static proc LM ( intvec exp ) |
---|
627 | " |
---|
628 | input : given exponent |
---|
629 | return: monom with this exponent... |
---|
630 | " |
---|
631 | { |
---|
632 | poly @f = 1; int @deg; |
---|
633 | for( int @i = 1; @i <= size(exp); @i ++ ) |
---|
634 | { |
---|
635 | @deg = exp[@i]; |
---|
636 | if ( @deg > 0 ) |
---|
637 | { |
---|
638 | @f = @f * (var(@i)^(@deg)); |
---|
639 | }; |
---|
640 | }; |
---|
641 | |
---|
642 | return (@f); |
---|
643 | }; |
---|
644 | |
---|
645 | /******************************************************/ |
---|
646 | // reduced centers - computation of "minimal" set of generators |
---|
647 | |
---|
648 | LIB "general.lib"; // for sort |
---|
649 | |
---|
650 | /******************************************************/ |
---|
651 | static proc zSort ( list @z ) |
---|
652 | " |
---|
653 | Sort elements of a list of polynoms, |
---|
654 | and normalize |
---|
655 | " |
---|
656 | { |
---|
657 | int n = size(@z); |
---|
658 | if ( n < 1 ) |
---|
659 | {// empty list |
---|
660 | return( @z ); |
---|
661 | }; |
---|
662 | |
---|
663 | if ( n == 1 ) |
---|
664 | { |
---|
665 | if ( @z[1] == 0 ) // if zero => empty list |
---|
666 | { |
---|
667 | return(list()); |
---|
668 | }; |
---|
669 | |
---|
670 | @z[1] = @z[1] * ( 1/leadcoef(@z[1]) ) ; |
---|
671 | |
---|
672 | return( @z ); |
---|
673 | }; |
---|
674 | |
---|
675 | int i = 1; |
---|
676 | |
---|
677 | while ( i<=n ) |
---|
678 | { |
---|
679 | if (size( @z[i] ) != 0) |
---|
680 | { |
---|
681 | break; |
---|
682 | }; |
---|
683 | i++; |
---|
684 | }; |
---|
685 | |
---|
686 | if ( i > n ) |
---|
687 | { // all zeroes |
---|
688 | return(list()); |
---|
689 | }; |
---|
690 | |
---|
691 | ideal id = @z[i]; |
---|
692 | i++; |
---|
693 | |
---|
694 | while ( i<= n ) |
---|
695 | { |
---|
696 | if( @z[i] != 0 ) |
---|
697 | { |
---|
698 | id=@z[i],id; |
---|
699 | }; |
---|
700 | i++; |
---|
701 | }; |
---|
702 | |
---|
703 | // can sort only ideal generators:-( |
---|
704 | list srt = sort(id); |
---|
705 | |
---|
706 | poly p; |
---|
707 | // convert srt[1]:ideal to result:list. |
---|
708 | list result = list(); |
---|
709 | for ( i = size(srt[1]); i>0; i-- ) |
---|
710 | { |
---|
711 | p = srt[1][i]; |
---|
712 | if ( p == 0 ) |
---|
713 | { |
---|
714 | i --; |
---|
715 | continue; |
---|
716 | }; |
---|
717 | p = p* (1/leadcoef(p)); // normalize |
---|
718 | result = list(p) + result; |
---|
719 | }; |
---|
720 | |
---|
721 | // "OUT SORT::"; |
---|
722 | // result; |
---|
723 | |
---|
724 | return ( result ); |
---|
725 | }; |
---|
726 | |
---|
727 | |
---|
728 | |
---|
729 | /******************************************************/ |
---|
730 | static proc zRefine ( list @z ) |
---|
731 | " |
---|
732 | kill terms by leading monomials... |
---|
733 | Note: based on myCoeff |
---|
734 | " |
---|
735 | { |
---|
736 | BCall("zRefine", "list"); |
---|
737 | |
---|
738 | int @i, @j; |
---|
739 | |
---|
740 | poly @f = 0; // leadmonom |
---|
741 | |
---|
742 | poly @ff = 0; // polyes |
---|
743 | poly @gg = 0; |
---|
744 | number @nf, @ng; |
---|
745 | |
---|
746 | int flag = 1; |
---|
747 | |
---|
748 | while ( flag == 1 ) |
---|
749 | { |
---|
750 | flag = 0; |
---|
751 | |
---|
752 | @z = zSort(@z); // sort out, < ... |
---|
753 | |
---|
754 | if( size(@z) < 2 ) |
---|
755 | { |
---|
756 | return (@z); |
---|
757 | }; |
---|
758 | |
---|
759 | for ( @i = size(@z); @i > 0; @i -- ) // at 1st the biggest ... then smallest... |
---|
760 | { |
---|
761 | |
---|
762 | @ff = @z[@i]; |
---|
763 | |
---|
764 | if( size(@ff) == 0 ) |
---|
765 | { |
---|
766 | @z = delete ( @z , @i ); |
---|
767 | @i --; |
---|
768 | continue; |
---|
769 | }; |
---|
770 | |
---|
771 | @ff = @ff*(1/leadcoef(@ff)); |
---|
772 | @z[@i] = @ff; |
---|
773 | @f = leadmonom(@ff); |
---|
774 | |
---|
775 | for ( @j = (@i-1); (@j>0); @j -- ) |
---|
776 | { |
---|
777 | @gg = @z[@j]; |
---|
778 | @ng = myCoeff(@gg, @f); // leads? |
---|
779 | if( @ng!=0 ) |
---|
780 | { |
---|
781 | @z[@j] = @gg - @ng * @ff; |
---|
782 | flag = 1; |
---|
783 | }; |
---|
784 | }; |
---|
785 | for ( @j = (@i+1); (@j<=size(@z)); @j ++ ) |
---|
786 | { |
---|
787 | @gg = @z[@j]; |
---|
788 | @ng = myCoeff(@gg, @f); |
---|
789 | if( @ng!=0 ) |
---|
790 | { |
---|
791 | @z[@j] = @gg - @ng * @ff; |
---|
792 | flag = 1; |
---|
793 | }; |
---|
794 | }; |
---|
795 | |
---|
796 | }; |
---|
797 | }; |
---|
798 | |
---|
799 | ECall("zRefine", "list"); |
---|
800 | |
---|
801 | return ( @z ); |
---|
802 | |
---|
803 | }; |
---|
804 | |
---|
805 | |
---|
806 | /******************************************************/ |
---|
807 | // procedures for building "bad leadmonomials" set |
---|
808 | |
---|
809 | |
---|
810 | /******************************************************/ |
---|
811 | static proc checkPolyUniq( list l, poly p ) |
---|
812 | " |
---|
813 | check whether p sits already in l, assume l to be size-sorted < |
---|
814 | return: -1 if present |
---|
815 | 1 if we need to add |
---|
816 | " |
---|
817 | { |
---|
818 | BCall( "checkPolyUniq", "{ " + string(l) + " }, " + string(p) ); |
---|
819 | int n = size(l); |
---|
820 | int i = 1; |
---|
821 | |
---|
822 | int s = size(p); |
---|
823 | |
---|
824 | while( i<= n ) |
---|
825 | { |
---|
826 | if ( size(l[i]) >= s ) |
---|
827 | { |
---|
828 | break; |
---|
829 | }; |
---|
830 | |
---|
831 | i ++; |
---|
832 | }; |
---|
833 | |
---|
834 | // now: size(l[i]) >= s |
---|
835 | while( i<= n ) |
---|
836 | { |
---|
837 | if ( size(l[i]) == s ) |
---|
838 | { |
---|
839 | break; |
---|
840 | |
---|
841 | }; |
---|
842 | if ( l[i] == p ) |
---|
843 | { |
---|
844 | ECall( "checkPolyUniq", -1 ); |
---|
845 | return (-1); |
---|
846 | }; |
---|
847 | i ++; |
---|
848 | }; |
---|
849 | |
---|
850 | ECall( "checkPolyUniq", 1 ); |
---|
851 | return ( 1 ); |
---|
852 | |
---|
853 | }; |
---|
854 | |
---|
855 | /******************************************************/ |
---|
856 | static proc addPolyUniq( list l, poly p ) |
---|
857 | " |
---|
858 | add p into l uniquely, and keep l size-sorted < |
---|
859 | " |
---|
860 | { |
---|
861 | BCall( "addPolyUniq", " { " + string(l) + " }, " + string(p) ); |
---|
862 | |
---|
863 | int n = size(l); |
---|
864 | |
---|
865 | if ( n == 0 ) |
---|
866 | { |
---|
867 | l = list(p); |
---|
868 | |
---|
869 | ECall( "addPolyUniq", l ); |
---|
870 | |
---|
871 | return (l); |
---|
872 | }; |
---|
873 | |
---|
874 | int s = size(p); |
---|
875 | |
---|
876 | int i = 1; |
---|
877 | while( i<= n ) |
---|
878 | { |
---|
879 | if( size(l[i]) > s ) |
---|
880 | { |
---|
881 | l = insert( l, p, i-1 ) ; |
---|
882 | break; |
---|
883 | }; |
---|
884 | |
---|
885 | if( size(l[i]) == s ) |
---|
886 | { |
---|
887 | if ( l[i] == p ) |
---|
888 | { |
---|
889 | break; |
---|
890 | }; |
---|
891 | }; |
---|
892 | |
---|
893 | i++; |
---|
894 | }; |
---|
895 | |
---|
896 | if( i > n ) |
---|
897 | { |
---|
898 | l = l + list(p); |
---|
899 | }; |
---|
900 | |
---|
901 | ECall( "addPolyUniq", l ); |
---|
902 | return(l); |
---|
903 | }; |
---|
904 | |
---|
905 | |
---|
906 | /******************************************************/ |
---|
907 | static proc mergePolysUniq( list a, list b ) |
---|
908 | " |
---|
909 | merge lists uniq |
---|
910 | " |
---|
911 | { |
---|
912 | BCall( "mergePolysUniq", "{ " + string(a) + " }, { " + string(b) + "} " ); |
---|
913 | |
---|
914 | for( int i = 1; i <= size(b); i ++ ) |
---|
915 | { |
---|
916 | a = addPolyUniq(a, b[i]); |
---|
917 | }; |
---|
918 | |
---|
919 | ECall( "mergePolysUniq", a ); |
---|
920 | |
---|
921 | return (a); |
---|
922 | }; |
---|
923 | |
---|
924 | |
---|
925 | /******************************************************/ |
---|
926 | static proc sortPolysUniq( list a ) |
---|
927 | " |
---|
928 | sort list uniq |
---|
929 | " |
---|
930 | { |
---|
931 | BCall( "sortPolysUniq", a ); |
---|
932 | |
---|
933 | if( size(a) < 2 ) |
---|
934 | { |
---|
935 | ECall( "sortPolysUniq", a ); |
---|
936 | return(a); |
---|
937 | }; |
---|
938 | |
---|
939 | list b = list(a[1]); |
---|
940 | |
---|
941 | for( int i = 2; i <= size(a); i ++ ) |
---|
942 | { |
---|
943 | b = addPolyUniq(b, a[i]); |
---|
944 | }; |
---|
945 | |
---|
946 | ECall( "sortPolysUniq", b ); |
---|
947 | |
---|
948 | return (b); |
---|
949 | }; |
---|
950 | |
---|
951 | /******************************************************/ |
---|
952 | static proc addRecordUniq ( list leadD, list newD, intvec texp, poly tm, int kind, list prs ) |
---|
953 | " |
---|
954 | if kind = 0 => for PBW - no products |
---|
955 | if kind = 1 => with products |
---|
956 | " |
---|
957 | { |
---|
958 | BCall ("addRecordUniq", "{old leads}, {new lead}, [" + string(texp) + "], " + string(tm) + ", " + string(kind) + ", {" + string(prs) + "}"); |
---|
959 | |
---|
960 | int i; |
---|
961 | int f_add = 0; |
---|
962 | |
---|
963 | prs = sortPolysUniq(prs); |
---|
964 | |
---|
965 | // trick: |
---|
966 | // check for presens of a monomial @tm in current index poly of @leads (=> in list @leads) |
---|
967 | // if char = 0 then new_size > (if not present) or = (if already present) |
---|
968 | // if char > 0 then new_size > (if not present) or =< (if already present) |
---|
969 | // !!!!! |
---|
970 | if( size(tm + leadD[2]) > size(leadD[2]) ) |
---|
971 | { |
---|
972 | f_add = 1; |
---|
973 | } else |
---|
974 | { |
---|
975 | if ( kind != 0 ) |
---|
976 | { |
---|
977 | for ( i = 1; i<= size(leadD[1]); i++ ) |
---|
978 | { |
---|
979 | if ( leadD[1][i][2] == tm ) |
---|
980 | { |
---|
981 | for ( i = size(prs); i>0; i-- ) |
---|
982 | { |
---|
983 | f_add = checkPolyUniq( leadD[1][i][3], prs[i] ); |
---|
984 | if( f_add == -1 ) |
---|
985 | { |
---|
986 | prs = delete(prs, i); |
---|
987 | }; |
---|
988 | }; |
---|
989 | |
---|
990 | break; |
---|
991 | }; |
---|
992 | }; |
---|
993 | }; |
---|
994 | }; |
---|
995 | |
---|
996 | if ( f_add == 1 ) |
---|
997 | { |
---|
998 | |
---|
999 | list newlist ; |
---|
1000 | if ( kind != 0 ) |
---|
1001 | { |
---|
1002 | newlist = list ( list ( texp, tm, prs ) ); |
---|
1003 | } else |
---|
1004 | { |
---|
1005 | newlist = list ( list ( texp, tm ) ); |
---|
1006 | }; |
---|
1007 | |
---|
1008 | |
---|
1009 | if ( size(newD[1]) == 0 ) |
---|
1010 | { |
---|
1011 | newD[1] = newlist; |
---|
1012 | newD[2] = tm; |
---|
1013 | } else |
---|
1014 | { |
---|
1015 | |
---|
1016 | if( size(tm + newD[2]) > size(newD[2]) ) |
---|
1017 | { |
---|
1018 | newD[1] = newD[1] + newlist; |
---|
1019 | newD[2] = newD[2] + tm; |
---|
1020 | } else |
---|
1021 | { |
---|
1022 | if ( kind != 0 ) |
---|
1023 | { |
---|
1024 | for ( i = 1; i<= size(newD[1]); i++ ) |
---|
1025 | { |
---|
1026 | if ( newD[1][i][2] == tm ) |
---|
1027 | { |
---|
1028 | newD[1][i][3] = mergePolysUniq( newD[1][i][3], prs ); |
---|
1029 | break; |
---|
1030 | }; |
---|
1031 | }; |
---|
1032 | }; |
---|
1033 | }; |
---|
1034 | }; |
---|
1035 | }; |
---|
1036 | |
---|
1037 | ECall("addRecordUniq", "{new leads}"); |
---|
1038 | return (newD); |
---|
1039 | }; |
---|
1040 | |
---|
1041 | |
---|
1042 | /******************************************************/ |
---|
1043 | static proc mergeRecordsUniq ( list old, list new, int kind ) |
---|
1044 | { |
---|
1045 | BCall ("mergeRecordsUniq", "{old leads}, {new leads}, " + string(kind) ); |
---|
1046 | |
---|
1047 | if( size(new[1]) > 0 ) |
---|
1048 | { |
---|
1049 | if ( size (old[1]) == 0 ) |
---|
1050 | { |
---|
1051 | old = new; |
---|
1052 | } else |
---|
1053 | { |
---|
1054 | if ( kind != 0 ) |
---|
1055 | { |
---|
1056 | int io; |
---|
1057 | for ( int in = 1; in <= size( new[1] ); in ++ ) |
---|
1058 | { |
---|
1059 | if( size( new[1][in][2] + old[2] ) > size( old[2] ) ) |
---|
1060 | { |
---|
1061 | old[1] = old[1] + list(new[1][in]); |
---|
1062 | old[2] = old[2] + new[1][in][2]; |
---|
1063 | } else |
---|
1064 | { |
---|
1065 | for( io = 1; io <= size( old[1] ); io ++ ) |
---|
1066 | { |
---|
1067 | if( old[1][io][2] == new[1][in][2] ) |
---|
1068 | { |
---|
1069 | old[1][io][3] = mergePolysUniq( old[1][io][3], new[1][in][3] ); |
---|
1070 | break; |
---|
1071 | }; |
---|
1072 | }; |
---|
1073 | }; |
---|
1074 | }; |
---|
1075 | } else |
---|
1076 | { |
---|
1077 | old[1] = old[1] + new[1]; |
---|
1078 | old[2] = old[2] + new[2]; |
---|
1079 | }; |
---|
1080 | }; |
---|
1081 | }; |
---|
1082 | |
---|
1083 | ECall ("mergeRecordsUniq", "{new leads}"); |
---|
1084 | return (old); |
---|
1085 | }; |
---|
1086 | |
---|
1087 | |
---|
1088 | |
---|
1089 | /******************************************************/ |
---|
1090 | static proc init_bads(int @deg) |
---|
1091 | " |
---|
1092 | initialization of an empty 'badleads' list |
---|
1093 | " |
---|
1094 | { |
---|
1095 | list @l = list(); |
---|
1096 | for( int @i = 0; @i <= @deg; @i ++ ) |
---|
1097 | { |
---|
1098 | @l[index(@i)] = list( list() , 0); |
---|
1099 | // new items: |
---|
1100 | // { |
---|
1101 | // list of leads, - empty list |
---|
1102 | // sum poly "index", - poly(0) |
---|
1103 | // } |
---|
1104 | }; |
---|
1105 | return (@l); |
---|
1106 | }; |
---|
1107 | |
---|
1108 | /******************************************************/ |
---|
1109 | static proc zAddBad ( list @leads, list @newz, int @maxDeg, int kind ) |
---|
1110 | " |
---|
1111 | input: |
---|
1112 | @leads: graded by deg list of |
---|
1113 | rec: |
---|
1114 | { |
---|
1115 | [1] - list of |
---|
1116 | rec: |
---|
1117 | { |
---|
1118 | [1] - leadexp. |
---|
1119 | [2] - leadmonom([1]) |
---|
1120 | if kind != 0 (for zReduce) => |
---|
1121 | [3] - !list! of all possible products which give this leadexp |
---|
1122 | }, |
---|
1123 | [2] - summ of all in [1] ('index' poly) |
---|
1124 | } |
---|
1125 | @newz: new elements for adding to @leads |
---|
1126 | @maxDeg: maximal degree |
---|
1127 | |
---|
1128 | if kind != 0 => keeps also list of all possible products which give leadexp of a record |
---|
1129 | |
---|
1130 | return: |
---|
1131 | updated @leads list |
---|
1132 | |
---|
1133 | " |
---|
1134 | { |
---|
1135 | BCall ("zAddBad", "{old leads}, { " + string(@newz) + " }, " + string(@maxDeg) + ", " + string(kind) + ""); |
---|
1136 | |
---|
1137 | int @newSize = size(@newz); |
---|
1138 | if ( @newSize < 1 ) |
---|
1139 | { |
---|
1140 | return (@leads); |
---|
1141 | }; |
---|
1142 | |
---|
1143 | int @i, @j, @k, @dd, @deg; |
---|
1144 | |
---|
1145 | |
---|
1146 | poly @m, @tm, @ttm; |
---|
1147 | intvec @exp, @texp, @ttexp; |
---|
1148 | |
---|
1149 | list @newzz; |
---|
1150 | int @size, @newdeg, @mydeg; |
---|
1151 | |
---|
1152 | poly @sum_old, @sum_new; |
---|
1153 | |
---|
1154 | poly a, b, @temp; |
---|
1155 | |
---|
1156 | /* |
---|
1157 | if kind = 0 => for PBW - no products |
---|
1158 | if kind = 1 => for zReduce - with products |
---|
1159 | */ |
---|
1160 | |
---|
1161 | poly pr = 0; int pri; list prs; |
---|
1162 | |
---|
1163 | for ( @i = @newSize; @i > 0; @i -- ) |
---|
1164 | {// for every new poly (@newz[@i]) do |
---|
1165 | |
---|
1166 | @m = leadmonom( @newz[@i] ); |
---|
1167 | @deg = deg(@m); |
---|
1168 | @exp = leadexp(@m); |
---|
1169 | |
---|
1170 | // newzz void new list |
---|
1171 | @newzz = init_bads(@maxDeg ); |
---|
1172 | |
---|
1173 | for( @mydeg = @deg; @mydeg <= @maxDeg; @mydeg = @mydeg + @deg ) |
---|
1174 | {// adding all possiblities for @newz[@i]^@j; |
---|
1175 | |
---|
1176 | if ( @mydeg > @deg ) |
---|
1177 | { |
---|
1178 | @texp = @texp + @exp; |
---|
1179 | @tm = LM ( @texp ); |
---|
1180 | if ( kind != 0) |
---|
1181 | { |
---|
1182 | pr = QNF_poly( pr * @newz[@i] ); // degrees must be there!!! |
---|
1183 | }; |
---|
1184 | } else |
---|
1185 | { |
---|
1186 | @texp = @exp; |
---|
1187 | @tm = @m; |
---|
1188 | if ( kind != 0) |
---|
1189 | { |
---|
1190 | pr = @newz[@i]; |
---|
1191 | }; |
---|
1192 | }; |
---|
1193 | |
---|
1194 | @temp = QNF_poly(@tm); |
---|
1195 | if( @temp != @tm ) |
---|
1196 | { |
---|
1197 | break; |
---|
1198 | }; |
---|
1199 | |
---|
1200 | /*!!*/ @newzz[index(@mydeg)] = |
---|
1201 | /*!!*/ addRecordUniq( @leads[index(@mydeg)], @newzz[index(@mydeg)], @texp, @tm, kind, list(pr) ); |
---|
1202 | |
---|
1203 | for ( @dd = 1; (@dd <= @maxDeg) and ((@dd + @mydeg) <= @maxDeg ); @dd ++ ) |
---|
1204 | { // for every good "deg" |
---|
1205 | |
---|
1206 | @newdeg = @dd + @mydeg; // any deg should be additive!!! |
---|
1207 | |
---|
1208 | for ( @k = size(@leads[index(@dd)][1]); @k > 0; @k -- ) |
---|
1209 | { |
---|
1210 | |
---|
1211 | @ttexp = (@leads[index(@dd)][1][@k][1]) + @texp; |
---|
1212 | @ttm = LM (@ttexp); |
---|
1213 | |
---|
1214 | if ( kind != 0 ) |
---|
1215 | { |
---|
1216 | prs = list(); |
---|
1217 | |
---|
1218 | for( pri = 1; pri <= size(@leads[index(@dd)][1][@k][3]); pri++) |
---|
1219 | { |
---|
1220 | // to do products into list and add one list !!! |
---|
1221 | a = QNF_poly( pr*@leads[index(@dd)][1][@k][3][pri]); |
---|
1222 | b = QNF_poly( @leads[index(@dd)][1][@k][3][pri]*pr); |
---|
1223 | |
---|
1224 | prs= prs + list(a); |
---|
1225 | |
---|
1226 | if ( a!=b ) |
---|
1227 | { |
---|
1228 | prs= prs + list(b); |
---|
1229 | }; |
---|
1230 | }; |
---|
1231 | |
---|
1232 | } else |
---|
1233 | { |
---|
1234 | prs = list(pr); |
---|
1235 | }; |
---|
1236 | |
---|
1237 | @temp = QNF_poly(@ttm); |
---|
1238 | if( @temp == @ttm ) |
---|
1239 | { |
---|
1240 | |
---|
1241 | /*!!*/ @newzz[index(@newdeg)] = |
---|
1242 | /*!!*/ addRecordUniq( @leads[index(@newdeg)], @newzz[index(@newdeg)], @ttexp, @ttm, kind, prs ); |
---|
1243 | }; |
---|
1244 | |
---|
1245 | |
---|
1246 | }; |
---|
1247 | }; // for |
---|
1248 | |
---|
1249 | if ( @deg == 0 ) |
---|
1250 | { |
---|
1251 | break; |
---|
1252 | }; |
---|
1253 | }; |
---|
1254 | |
---|
1255 | for ( @mydeg = 0; @mydeg <= @maxDeg; @mydeg ++ ) |
---|
1256 | { // adding currently generated to result |
---|
1257 | @leads[index(@mydeg)] = mergeRecordsUniq ( @leads[index(@mydeg)], @newzz[index(@mydeg)], kind ); |
---|
1258 | }; |
---|
1259 | |
---|
1260 | }; |
---|
1261 | |
---|
1262 | ECall ("zAddBad", "{new leads}"); |
---|
1263 | |
---|
1264 | return (@leads); |
---|
1265 | }; |
---|
1266 | |
---|
1267 | |
---|
1268 | /******************************************************/ |
---|
1269 | // procedure for reducing a given poly wrt already calculated "badleadings" |
---|
1270 | |
---|
1271 | /******************************************************/ |
---|
1272 | static proc zReducePoly ( list leads, poly new, poly anfang ) |
---|
1273 | " |
---|
1274 | reduce poly new wrt found leads, |
---|
1275 | return: list of all possible reductions... |
---|
1276 | " |
---|
1277 | { |
---|
1278 | poly temp = new; |
---|
1279 | poly rest = anfang; |
---|
1280 | |
---|
1281 | poly lm; |
---|
1282 | number lc; |
---|
1283 | |
---|
1284 | list LEADS; |
---|
1285 | |
---|
1286 | int i, n; |
---|
1287 | list prs; |
---|
1288 | |
---|
1289 | while ( size(temp) > 0 ) |
---|
1290 | { |
---|
1291 | // do for every term... beginning with leading... 'till smallest |
---|
1292 | lm = leadmonom( temp ); |
---|
1293 | |
---|
1294 | LEADS = leads[index(deg( lm ))]; // currently users bad leading products... |
---|
1295 | |
---|
1296 | if( size(LEADS[2] + lm ) <= size(LEADS[2]) ) |
---|
1297 | { // need to reduce, since leacmonom already in LEADS |
---|
1298 | |
---|
1299 | for ( i = 1; i <= size(LEADS[1]); i++ ) |
---|
1300 | { |
---|
1301 | if( LEADS[1][i][2] == lm ) |
---|
1302 | { |
---|
1303 | |
---|
1304 | lc = leadcoef( temp ); // no need be the unit |
---|
1305 | |
---|
1306 | prs = LEADS[1][i][3]; // shouldbe generated by zAddBad with kind == 1 |
---|
1307 | n = size(prs) ; |
---|
1308 | |
---|
1309 | if ( n == 1 ) |
---|
1310 | { // no recursion |
---|
1311 | |
---|
1312 | temp = leadcoef(prs[1]) * temp - lc * prs[1]; // leadmonom goes down |
---|
1313 | |
---|
1314 | } else |
---|
1315 | { // do recursion |
---|
1316 | |
---|
1317 | list result = list(); |
---|
1318 | poly newnew; |
---|
1319 | int f_addrest = 0; |
---|
1320 | |
---|
1321 | for( int pri = 1; pri <= n ; pri ++ ) |
---|
1322 | { |
---|
1323 | newnew = leadcoef(prs[pri]) * temp - lc * prs[pri]; // leadmonom goes down |
---|
1324 | |
---|
1325 | if( size( newnew ) > 0 ) |
---|
1326 | { |
---|
1327 | result = result + zReducePoly(leads, newnew, rest); |
---|
1328 | } else |
---|
1329 | { |
---|
1330 | f_addrest = 1; |
---|
1331 | }; |
---|
1332 | }; |
---|
1333 | |
---|
1334 | if ( f_addrest == 1 ) |
---|
1335 | { |
---|
1336 | result = result + list(rest); |
---|
1337 | }; |
---|
1338 | return ( result ); |
---|
1339 | |
---|
1340 | }; |
---|
1341 | break; |
---|
1342 | }; |
---|
1343 | }; |
---|
1344 | |
---|
1345 | } else |
---|
1346 | { // no such leadmonom in leads |
---|
1347 | |
---|
1348 | rest = rest + lead ( temp ); |
---|
1349 | temp = temp - lead ( temp ); // leadcoeff goes down |
---|
1350 | }; |
---|
1351 | }; |
---|
1352 | |
---|
1353 | return (list(rest)); |
---|
1354 | |
---|
1355 | }; |
---|
1356 | |
---|
1357 | |
---|
1358 | /******************************************************/ |
---|
1359 | static proc zCancel ( list @tt, list @leads ) |
---|
1360 | " |
---|
1361 | just kill entries of plane PBW base with leading monomials in @leads... |
---|
1362 | " |
---|
1363 | { |
---|
1364 | int @i; |
---|
1365 | |
---|
1366 | if ( (size(@tt) == 0) ) |
---|
1367 | { |
---|
1368 | return (@tt); |
---|
1369 | }; |
---|
1370 | |
---|
1371 | list result = list(); |
---|
1372 | poly g, f; |
---|
1373 | |
---|
1374 | for ( @i = size(@tt); @i > 0 ; @i -- ) |
---|
1375 | // for all PBW entries: |
---|
1376 | { |
---|
1377 | g = leadmonom(@tt[@i][1]); // pbw entry |
---|
1378 | f = @leads[index(deg(g))][2]; // 'index' poly from @leads |
---|
1379 | |
---|
1380 | if ( size(f + g) > size(f) ) // if g not in @leads (works only for monomials) |
---|
1381 | { |
---|
1382 | result = list( @tt[@i] ) + result; |
---|
1383 | }; |
---|
1384 | }; |
---|
1385 | |
---|
1386 | return (result); |
---|
1387 | |
---|
1388 | }; |
---|
1389 | |
---|
1390 | /******************************************************/ |
---|
1391 | // procedure for computing "minimal" set of generators... |
---|
1392 | |
---|
1393 | /******************************************************/ |
---|
1394 | static proc zReduce( list @z ) |
---|
1395 | " |
---|
1396 | reduce a set @z - base of Center as V.S. |
---|
1397 | into a minimal set!! |
---|
1398 | " |
---|
1399 | { |
---|
1400 | BCall( "zReduce", @z ); |
---|
1401 | |
---|
1402 | @z = zRefine ( @z ); |
---|
1403 | int n = size( @z ); |
---|
1404 | |
---|
1405 | if ( n < 2 ) |
---|
1406 | { |
---|
1407 | return (@z); |
---|
1408 | }; |
---|
1409 | |
---|
1410 | int d = maxDeg( @z ); |
---|
1411 | |
---|
1412 | if( d== -1 ) |
---|
1413 | { |
---|
1414 | return (@z); |
---|
1415 | }; |
---|
1416 | |
---|
1417 | list leads = init_bads( d ); |
---|
1418 | |
---|
1419 | list @red; |
---|
1420 | |
---|
1421 | int f_add; int j; |
---|
1422 | |
---|
1423 | list result = list(); |
---|
1424 | |
---|
1425 | poly p; |
---|
1426 | |
---|
1427 | for ( int i = 1; i <= n ; i++ ) |
---|
1428 | {// in this order... from least to maximal... |
---|
1429 | |
---|
1430 | p = @z[i]; |
---|
1431 | |
---|
1432 | @red = zReducePoly ( leads, p, 0 ); |
---|
1433 | |
---|
1434 | f_add = 1; |
---|
1435 | |
---|
1436 | for( j = 1; j <= size(@red); j++ ) |
---|
1437 | { |
---|
1438 | if ( @red[j] == 0 ) |
---|
1439 | { |
---|
1440 | f_add = 0; |
---|
1441 | break; // nothing new.... |
---|
1442 | }; |
---|
1443 | }; |
---|
1444 | |
---|
1445 | @red = zRefine( @red ); // there will be no zeroes after ??? |
---|
1446 | |
---|
1447 | if( size(@red) > 0 ) |
---|
1448 | { |
---|
1449 | leads = zAddBad( leads, @red, d, 1); |
---|
1450 | }; |
---|
1451 | |
---|
1452 | if ( f_add == 1 ) |
---|
1453 | { |
---|
1454 | // we got something new.... |
---|
1455 | result = result + list(@red); // ??? which one to add? => trying to add all |
---|
1456 | }; |
---|
1457 | }; |
---|
1458 | |
---|
1459 | ECall( "zReduce", result ); |
---|
1460 | |
---|
1461 | return (result); |
---|
1462 | }; |
---|
1463 | |
---|
1464 | |
---|
1465 | /******************************************************/ |
---|
1466 | // PBW procedures .... |
---|
1467 | |
---|
1468 | /* |
---|
1469 | PBW: graded by deg list of lists of PBW terms. |
---|
1470 | PBW term: list of 2 polynoms: [1] = PBW monom, [2] = Ad_p( [1] ), [3] = last variable in [1] (for speed) |
---|
1471 | */ |
---|
1472 | |
---|
1473 | /******************************************************/ |
---|
1474 | static proc PBW_init() |
---|
1475 | " |
---|
1476 | PBW[ index(0) ] = PBW part of degree 0: |
---|
1477 | record: |
---|
1478 | [1] : 1 // var(0) ;-) |
---|
1479 | [2] : 0 |
---|
1480 | [3] : 1 |
---|
1481 | " |
---|
1482 | { |
---|
1483 | BCall( "PBW_init" ); |
---|
1484 | return (list(list(1, 0, 1))); |
---|
1485 | }; |
---|
1486 | |
---|
1487 | |
---|
1488 | /******************************************************/ |
---|
1489 | static proc PBW_sys(list VARS, poly p) |
---|
1490 | " |
---|
1491 | calculate the array [1..NVARS] of records: |
---|
1492 | record[i] = |
---|
1493 | [1] = var(i) // i^{th} variable |
---|
1494 | [2] = p*[1]-[1]*p // [ p, var(i) ] |
---|
1495 | [3] = deg(var(i)) // i^{th} variable's deg |
---|
1496 | |
---|
1497 | ?? need debug ?? |
---|
1498 | " |
---|
1499 | { |
---|
1500 | |
---|
1501 | BCall( "PBW_sys" ); |
---|
1502 | |
---|
1503 | poly t; |
---|
1504 | for (int v = size(VARS); v>0; v -- ) |
---|
1505 | { |
---|
1506 | t = VARS[v]; |
---|
1507 | VARS[v] = list( t, QNF_poly( p*t - t*p ), deg(t) ) ; |
---|
1508 | }; |
---|
1509 | |
---|
1510 | return (VARS); |
---|
1511 | }; |
---|
1512 | |
---|
1513 | /******************************************************/ |
---|
1514 | static proc PBW_done( list PBW ) |
---|
1515 | " |
---|
1516 | collects all together, from graded lists into plane list. |
---|
1517 | |
---|
1518 | Note: also the last vars... |
---|
1519 | " |
---|
1520 | { |
---|
1521 | BCall( "PBW_done" ); |
---|
1522 | |
---|
1523 | list result = list(); |
---|
1524 | |
---|
1525 | int n = size(PBW); |
---|
1526 | for ( int i = 1; i <= n; i++) |
---|
1527 | { |
---|
1528 | result = result + PBW[i]; |
---|
1529 | }; |
---|
1530 | |
---|
1531 | return (result); |
---|
1532 | }; |
---|
1533 | |
---|
1534 | /******************************************************/ |
---|
1535 | static proc PBW_next ( list PBW, int k, list sys) |
---|
1536 | { |
---|
1537 | BCall( "PBW_next", k ); |
---|
1538 | |
---|
1539 | list temp; |
---|
1540 | int kk, nn, ii; |
---|
1541 | list result = list(); // PBW[ index(k) ] ought to be empty ?? |
---|
1542 | int N = size(sys); |
---|
1543 | |
---|
1544 | for ( int i = 1; i <= N; i ++ ) // for all vars |
---|
1545 | { |
---|
1546 | kk = (k - sys[i][3]); // any deg should be additive function wrt multiplication |
---|
1547 | if ( kk >= 0 ) |
---|
1548 | { |
---|
1549 | nn = size( PBW[ index(kk) ] ); |
---|
1550 | for ( ii = 1; ii <= nn; ii ++ ) |
---|
1551 | { |
---|
1552 | temp = PBW[ index(kk) ][ii]; |
---|
1553 | |
---|
1554 | if ( temp[3] <= i ) |
---|
1555 | { |
---|
1556 | temp[2] = temp[2]*sys[i][1] + temp[1]*sys[i][2]; // recursive [,] |
---|
1557 | temp[1] = temp[1]*sys[i][1]; |
---|
1558 | temp[3] = i; |
---|
1559 | |
---|
1560 | result = result + list ( temp ); |
---|
1561 | }; |
---|
1562 | }; |
---|
1563 | }; |
---|
1564 | }; |
---|
1565 | |
---|
1566 | result = QNF_list( result, 2 + 4); |
---|
1567 | ECall( "PBW_next", "list result" ); |
---|
1568 | return (result); |
---|
1569 | |
---|
1570 | }; |
---|
1571 | |
---|
1572 | /******************************************************/ |
---|
1573 | static proc PBW_base( int MaxDeg, poly p ) |
---|
1574 | { |
---|
1575 | BCall( "PBW_base", MaxDeg, p ); |
---|
1576 | |
---|
1577 | list PBW = list(); |
---|
1578 | list SYS = PBW_sys( my_vars(), p ); |
---|
1579 | |
---|
1580 | PBW[index(0)] = PBW_init(); |
---|
1581 | |
---|
1582 | for (int k = 1; k <= MaxDeg; k ++ ) |
---|
1583 | { |
---|
1584 | PBW[index(k)] = PBW_next( PBW, k, SYS ); |
---|
1585 | }; |
---|
1586 | |
---|
1587 | return (PBW_done( PBW )); |
---|
1588 | }; |
---|
1589 | |
---|
1590 | |
---|
1591 | /******************************************************/ |
---|
1592 | // standard center procedures... |
---|
1593 | |
---|
1594 | /******************************************************/ |
---|
1595 | static proc FSS (list @given, poly @BASE_POLY) |
---|
1596 | " |
---|
1597 | Gauss with computation of kernel v.s basis |
---|
1598 | to be optimized |
---|
1599 | " |
---|
1600 | { |
---|
1601 | BCall( "FSS", "list", "basepoly: ", @BASE_POLY); |
---|
1602 | |
---|
1603 | list @nums = list (); |
---|
1604 | intvec @ones; |
---|
1605 | int @j, @k, |
---|
1606 | @n, @v, |
---|
1607 | @a, @nn; |
---|
1608 | |
---|
1609 | @n = size( @BASE_POLY ); |
---|
1610 | @v = size( @given ); |
---|
1611 | |
---|
1612 | if ( (@v == 0) or (@n == 0) ) |
---|
1613 | { |
---|
1614 | return (@given); |
---|
1615 | }; |
---|
1616 | |
---|
1617 | matrix MD[@n][@v]; |
---|
1618 | |
---|
1619 | number @max; |
---|
1620 | poly @test; |
---|
1621 | poly @t; |
---|
1622 | |
---|
1623 | list LM = list(); // rec: { exp, poly } |
---|
1624 | |
---|
1625 | for( @k = @v; @k > 0; @k -- ) |
---|
1626 | { |
---|
1627 | LM[@k] = list(); |
---|
1628 | @t = @given[@k][2]; |
---|
1629 | LM[@k][2]= @t; |
---|
1630 | LM[@k][1]= leadexp( @t ); |
---|
1631 | }; |
---|
1632 | |
---|
1633 | |
---|
1634 | intvec @w; |
---|
1635 | for ( @a = 1 ; @a <= @n ; @a ++ ) |
---|
1636 | { |
---|
1637 | @w = leadexp(@BASE_POLY[@a]); |
---|
1638 | for( @k = @v; @k > 0; @k -- ) |
---|
1639 | { |
---|
1640 | if( @w == LM[@k][1] ) |
---|
1641 | { |
---|
1642 | @t = LM[@k][2]; |
---|
1643 | MD[@a, @k] = leadcoef( @t ); |
---|
1644 | @t = @t - lead ( @t ); |
---|
1645 | LM[@k][2] = @t; |
---|
1646 | |
---|
1647 | LM[@k][1] = leadexp(@t); |
---|
1648 | }; |
---|
1649 | }; |
---|
1650 | |
---|
1651 | }; |
---|
1652 | |
---|
1653 | int @x, @y; |
---|
1654 | number @min; |
---|
1655 | |
---|
1656 | number @div; |
---|
1657 | // @nums = list(); |
---|
1658 | |
---|
1659 | // Gauss |
---|
1660 | // print("Before Gauss, Matrix: "); |
---|
1661 | // print( MD ); |
---|
1662 | |
---|
1663 | |
---|
1664 | @x = 1; |
---|
1665 | @y = 1; |
---|
1666 | while ( (@y <= @n) && (@x <= @v)) |
---|
1667 | // main Gauss loop... |
---|
1668 | { |
---|
1669 | @min = leadcoef( MD[@y, @x] ); // curr diag. |
---|
1670 | |
---|
1671 | if ( @min == 0 ) |
---|
1672 | // if zero on diag... |
---|
1673 | { |
---|
1674 | @j = 0; |
---|
1675 | |
---|
1676 | // let's find the minimal |
---|
1677 | for ( @k = @y+1; @k <= @n; @k ++ ) |
---|
1678 | { |
---|
1679 | @max = leadcoef( MD[ @k, @x ] ); |
---|
1680 | if ( @max != 0 ) |
---|
1681 | { |
---|
1682 | @j = @k; |
---|
1683 | @min = @max; |
---|
1684 | // @k ++; |
---|
1685 | break; // this pure for |
---|
1686 | // continue; |
---|
1687 | // for find min |
---|
1688 | }; |
---|
1689 | }; // for (@k) // found minimal |
---|
1690 | |
---|
1691 | if ( @j == 0 ) |
---|
1692 | { |
---|
1693 | // das ist gut! // |
---|
1694 | @nums = @nums + list ( list ( @x, @y ) ); |
---|
1695 | @x ++; |
---|
1696 | continue; // while |
---|
1697 | } ; |
---|
1698 | |
---|
1699 | for ( @k = @x ; @k <= @v ; @k ++ ) |
---|
1700 | { |
---|
1701 | @t = MD[ @j, @k ]; |
---|
1702 | MD[ @j, @k ] = MD[ @y, @k ]; |
---|
1703 | MD[ @y, @k ] = @t; |
---|
1704 | }; |
---|
1705 | |
---|
1706 | }; // if diag === zero. |
---|
1707 | @ones[@y] = @x; |
---|
1708 | |
---|
1709 | if ( @min == 0 ) |
---|
1710 | { |
---|
1711 | // write_latex_str ( " ******************* ERROR ******************* " ); |
---|
1712 | quit; |
---|
1713 | }; |
---|
1714 | |
---|
1715 | |
---|
1716 | if ( @min != 1) // let's norm the row... to make the '1' ! |
---|
1717 | { |
---|
1718 | @min = 1 / @min; |
---|
1719 | for ( @k = @x; @k <= @v; @k++ ) |
---|
1720 | { |
---|
1721 | @max = leadcoef( MD[@y, @k] ); |
---|
1722 | |
---|
1723 | if ( @max == 0) |
---|
1724 | { |
---|
1725 | @k ++ ; |
---|
1726 | continue; // for : norming the row... |
---|
1727 | }; |
---|
1728 | |
---|
1729 | MD[@y, @k] = @max * @min; // here must be Field ... |
---|
1730 | }; |
---|
1731 | |
---|
1732 | // @min = 1; |
---|
1733 | |
---|
1734 | }; |
---|
1735 | |
---|
1736 | // here must be @min != 0; |
---|
1737 | for ( @k = 1; @k <= @n; @k++ ) |
---|
1738 | { |
---|
1739 | @max = leadcoef( MD[@k, @x] ); |
---|
1740 | |
---|
1741 | if ( ( @k == @y) || (@max == 0) ) |
---|
1742 | { |
---|
1743 | @k ++ ; |
---|
1744 | continue; |
---|
1745 | }; |
---|
1746 | |
---|
1747 | for ( @j = @x; @j <= @v ; @j ++ ) |
---|
1748 | { |
---|
1749 | MD[ @k, @j ] = MD[ @k, @j ] - @max * MD[ @y, @j ]; |
---|
1750 | }; |
---|
1751 | |
---|
1752 | }; //killing |
---|
1753 | |
---|
1754 | @x ++; |
---|
1755 | @y ++; |
---|
1756 | }; //main while. |
---|
1757 | /******************************************************/ |
---|
1758 | |
---|
1759 | // print("Gaussed Matrix: "); |
---|
1760 | // print( MD ); |
---|
1761 | |
---|
1762 | // computation of kernel's basis |
---|
1763 | |
---|
1764 | if ( @x <= @v ) |
---|
1765 | { |
---|
1766 | for ( @k = @x; @k <= @v ; @k ++ ) |
---|
1767 | { |
---|
1768 | @nums = @nums + list ( list ( @k, @n+1 ) ); |
---|
1769 | } |
---|
1770 | } |
---|
1771 | |
---|
1772 | // print("Nums: " ); |
---|
1773 | // print (@nums); |
---|
1774 | |
---|
1775 | list result = list(); |
---|
1776 | |
---|
1777 | // real calculations of the Base of a Ker as V.S. |
---|
1778 | |
---|
1779 | for ( @k = 1; @k <= size(@nums) ; @k ++ ) |
---|
1780 | { |
---|
1781 | @x = @nums[@k][1]; |
---|
1782 | @j = @nums[@k][2]; |
---|
1783 | |
---|
1784 | @t = @given[@x][1]; |
---|
1785 | |
---|
1786 | for ( @y = 1; @y < @j ; @y ++ ) |
---|
1787 | // for every "@x" column |
---|
1788 | { |
---|
1789 | @max = leadcoef( MD[@y, @x] ); |
---|
1790 | if ( (@max != 0) ) |
---|
1791 | { |
---|
1792 | @a = @ones[@y]; |
---|
1793 | @t = @t - @max * @given[@a][1]; |
---|
1794 | }; |
---|
1795 | }; |
---|
1796 | |
---|
1797 | result[@k] = @t; |
---|
1798 | }; |
---|
1799 | |
---|
1800 | ECall( "FSS", "list" ); |
---|
1801 | return ( result ); |
---|
1802 | }; |
---|
1803 | |
---|
1804 | /******************************************************/ |
---|
1805 | static proc reduce_one_per_row(list @given) |
---|
1806 | { |
---|
1807 | BCall( "reduce_one_per_row" ); |
---|
1808 | |
---|
1809 | int @k; |
---|
1810 | int @l = size (@given); |
---|
1811 | |
---|
1812 | if( @l == 0 ) |
---|
1813 | { |
---|
1814 | return (@given); |
---|
1815 | }; |
---|
1816 | |
---|
1817 | int @char = char(basering); |
---|
1818 | |
---|
1819 | intvec @marks; |
---|
1820 | |
---|
1821 | if ( @char == 0 ) |
---|
1822 | { |
---|
1823 | poly @t = poly(0); |
---|
1824 | }; |
---|
1825 | |
---|
1826 | list @unis = list (); |
---|
1827 | poly p; |
---|
1828 | |
---|
1829 | for ( @k = @l; @k > 0; @k -- ) |
---|
1830 | { |
---|
1831 | p = uni_poly( @given[@k][2] ); |
---|
1832 | @unis[@k] = p; |
---|
1833 | |
---|
1834 | if( p == 0 ) |
---|
1835 | { |
---|
1836 | @marks[@k] = 2; |
---|
1837 | } else |
---|
1838 | { |
---|
1839 | @marks[@k] = 1; |
---|
1840 | if ( @char == 0 ) |
---|
1841 | { |
---|
1842 | @t = @t + p; |
---|
1843 | }; |
---|
1844 | }; |
---|
1845 | }; |
---|
1846 | |
---|
1847 | if ( @char != 0 ) |
---|
1848 | { |
---|
1849 | def save = basering; |
---|
1850 | execute( "ring NewRingWithGoodField = (0), (" + varstr(save) + "), (" + ordstr(save) + "); "); |
---|
1851 | poly @t = 0; |
---|
1852 | |
---|
1853 | if(! defined (@unis) ) |
---|
1854 | { |
---|
1855 | list @unis = imap( save, @unis ); |
---|
1856 | }; |
---|
1857 | |
---|
1858 | for ( @k = @l; @k > 0; @k -- ) |
---|
1859 | { |
---|
1860 | if( @marks[@k] == 1 ) |
---|
1861 | { |
---|
1862 | @t = @t + @unis[@k]; |
---|
1863 | }; |
---|
1864 | }; |
---|
1865 | }; |
---|
1866 | |
---|
1867 | int @loop_size = size(@t); |
---|
1868 | poly @for_delete, @tt; |
---|
1869 | int @ll; |
---|
1870 | |
---|
1871 | while( @loop_size > 0 ) |
---|
1872 | { |
---|
1873 | @for_delete = poly(0); |
---|
1874 | @ll = size(@t); |
---|
1875 | |
---|
1876 | for ( @k = @ll; @k > 0; @k -- ) |
---|
1877 | { |
---|
1878 | if ( leadcoef(@t[@k]) == 1 ) |
---|
1879 | { |
---|
1880 | @for_delete = @for_delete + @t[@k]; |
---|
1881 | }; |
---|
1882 | }; |
---|
1883 | |
---|
1884 | @loop_size = size( @for_delete ); |
---|
1885 | |
---|
1886 | if ( @loop_size>0 ) |
---|
1887 | { |
---|
1888 | for( @k = @l ; @k > 0 ; @k -- ) |
---|
1889 | { |
---|
1890 | if ( @marks[@k] == 1) |
---|
1891 | { |
---|
1892 | @tt = @unis[@k]; |
---|
1893 | |
---|
1894 | if( size( @for_delete + @tt ) != ( size( @for_delete ) + size( @tt ) ) ) |
---|
1895 | { |
---|
1896 | @t = @t - @tt; |
---|
1897 | @marks[@k] = 0; |
---|
1898 | }; |
---|
1899 | }; |
---|
1900 | }; |
---|
1901 | }; |
---|
1902 | }; |
---|
1903 | |
---|
1904 | if ( @char != 0 ) |
---|
1905 | { |
---|
1906 | setring(save); |
---|
1907 | kill(NewRingWithGoodField); |
---|
1908 | }; |
---|
1909 | |
---|
1910 | list @reduced = list(); |
---|
1911 | |
---|
1912 | for ( @k = @l ; @k>0 ; @k --) |
---|
1913 | { |
---|
1914 | if (@marks[@k]==2) |
---|
1915 | { |
---|
1916 | @reduced = list ( @given[@k] ) + @reduced ; |
---|
1917 | } else |
---|
1918 | { |
---|
1919 | if (@marks[@k]==1) |
---|
1920 | { |
---|
1921 | @reduced = @reduced + list ( @given[@k] ); |
---|
1922 | }; |
---|
1923 | }; |
---|
1924 | }; |
---|
1925 | |
---|
1926 | ECall( "reduce_one_per_row", "structured list" ); |
---|
1927 | return (@reduced); |
---|
1928 | }; |
---|
1929 | |
---|
1930 | |
---|
1931 | |
---|
1932 | /******************************************************/ |
---|
1933 | static proc calc_base (list @AD_GIVEN) |
---|
1934 | " |
---|
1935 | sort out given 'pbw' into groups due to common monoms in images = ([2])s |
---|
1936 | " |
---|
1937 | { |
---|
1938 | BCall( "calc_base" ); |
---|
1939 | |
---|
1940 | list @AD_CALC = list(); |
---|
1941 | int @ll, @k, @j, @loop_size; |
---|
1942 | |
---|
1943 | poly @t = poly(0); |
---|
1944 | poly @tt, @for_delete; |
---|
1945 | |
---|
1946 | int @t_size; |
---|
1947 | |
---|
1948 | list @GR, @GR_TEMP; |
---|
1949 | |
---|
1950 | int @number = size(@AD_GIVEN); |
---|
1951 | |
---|
1952 | while ( @number > 0 ) |
---|
1953 | { |
---|
1954 | @tt = @AD_GIVEN[@number][2]; |
---|
1955 | if ( size (@tt) == 0) |
---|
1956 | { |
---|
1957 | @AD_CALC = @AD_CALC + list ( @AD_GIVEN[@number][1] ); |
---|
1958 | } else |
---|
1959 | { |
---|
1960 | @t = uni_poly( @tt ); |
---|
1961 | @t_size = size(@t); |
---|
1962 | |
---|
1963 | @GR_TEMP = list (); |
---|
1964 | @GR_TEMP[1] = @t; |
---|
1965 | @GR_TEMP[2] = list ( @AD_GIVEN[@number] ); |
---|
1966 | |
---|
1967 | @loop_size = size(@GR); |
---|
1968 | if ( @loop_size == 0 ) |
---|
1969 | { |
---|
1970 | @GR = list(@GR_TEMP); |
---|
1971 | } else |
---|
1972 | { |
---|
1973 | for ( @k = @loop_size; @k > 0 ; @k -- ) |
---|
1974 | { |
---|
1975 | @tt = @GR[@k][1]; |
---|
1976 | if ( size( @t + @tt ) != ( @t_size + size(@tt) ) ) |
---|
1977 | // whether @tt and @i intersencts? ( will not work in char == 2 !!!) |
---|
1978 | { |
---|
1979 | |
---|
1980 | if ( char(basering) == 0 ) |
---|
1981 | { |
---|
1982 | @GR_TEMP[1] = @GR_TEMP[1] + @tt; |
---|
1983 | } else |
---|
1984 | { |
---|
1985 | @GR_TEMP[1] = uni_poly( @GR_TEMP[1] + @tt ); |
---|
1986 | }; |
---|
1987 | |
---|
1988 | @GR_TEMP[2] = @GR_TEMP[2] + @GR[@k][2]; |
---|
1989 | @GR = delete ( @GR, @k ); |
---|
1990 | }; |
---|
1991 | }; |
---|
1992 | @GR = @GR + list(@GR_TEMP); |
---|
1993 | }; |
---|
1994 | }; |
---|
1995 | @number --; |
---|
1996 | }; // main while |
---|
1997 | |
---|
1998 | list @res; |
---|
1999 | |
---|
2000 | for ( @k = size(@GR); @k > 0 ; @k -- ) |
---|
2001 | { |
---|
2002 | if ( size (@GR[@k][2]) > 1 ) // ! zeroes in AD_CALC so here must be non zero |
---|
2003 | { |
---|
2004 | @res = FSS ( @GR[@k][2], uni_poly(@GR[@k][1])); |
---|
2005 | |
---|
2006 | if ( size (@res) > 0 ) |
---|
2007 | { |
---|
2008 | @AD_CALC = @AD_CALC + @res; |
---|
2009 | }; |
---|
2010 | }; |
---|
2011 | }; |
---|
2012 | |
---|
2013 | ECall( "calc_base" ); |
---|
2014 | return (@AD_CALC); |
---|
2015 | |
---|
2016 | }; |
---|
2017 | |
---|
2018 | |
---|
2019 | /******************************************************/ |
---|
2020 | static proc ApplyAd( list l, list # ) |
---|
2021 | " |
---|
2022 | Apply Ad_{#} to every element of a list of polynomils |
---|
2023 | " |
---|
2024 | { |
---|
2025 | BCall( "ApplyAd", l, # ); |
---|
2026 | |
---|
2027 | if( (size(l) == 0) or (size(#) == 0) ) |
---|
2028 | { |
---|
2029 | return (l); |
---|
2030 | }; |
---|
2031 | |
---|
2032 | poly f, t; |
---|
2033 | |
---|
2034 | if ( typeof(#[1]) == "int") |
---|
2035 | { |
---|
2036 | int next = #[1]; |
---|
2037 | f = my_var (next); |
---|
2038 | } else |
---|
2039 | { |
---|
2040 | if ( typeof(#[1]) == "poly") |
---|
2041 | { |
---|
2042 | f = #[1]; |
---|
2043 | } else |
---|
2044 | { |
---|
2045 | print("Error: cannot differentiate with '" + string(#)+"'"); |
---|
2046 | return (); |
---|
2047 | }; |
---|
2048 | }; |
---|
2049 | |
---|
2050 | for ( int k = size(l); k > 0; k --) |
---|
2051 | { |
---|
2052 | t = l[k]; |
---|
2053 | l[k] = list( t, f * t - t * f ); |
---|
2054 | }; |
---|
2055 | |
---|
2056 | ECall("ApplyAd", l ); |
---|
2057 | return (l); |
---|
2058 | }; |
---|
2059 | |
---|
2060 | |
---|
2061 | /******************************************************/ |
---|
2062 | static proc calc_k_base (list l, list #) |
---|
2063 | " |
---|
2064 | calculation of Kernel of an Ad operator as a K-Vector Space |
---|
2065 | " |
---|
2066 | { |
---|
2067 | BCall( "calc_k_base", "list, " + string(#) ); |
---|
2068 | |
---|
2069 | list t = reduce_one_per_row( l, 0); // optimization (a1) |
---|
2070 | return( QNF_list ( ApplyAd( calc_base(t), # ), 2) ); // calculation of groupps (a2) + gauss. |
---|
2071 | |
---|
2072 | }; |
---|
2073 | |
---|
2074 | LIB "poly.lib"; // for content in proc makeIdeal |
---|
2075 | |
---|
2076 | /******************************************************/ |
---|
2077 | static proc makeIdeal ( list l ) |
---|
2078 | " |
---|
2079 | return: ideal: where the generators are polynomials from list, without 1 and zeroes |
---|
2080 | " |
---|
2081 | { |
---|
2082 | ideal I; poly p; |
---|
2083 | |
---|
2084 | for( int i = 1; i <= size(l); i++ ) |
---|
2085 | { |
---|
2086 | if ( typeof( l[i] ) == "list" ) |
---|
2087 | { |
---|
2088 | p = l[i][1]; // just take the 1st polynom... |
---|
2089 | } else |
---|
2090 | { |
---|
2091 | p = l[i]; |
---|
2092 | }; |
---|
2093 | |
---|
2094 | p = cleardenom( p* (1/content(p)) ); |
---|
2095 | if ( (p != 1) and (p != 0) ) |
---|
2096 | { |
---|
2097 | I = I, p; |
---|
2098 | }; |
---|
2099 | }; |
---|
2100 | |
---|
2101 | I = simplify ( I, 2 ); // no zeroes |
---|
2102 | |
---|
2103 | return(I); |
---|
2104 | }; |
---|
2105 | |
---|
2106 | /******************************************************/ |
---|
2107 | static proc inCenter_poly( poly p ) |
---|
2108 | " |
---|
2109 | if p in center => return 1 |
---|
2110 | otherwise return 0 |
---|
2111 | " |
---|
2112 | { |
---|
2113 | poly t; |
---|
2114 | for (int k = nvars(basering); k>0; k-- ) |
---|
2115 | { |
---|
2116 | t = var(k); |
---|
2117 | if ( QNF_poly(t * p - p * t) != 0 ) |
---|
2118 | { |
---|
2119 | if( toprint() ) |
---|
2120 | { |
---|
2121 | "POLY: ", string (p), " is NOT in center"; |
---|
2122 | }; |
---|
2123 | return (0); |
---|
2124 | }; |
---|
2125 | }; |
---|
2126 | |
---|
2127 | if( toprint() ) |
---|
2128 | { |
---|
2129 | "POLY: ", string (p), " is in center"; |
---|
2130 | }; |
---|
2131 | return (1); |
---|
2132 | }; |
---|
2133 | |
---|
2134 | /******************************************************/ |
---|
2135 | static proc inCenter_list(def l) |
---|
2136 | { |
---|
2137 | for ( int @i = 1; @i <= size(l); @i++ ) |
---|
2138 | { |
---|
2139 | if ( typeof(l[@i])=="poly" ) |
---|
2140 | { |
---|
2141 | if (! inCenter_poly(l[@i]) ) |
---|
2142 | { |
---|
2143 | return(0); |
---|
2144 | }; |
---|
2145 | |
---|
2146 | } else |
---|
2147 | { |
---|
2148 | if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") ) |
---|
2149 | { |
---|
2150 | if (! inCenter_list(l[@i]) ) |
---|
2151 | { |
---|
2152 | return(0); |
---|
2153 | }; |
---|
2154 | }; |
---|
2155 | }; |
---|
2156 | }; |
---|
2157 | return(1); |
---|
2158 | }; |
---|
2159 | |
---|
2160 | /******************************************************/ |
---|
2161 | static proc inCentralizer_poly( poly p, poly f ) |
---|
2162 | " |
---|
2163 | if p in centralizer(f) => return 1 |
---|
2164 | otherwise return 0 |
---|
2165 | " |
---|
2166 | { |
---|
2167 | if ( QNF_poly(f * p - p * f) != 0 ) |
---|
2168 | { |
---|
2169 | if( toprint() ) |
---|
2170 | { |
---|
2171 | "POLY: ", string (p), " is NOT in centralizer(f)"; |
---|
2172 | }; |
---|
2173 | return (0); |
---|
2174 | }; |
---|
2175 | |
---|
2176 | if( toprint() ) |
---|
2177 | { |
---|
2178 | "POLY: ", string (p), " is in centralizer(f)"; |
---|
2179 | }; |
---|
2180 | return (1); |
---|
2181 | }; |
---|
2182 | |
---|
2183 | /******************************************************/ |
---|
2184 | static proc inCentralizer_list( def l, poly f ) |
---|
2185 | { |
---|
2186 | for ( int @i = 1; @i <= size(l); @i++ ) |
---|
2187 | { |
---|
2188 | if ( typeof(l[@i])=="poly" ) |
---|
2189 | { |
---|
2190 | if (! inCentralizer_poly(l[@i], f) ) |
---|
2191 | { |
---|
2192 | return(0); |
---|
2193 | }; |
---|
2194 | |
---|
2195 | } else |
---|
2196 | { |
---|
2197 | if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") ) |
---|
2198 | { |
---|
2199 | if (! inCentralizer_list(l[@i], f) ) |
---|
2200 | { |
---|
2201 | return(0); |
---|
2202 | }; |
---|
2203 | }; |
---|
2204 | }; |
---|
2205 | }; |
---|
2206 | return(1); |
---|
2207 | }; |
---|
2208 | |
---|
2209 | |
---|
2210 | /******************************************************/ |
---|
2211 | /******************************************************/ |
---|
2212 | // main algorithms |
---|
2213 | |
---|
2214 | |
---|
2215 | /******************************************************/ |
---|
2216 | /******************************************************/ |
---|
2217 | // center's computation |
---|
2218 | |
---|
2219 | /******************************************************/ |
---|
2220 | /* static */ proc center_min_iterative( int MaxDeg, list # ) |
---|
2221 | " |
---|
2222 | computes the 'minimal' set of central elements (of deg <= MaxDeg) in iterative way |
---|
2223 | Note: based on calc_k_base, zAddBad, zRefine, zCancel, PBW_* |
---|
2224 | " |
---|
2225 | { |
---|
2226 | BCall("center_min_iterative", MaxDeg, #); |
---|
2227 | |
---|
2228 | int n = myInt(#); |
---|
2229 | int m = ( MaxDeg < 0 ); // no bound on Degree |
---|
2230 | |
---|
2231 | int MinDeg = 6; // starting guess for MaxDeg |
---|
2232 | int Delta = 4; // increment of MaxDeg |
---|
2233 | |
---|
2234 | if( m ) |
---|
2235 | { |
---|
2236 | // minimal guess |
---|
2237 | MaxDeg = MinDeg; |
---|
2238 | }; |
---|
2239 | |
---|
2240 | list @q; int @i; int N = nvars(basering); |
---|
2241 | |
---|
2242 | my_var_init(); // setup for my_var(n) |
---|
2243 | QNF_init(); |
---|
2244 | |
---|
2245 | list PBW = list(); |
---|
2246 | list PBW_PLAIN = list(); |
---|
2247 | |
---|
2248 | PBW[ index(0) ] = PBW_init(); |
---|
2249 | |
---|
2250 | list SYS = PBW_sys( my_vars(), my_var(1) ); |
---|
2251 | |
---|
2252 | |
---|
2253 | list @z = list (); // center list |
---|
2254 | list @l = init_bads( MaxDeg ); // verbotten leadexps... |
---|
2255 | |
---|
2256 | @q = PBW[ index(0) ]; |
---|
2257 | |
---|
2258 | int k = 1; |
---|
2259 | while( k <= ( MaxDeg+1 ) ) |
---|
2260 | { |
---|
2261 | for ( @i = 2; @i <= N; @i ++ ) |
---|
2262 | { |
---|
2263 | @q = calc_k_base (@q, my_var(@i)); |
---|
2264 | }; |
---|
2265 | |
---|
2266 | @q = zRefine (calc_k_base(@q)); // new center! |
---|
2267 | |
---|
2268 | |
---|
2269 | if ( size(@q) > 0 ) |
---|
2270 | { |
---|
2271 | @z = @z + @q; // computed central elements |
---|
2272 | |
---|
2273 | if( (n > 0) and (size(@z) > n) ) |
---|
2274 | { |
---|
2275 | break; // found all central elements |
---|
2276 | }; |
---|
2277 | }; |
---|
2278 | |
---|
2279 | if( k == ( MaxDeg+1 ) ) |
---|
2280 | { |
---|
2281 | if( (n == 0) and ( !m ) ) |
---|
2282 | { |
---|
2283 | break; // that's all |
---|
2284 | }; |
---|
2285 | |
---|
2286 | MaxDeg = MaxDeg + Delta; |
---|
2287 | |
---|
2288 | // renew bad list |
---|
2289 | @l = init_bads( MaxDeg ); |
---|
2290 | @l = zAddBad( @l, @z, MaxDeg, 0 ); |
---|
2291 | |
---|
2292 | } else |
---|
2293 | { |
---|
2294 | |
---|
2295 | if ( size(@q) > 0 ) |
---|
2296 | { |
---|
2297 | @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' ! |
---|
2298 | }; |
---|
2299 | |
---|
2300 | }; |
---|
2301 | |
---|
2302 | PBW[index(k)] = PBW_next( PBW, k, SYS ); |
---|
2303 | |
---|
2304 | PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l ); |
---|
2305 | |
---|
2306 | @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d] |
---|
2307 | |
---|
2308 | k ++; |
---|
2309 | }; |
---|
2310 | |
---|
2311 | QNF_done(); |
---|
2312 | my_var_done(); |
---|
2313 | |
---|
2314 | return (makeIdeal(@z)); |
---|
2315 | }; |
---|
2316 | |
---|
2317 | /******************************************************/ |
---|
2318 | static proc center_vectorspace( int MaxDeg ) |
---|
2319 | " |
---|
2320 | pure calculation of center as a finitely dimensional Vector Space (deg <= MaxDeg ) |
---|
2321 | " |
---|
2322 | { |
---|
2323 | my_var_init(); |
---|
2324 | |
---|
2325 | int N = nvars( basering ); |
---|
2326 | list P = PBW_base( MaxDeg, my_var(1) ); |
---|
2327 | |
---|
2328 | for( int v = 2; v <= N; v++ ) |
---|
2329 | { |
---|
2330 | P = calc_k_base( P, my_var(v) ); |
---|
2331 | }; |
---|
2332 | |
---|
2333 | my_var_done(); |
---|
2334 | |
---|
2335 | return( calc_k_base ( P ) ); |
---|
2336 | }; |
---|
2337 | |
---|
2338 | /******************************************************/ |
---|
2339 | /* static */ proc center_min_vectorspace( int MaxDeg ) |
---|
2340 | " |
---|
2341 | computes the 'minimal' set of central elements (of deg <= MaxDeg) |
---|
2342 | by reducing the set of it's generators as vector space |
---|
2343 | |
---|
2344 | Note: based on center_vectorspace. |
---|
2345 | Note: reduction by zReduce. |
---|
2346 | " |
---|
2347 | { |
---|
2348 | |
---|
2349 | QNF_init(); |
---|
2350 | ideal res = makeIdeal( zReduce( center_vectorspace( MaxDeg))); |
---|
2351 | QNF_done(); |
---|
2352 | |
---|
2353 | return( res ); |
---|
2354 | }; |
---|
2355 | |
---|
2356 | |
---|
2357 | /******************************************************/ |
---|
2358 | /******************************************************/ |
---|
2359 | // centralizer's computations |
---|
2360 | |
---|
2361 | /******************************************************/ |
---|
2362 | static proc centralizer_vectorspace( poly p, int MaxDeg ) |
---|
2363 | { |
---|
2364 | return( calc_k_base( PBW_base( MaxDeg, p ))); |
---|
2365 | }; |
---|
2366 | |
---|
2367 | /******************************************************/ |
---|
2368 | /* static */ proc centralizer_min_vectorspace( poly p, int MaxDeg ) |
---|
2369 | { |
---|
2370 | |
---|
2371 | QNF_init(); |
---|
2372 | ideal res = makeIdeal( zReduce( centralizer_vectorspace( p, MaxDeg))); |
---|
2373 | QNF_done(); |
---|
2374 | |
---|
2375 | return( res ); |
---|
2376 | }; |
---|
2377 | |
---|
2378 | /******************************************************/ |
---|
2379 | /* static */ proc centralizer_min_iterative( poly p, int MaxDeg, list # ) |
---|
2380 | " |
---|
2381 | computes the 'minimal' set of elements (of deg <= MaxDeg) generating centralizer of p in iterative way |
---|
2382 | Note: based on calc_k_base |
---|
2383 | |
---|
2384 | !!! NEED DEBUG !!! |
---|
2385 | Note: no proof that it is really centralizer and moreover 'minimal' centralizer |
---|
2386 | " |
---|
2387 | { |
---|
2388 | QNF_init(); |
---|
2389 | |
---|
2390 | int n = myInt(#); |
---|
2391 | int m = (MaxDeg < 0); |
---|
2392 | |
---|
2393 | int MinDeg = 6; // starting guess for MaxDeg |
---|
2394 | int Delta = 4; // increment of MaxDeg |
---|
2395 | |
---|
2396 | if( m ) |
---|
2397 | { |
---|
2398 | // minimal guess |
---|
2399 | MaxDeg = MinDeg; |
---|
2400 | }; |
---|
2401 | |
---|
2402 | list @q; |
---|
2403 | |
---|
2404 | list PBW = list(); |
---|
2405 | list PBW_PLAIN = list(); |
---|
2406 | |
---|
2407 | PBW[ index(0) ] = PBW_init(); |
---|
2408 | |
---|
2409 | list SYS = PBW_sys( my_vars(), p ); |
---|
2410 | |
---|
2411 | list @z = list (); // result list |
---|
2412 | list @l = init_bads( MaxDeg ); // verbotten loeadexps... |
---|
2413 | |
---|
2414 | @q = PBW[ index(0) ]; |
---|
2415 | |
---|
2416 | for (int k = 1; k <= ( MaxDeg+1 ); k ++ ) |
---|
2417 | { |
---|
2418 | @q = zRefine( calc_k_base(@q), 1 ); |
---|
2419 | |
---|
2420 | if ( size(@q) > 0 ) |
---|
2421 | { |
---|
2422 | @z = @z + @q; // computed desired elements |
---|
2423 | |
---|
2424 | if( (n > 0) and (size(@z) > n) ) |
---|
2425 | { |
---|
2426 | break; // found needed elements |
---|
2427 | }; |
---|
2428 | }; |
---|
2429 | |
---|
2430 | if( k == ( MaxDeg+1 ) ) |
---|
2431 | { |
---|
2432 | if( (n == 0) or ( !m ) ) |
---|
2433 | { |
---|
2434 | break; // that's all |
---|
2435 | }; |
---|
2436 | |
---|
2437 | MaxDeg = MaxDeg + Delta; |
---|
2438 | |
---|
2439 | // renew bad list |
---|
2440 | @l = init_bads( MaxDeg ); |
---|
2441 | @l = zAddBad( @l, @z, MaxDeg, 0 ); |
---|
2442 | |
---|
2443 | } else |
---|
2444 | { |
---|
2445 | |
---|
2446 | if ( size(@q) > 0 ) |
---|
2447 | { |
---|
2448 | @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' ! |
---|
2449 | }; |
---|
2450 | |
---|
2451 | }; |
---|
2452 | |
---|
2453 | PBW[index(k)] = PBW_next( PBW, k, SYS ); |
---|
2454 | PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l ); |
---|
2455 | @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d] |
---|
2456 | |
---|
2457 | }; |
---|
2458 | |
---|
2459 | QNF_done(); |
---|
2460 | return (makeIdeal(@z)); |
---|
2461 | }; |
---|
2462 | |
---|
2463 | |
---|
2464 | |
---|
2465 | /******************************************************/ |
---|
2466 | /******************************************************/ |
---|
2467 | // main aliases |
---|
2468 | |
---|
2469 | /******************************************************/ |
---|
2470 | proc inCenter( def a ) |
---|
2471 | "USAGE: inCenter(a); a poly/list/ideal |
---|
2472 | RETURN: integer (1 if a in the center, 0 otherwise) |
---|
2473 | EXAMPLE: example inCenter; shows examples" |
---|
2474 | { |
---|
2475 | QNF_init(); |
---|
2476 | def res; |
---|
2477 | |
---|
2478 | if ( typeof(a) == "poly" ) |
---|
2479 | { |
---|
2480 | res = inCenter_poly(a); |
---|
2481 | } else |
---|
2482 | { |
---|
2483 | if ( (typeof(a)=="list") or (typeof(a)=="ideal") ) |
---|
2484 | { |
---|
2485 | res = inCenter_list(a); |
---|
2486 | } else |
---|
2487 | { |
---|
2488 | res = a; |
---|
2489 | }; |
---|
2490 | }; |
---|
2491 | |
---|
2492 | QNF_done(); |
---|
2493 | return (res); |
---|
2494 | } |
---|
2495 | example |
---|
2496 | { |
---|
2497 | "EXAMPLE:";echo=2; |
---|
2498 | ring r=0,(x,y,z),dp; |
---|
2499 | matrix D[3][3]=0; |
---|
2500 | D[1,2]=-z; |
---|
2501 | D[1,3]=2*x; |
---|
2502 | D[2,3]=-2*y; |
---|
2503 | ncalgebra(1,D); // this is U(sl_2) |
---|
2504 | poly p=4*x*y+z^2-2*z; |
---|
2505 | inCenter(p); |
---|
2506 | poly f=4*x*y; |
---|
2507 | inCenter(f); |
---|
2508 | list l= list( 1, p, p^2, p^3); |
---|
2509 | inCenter(l); |
---|
2510 | ideal I= p, f; |
---|
2511 | inCenter(I); |
---|
2512 | }; |
---|
2513 | |
---|
2514 | |
---|
2515 | /******************************************************************************/ |
---|
2516 | proc inCentralizer( def a, poly f ) |
---|
2517 | "USAGE: inCentralizer(a, f); a poly/list/ideal, f poly |
---|
2518 | RETURN: integer (1 if a in the centralizer(f), 0 otherwise) |
---|
2519 | EXAMPLE: example inCentralizer; shows examples" |
---|
2520 | { |
---|
2521 | QNF_init(); |
---|
2522 | def res; |
---|
2523 | |
---|
2524 | if ( typeof(a) == "poly" ) |
---|
2525 | { |
---|
2526 | res = inCentralizer_poly(a, f); |
---|
2527 | } else |
---|
2528 | { |
---|
2529 | if ( (typeof(a)=="list") or (typeof(a)=="ideal") ) |
---|
2530 | { |
---|
2531 | res = inCentralizer_list(a, f); |
---|
2532 | } else |
---|
2533 | { |
---|
2534 | res = a; |
---|
2535 | }; |
---|
2536 | }; |
---|
2537 | |
---|
2538 | QNF_done(); |
---|
2539 | return (res); |
---|
2540 | } |
---|
2541 | example |
---|
2542 | { |
---|
2543 | "EXAMPLE:";echo=2; |
---|
2544 | ring r=0,(x,y,z),dp; |
---|
2545 | matrix D[3][3]=0; |
---|
2546 | D[1,2]=-z; |
---|
2547 | ncalgebra(1,D); // the Heisenberg algebra |
---|
2548 | poly f = x^2; |
---|
2549 | poly a = z; // lies in center |
---|
2550 | poly b = y^2; |
---|
2551 | inCentralizer(a, f); |
---|
2552 | inCentralizer(b, f); |
---|
2553 | list l = list(1, a); |
---|
2554 | inCentralizer(l, f); |
---|
2555 | ideal I = a, b; |
---|
2556 | inCentralizer(I, f); |
---|
2557 | }; |
---|
2558 | |
---|
2559 | /******************************************************/ |
---|
2560 | proc center(int MaxDeg, list # ) |
---|
2561 | "USAGE: center(MaxDeg[, N]); int MaxDeg, int N |
---|
2562 | RETURN: ideal generated by found elements |
---|
2563 | NOTE: computes the 'minimal' set of central elements. |
---|
2564 | Since in general algorithms knows nothing about the number and degrees of |
---|
2565 | desired polynomials one have to specify any kind of termination condition: |
---|
2566 | 1. MaxDeg - maximal degree of desired elements or/and |
---|
2567 | 2. n - the minimal number of desired elements to find. |
---|
2568 | SEE ALSO: centralizer; inCenter |
---|
2569 | EXAMPLE: example center; shows an example |
---|
2570 | " |
---|
2571 | { |
---|
2572 | if ( myInt(#) > 0 ) // given a number of central elements to compute |
---|
2573 | { |
---|
2574 | return ( center_min_iterative( MaxDeg, # ) ); |
---|
2575 | }; |
---|
2576 | |
---|
2577 | if( MaxDeg >= 0 ) // standard computation - the fastest one (often) |
---|
2578 | { |
---|
2579 | // return ( center_min_iterative( MaxDeg, # ) ); |
---|
2580 | return ( center_min_vectorspace ( MaxDeg ) ); |
---|
2581 | }; |
---|
2582 | |
---|
2583 | print( "Error: wrong arguments." ); |
---|
2584 | return(); |
---|
2585 | |
---|
2586 | } |
---|
2587 | example |
---|
2588 | { "EXAMPLE:"; echo = 2; |
---|
2589 | ring a=0,(x,y,z),dp; |
---|
2590 | matrix D[3][3]=0; |
---|
2591 | D[1,2]=-z; |
---|
2592 | D[1,3]=2*x; |
---|
2593 | D[2,3]=-2*y; |
---|
2594 | ncalgebra(1,D); // this is U(sl_2) |
---|
2595 | ideal Z = center(2); // find all central elements of degree <= 2 |
---|
2596 | Z; |
---|
2597 | inCenter(Z); |
---|
2598 | ideal ZZ = center(-1, 1 ); // find the first non trivial central element |
---|
2599 | ZZ; ""; |
---|
2600 | inCenter(ZZ); |
---|
2601 | }; |
---|
2602 | |
---|
2603 | |
---|
2604 | /******************************************************/ |
---|
2605 | proc centralizer( poly p, int MaxDeg, list # ) |
---|
2606 | "USAGE: centralizer(F, MaxDeg[, N]); poly F, int MaxDeg, int N |
---|
2607 | RETURN: ideal generated by found elements |
---|
2608 | NOTE: computes the 'minimal' set of elements of centralizer(F). |
---|
2609 | Since in general algorithms knows nothing about the number and degrees of |
---|
2610 | desired polynomials one have to specify any kind of termination condition: |
---|
2611 | 1. MaxDeg - maximal degree of desired elements or/and |
---|
2612 | 2. n - the minimal number of desired elements to find. |
---|
2613 | SEE ALSO: center; inCentralizer |
---|
2614 | EXAMPLE: example centralizer; shows an example" |
---|
2615 | { |
---|
2616 | |
---|
2617 | if( myInt(#) > 0 ) |
---|
2618 | { |
---|
2619 | return( centralizer_min_iterative( p, MaxDeg, # ) ); |
---|
2620 | }; |
---|
2621 | |
---|
2622 | if( MaxDeg >= 0 ) |
---|
2623 | { |
---|
2624 | // return( centralizer_min_iterative( p, MaxDeg, # ) ); |
---|
2625 | return( centralizer_min_vectorspace( p, MaxDeg ) ); |
---|
2626 | }; |
---|
2627 | |
---|
2628 | print( "Error: wrong arguments." ); |
---|
2629 | return(); |
---|
2630 | |
---|
2631 | } |
---|
2632 | example |
---|
2633 | { "EXAMPLE:"; echo = 2; |
---|
2634 | ring a=0,(x,y,z),dp; |
---|
2635 | matrix D[3][3]=0; |
---|
2636 | D[1,2]=-z; |
---|
2637 | D[1,3]=2*x; |
---|
2638 | D[2,3]=-2*y; |
---|
2639 | ncalgebra(1,D); // this is U(sl_2) |
---|
2640 | poly f = 4*x*y+z^2-2*z; // central polynomial |
---|
2641 | f; |
---|
2642 | ideal c = centralizer(f, 2); // find all elements of degree <= 2 which lies in centralizer of f |
---|
2643 | c; |
---|
2644 | inCentralizer(c, f); |
---|
2645 | ideal cc = centralizer(f, -1, 2 ); // find at least first two non trivial elements of the centralizer of f |
---|
2646 | cc; |
---|
2647 | inCentralizer(cc, f); |
---|
2648 | poly g = z^2-2*z; // any polynomial |
---|
2649 | g; ""; |
---|
2650 | c = centralizer(g, 2); // find all elements of degree <= 2 which lies in centralizer of g |
---|
2651 | c; ""; |
---|
2652 | inCentralizer(c, g); |
---|
2653 | cc = centralizer(g, -1, 2 ); // find at least first two non trivial elements of the centralizer of g |
---|
2654 | cc; ""; |
---|
2655 | inCentralizer(cc, g); |
---|
2656 | }; |
---|