1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /***************************************************************
5 * File: sca.cc
6 * Purpose: supercommutative kernel procedures
7 * Author: motsak (Oleksandr Motsak)
8 * Created: 2006/12/18
9 *******************************************************************/
10
11 // set it here if needed.
12 #define OUTPUT 0
13 #define MYTEST 0
14
15 #if MYTEST
16 #define OM_CHECK 4
17 #define OM_TRACK 5
18 #endif
19
20 // #define PDEBUG 2
21
22
23
24 #include "misc/auxiliary.h"
25
26 #ifdef HAVE_PLURAL
27
28 // for
29 #define PLURAL_INTERNAL_DECLARATIONS
30
31 #include "polys/nc/sca.h"
32 #include "polys/nc/nc.h"
33 #include "polys/nc/gb_hack.h"
34
35 #include "coeffs/numbers.h"
36
37 #include "misc/options.h"
38
39 #include "polys/monomials/p_polys.h"
40
41 #include "polys/simpleideals.h"
42 #include "misc/intvec.h"
43
44 #include "polys/monomials/ring.h"
45 #include "polys/kbuckets.h"
46 #include "polys/sbuckets.h"
47
48 #include "polys/prCopy.h"
49
50 #include "polys/operations/p_Mult_q.h"
51 #include "polys/templates/p_MemAdd.h"
52
53 #include "polys/weight.h"
54
55 // poly functions defined in p_Procs :
56
57 // return pPoly * pMonom; preserve pPoly and pMonom.
58 poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing, poly &);
59
60 // return pMonom * pPoly; preserve pPoly and pMonom.
61 static poly sca_pp_mm_Mult(const poly pPoly, const poly pMonom, const ring rRing);
62
63 // return pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
64 poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing);
65
66 // return pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
67 static poly sca_p_mm_Mult(poly pPoly, const poly pMonom, const ring rRing);
68
69
70 // compute the spolynomial of p1 and p2
71 poly sca_SPoly(const poly p1, const poly p2, const ring r);
72 poly sca_ReduceSpoly(const poly p1, poly p2, const ring r);
73
74
75 ////////////////////////////////////////////////////////////////////////////////////////////////////
76 // Super Commutative Algebra extension by Oleksandr
77 ////////////////////////////////////////////////////////////////////////////////////////////////////
78
79 // returns the sign of: lm(pMonomM) * lm(pMonomMM),
80 // preserves input, may return +/-1, 0
sca_Sign_mm_Mult_mm(const poly pMonomM,const poly pMonomMM,const ring rRing)81 static inline int sca_Sign_mm_Mult_mm( const poly pMonomM, const poly pMonomMM, const ring rRing )
82 {
83 #ifdef PDEBUG
84 p_Test(pMonomM, rRing);
85 p_Test(pMonomMM, rRing);
86 #endif
87
88 const short iFirstAltVar = scaFirstAltVar(rRing);
89 const short iLastAltVar = scaLastAltVar(rRing);
90
91 REGISTER unsigned int tpower = 0;
92 REGISTER unsigned int cpower = 0;
93
94 for( REGISTER short j = iLastAltVar; j >= iFirstAltVar; j-- )
95 {
96 const unsigned int iExpM = p_GetExp(pMonomM, j, rRing);
97 const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
98
99 #ifdef PDEBUG
100 assume( iExpM <= 1);
101 assume( iExpMM <= 1);
102 #endif
103
104 if( iExpMM != 0 ) // TODO: think about eliminating there if-s...
105 {
106 if( iExpM != 0 )
107 {
108 return 0; // lm(pMonomM) * lm(pMonomMM) == 0
109 }
110 tpower ^= cpower; // compute degree of (-1).
111 }
112 cpower ^= iExpM;
113 }
114
115 #ifdef PDEBUG
116 assume(tpower <= 1);
117 #endif
118
119 // 1 => -1 // degree is odd => negate coeff.
120 // 0 => 1
121
122 return(1 - (tpower << 1) );
123 }
124
125
126
127
128 // returns and changes pMonomM: lt(pMonomM) = lt(pMonomM) * lt(pMonomMM),
129 // preserves pMonomMM. may return NULL!
130 // if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomM) * lt(pMonomMM)
131 // if result == NULL => pMonomM MUST BE deleted manually!
sca_m_Mult_mm(poly pMonomM,const poly pMonomMM,const ring rRing)132 static inline poly sca_m_Mult_mm( poly pMonomM, const poly pMonomMM, const ring rRing )
133 {
134 #ifdef PDEBUG
135 p_Test(pMonomM, rRing);
136 p_Test(pMonomMM, rRing);
137 #endif
138
139 const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
140 const unsigned int iLastAltVar = scaLastAltVar(rRing);
141
142 REGISTER unsigned int tpower = 0;
143 REGISTER unsigned int cpower = 0;
144
145 for( REGISTER unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
146 {
147 const unsigned int iExpM = p_GetExp(pMonomM, j, rRing);
148 const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
149
150 #ifdef PDEBUG
151 assume( iExpM <= 1);
152 assume( iExpMM <= 1);
153 #endif
154
155 if( iExpMM != 0 )
156 {
157 if( iExpM != 0 ) // result is zero!
158 {
159 return NULL; // we do nothing with pMonomM in this case!
160 }
161
162 tpower ^= cpower; // compute degree of (-1).
163 }
164
165 cpower ^= iExpM;
166 }
167
168 #ifdef PDEBUG
169 assume(tpower <= 1);
170 #endif
171
172 p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
173
174 number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
175
176 if( (tpower) != 0 ) // degree is odd => negate coeff.
177 nCoeffM = n_InpNeg(nCoeffM, rRing->cf); // negate nCoeff (will destroy the original number)
178
179 const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
180
181 number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing->cf); // new number!
182
183 p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
184
185 #ifdef PDEBUG
186 p_LmTest(pMonomM, rRing);
187 #endif
188
189 return(pMonomM);
190 }
191
192 // returns and changes pMonomM: lt(pMonomM) = lt(pMonomMM) * lt(pMonomM),
193 // preserves pMonomMM. may return NULL!
194 // if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomMM) * lt(pMonomM)
195 // if result == NULL => pMonomM MUST BE deleted manually!
sca_mm_Mult_m(const poly pMonomMM,poly pMonomM,const ring rRing)196 static inline poly sca_mm_Mult_m( const poly pMonomMM, poly pMonomM, const ring rRing )
197 {
198 #ifdef PDEBUG
199 p_Test(pMonomM, rRing);
200 p_Test(pMonomMM, rRing);
201 #endif
202
203 const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
204 const unsigned int iLastAltVar = scaLastAltVar(rRing);
205
206 REGISTER unsigned int tpower = 0;
207 REGISTER unsigned int cpower = 0;
208
209 for( REGISTER unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
210 {
211 const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
212 const unsigned int iExpM = p_GetExp(pMonomM, j, rRing);
213
214 #ifdef PDEBUG
215 assume( iExpM <= 1);
216 assume( iExpMM <= 1);
217 #endif
218
219 if( iExpM != 0 )
220 {
221 if( iExpMM != 0 ) // result is zero!
222 {
223 return NULL; // we do nothing with pMonomM in this case!
224 }
225
226 tpower ^= cpower; // compute degree of (-1).
227 }
228
229 cpower ^= iExpMM;
230 }
231
232 #ifdef PDEBUG
233 assume(tpower <= 1);
234 #endif
235
236 p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
237
238 number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
239
240 if( (tpower) != 0 ) // degree is odd => negate coeff.
241 nCoeffM = n_InpNeg(nCoeffM, rRing->cf); // negate nCoeff (will destroy the original number), creates new number!
242
243 const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
244
245 number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing->cf); // new number!
246
247 p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
248
249 #ifdef PDEBUG
250 p_LmTest(pMonomM, rRing);
251 #endif
252
253 return(pMonomM);
254 }
255
256
257
258 // returns: result = lt(pMonom1) * lt(pMonom2),
259 // preserves pMonom1, pMonom2. may return NULL!
260 // if result != NULL => next(result) = NULL, lt(result) = lt(pMonom1) * lt(pMonom2)
sca_mm_Mult_mm(poly pMonom1,const poly pMonom2,const ring rRing)261 static inline poly sca_mm_Mult_mm( poly pMonom1, const poly pMonom2, const ring rRing )
262 {
263 #ifdef PDEBUG
264 p_Test(pMonom1, rRing);
265 p_Test(pMonom2, rRing);
266 #endif
267
268 const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
269 const unsigned int iLastAltVar = scaLastAltVar(rRing);
270
271 REGISTER unsigned int tpower = 0;
272 REGISTER unsigned int cpower = 0;
273
274 for( REGISTER unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
275 {
276 const unsigned int iExp1 = p_GetExp(pMonom1, j, rRing);
277 const unsigned int iExp2 = p_GetExp(pMonom2, j, rRing);
278
279 #ifdef PDEBUG
280 assume( iExp1 <= 1);
281 assume( iExp2 <= 1);
282 #endif
283
284 if( iExp2 != 0 )
285 {
286 if( iExp1 != 0 ) // result is zero!
287 {
288 return NULL;
289 }
290 tpower ^= cpower; // compute degree of (-1).
291 }
292 cpower ^= iExp1;
293 }
294
295 #ifdef PDEBUG
296 assume(cpower <= 1);
297 #endif
298
299 poly pResult;
300 p_AllocBin(pResult,rRing->PolyBin,rRing);
301
302 p_ExpVectorSum(pResult, pMonom1, pMonom2, rRing); // "exponents" are additive!!!
303
304 pNext(pResult) = NULL;
305
306 const number nCoeff1 = p_GetCoeff(pMonom1, rRing); // no new copy!
307 const number nCoeff2 = p_GetCoeff(pMonom2, rRing); // no new copy!
308
309 number nCoeff = n_Mult(nCoeff1, nCoeff2, rRing->cf); // new number!
310
311 if( (tpower) != 0 ) // degree is odd => negate coeff.
312 nCoeff = n_InpNeg(nCoeff, rRing->cf); // negate nCoeff (will destroy the original number)
313
314 p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
315
316 #ifdef PDEBUG
317 p_LmTest(pResult, rRing);
318 #endif
319
320 return(pResult);
321 }
322
323 // returns: result = x_i * lt(pMonom),
324 // preserves pMonom. may return NULL!
325 // if result != NULL => next(result) = NULL, lt(result) = x_i * lt(pMonom)
sca_xi_Mult_mm(short i,const poly pMonom,const ring rRing)326 static inline poly sca_xi_Mult_mm(short i, const poly pMonom, const ring rRing)
327 {
328 #ifdef PDEBUG
329 p_Test(pMonom, rRing);
330 #endif
331
332 assume( i <= scaLastAltVar(rRing));
333 assume( scaFirstAltVar(rRing) <= i );
334
335 if( p_GetExp(pMonom, i, rRing) != 0 ) // => result is zero!
336 return NULL;
337
338 const short iFirstAltVar = scaFirstAltVar(rRing);
339
340 REGISTER unsigned int cpower = 0;
341
342 for( REGISTER short j = iFirstAltVar; j < i ; j++ )
343 cpower ^= p_GetExp(pMonom, j, rRing);
344
345 #ifdef PDEBUG
346 assume(cpower <= 1);
347 #endif
348
349 poly pResult = p_LmInit(pMonom, rRing);
350
351 p_SetExp(pResult, i, 1, rRing); // pResult*=X_i &&
352 p_Setm(pResult, rRing); // addjust degree after previous step!
353
354 number nCoeff = n_Copy(pGetCoeff(pMonom), rRing->cf); // new number!
355
356 if( cpower != 0 ) // degree is odd => negate coeff.
357 nCoeff = n_InpNeg(nCoeff, rRing->cf); // negate nCoeff (will destroy the original number)
358
359 p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
360
361 #ifdef PDEBUG
362 p_LmTest(pResult, rRing);
363 #endif
364
365 return(pResult);
366 }
367
368 //-----------------------------------------------------------------------------------//
369
370 // return poly = pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
sca_p_Mult_mm(poly pPoly,const poly pMonom,const ring rRing)371 poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing)
372 {
373 assume( rIsSCA(rRing) );
374
375 #ifdef PDEBUG
376 // PrintS("sca_p_Mult_mm\n"); // !
377
378 p_Test(pPoly, rRing);
379 p_Test(pMonom, rRing);
380 #endif
381
382 if( pPoly == NULL )
383 return NULL;
384
385 assume(pMonom !=NULL);
386 //if( pMonom == NULL )
387 //{
388 // // pPoly != NULL =>
389 // p_Delete( &pPoly, rRing );
390 // return NULL;
391 //}
392
393 const int iComponentMonomM = p_GetComp(pMonom, rRing);
394
395 poly p = pPoly; poly* ppPrev = &pPoly;
396
397 loop
398 {
399 #ifdef PDEBUG
400 p_Test(p, rRing);
401 #endif
402 const int iComponent = p_GetComp(p, rRing);
403
404 if( iComponent!=0 )
405 {
406 if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
407 {
408 // REPORT_ERROR
409 Werror("sca_p_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
410 // what should we do further?!?
411
412 p_Delete( &pPoly, rRing); // delete the result AND rest
413 return NULL;
414 }
415 #ifdef PDEBUG
416 if(iComponentMonomM==0 )
417 {
418 dReportError("sca_p_Mult_mm: Multiplication in the left module from the right");
419 }
420 #endif
421 }
422
423 // terms will be in the same order because of quasi-ordering!
424 poly v = sca_m_Mult_mm(p, pMonom, rRing);
425
426 if( v != NULL )
427 {
428 ppPrev = &pNext(p); // fixed!
429
430 // *p is changed if v != NULL ( p == v )
431 pIter(p);
432
433 if( p == NULL )
434 break;
435 }
436 else
437 { // Upps! Zero!!! we must kill this term!
438
439 //
440 p = p_LmDeleteAndNext(p, rRing);
441
442 *ppPrev = p;
443
444 if( p == NULL )
445 break;
446 }
447 }
448
449 #ifdef PDEBUG
450 p_Test(pPoly,rRing);
451 #endif
452
453 return(pPoly);
454 }
455
456 // return new poly = pPoly * pMonom; preserve pPoly and pMonom.
sca_pp_Mult_mm(const poly pPoly,const poly pMonom,const ring rRing)457 poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing)
458 {
459 assume( rIsSCA(rRing) );
460
461 #ifdef PDEBUG
462 // PrintS("sca_pp_Mult_mm\n"); // !
463
464 p_Test(pPoly, rRing);
465 p_Test(pMonom, rRing);
466 #endif
467
468 if (/*(*/ pPoly == NULL /*)*/) /*|| ( pMonom == NULL )*/
469 return NULL;
470
471 const int iComponentMonomM = p_GetComp(pMonom, rRing);
472
473 poly pResult = NULL;
474 poly* ppPrev = &pResult;
475
476 for( poly p = pPoly; p!= NULL; pIter(p) )
477 {
478 #ifdef PDEBUG
479 p_Test(p, rRing);
480 #endif
481 const int iComponent = p_GetComp(p, rRing);
482
483 if( iComponent!=0 )
484 {
485 if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
486 {
487 // REPORT_ERROR
488 Werror("sca_pp_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
489 // what should we do further?!?
490
491 p_Delete( &pResult, rRing); // delete the result
492 return NULL;
493 }
494
495 #ifdef PDEBUG
496 if(iComponentMonomM==0 )
497 {
498 dReportError("sca_pp_Mult_mm: Multiplication in the left module from the right");
499 }
500 #endif
501 }
502
503 // terms will be in the same order because of quasi-ordering!
504 poly v = sca_mm_Mult_mm(p, pMonom, rRing);
505
506 if( v != NULL )
507 {
508 *ppPrev = v;
509 ppPrev = &pNext(v);
510 }
511 }
512
513 #ifdef PDEBUG
514 p_Test(pResult,rRing);
515 #endif
516
517 return(pResult);
518 }
519
520 //-----------------------------------------------------------------------------------//
521
522 // return x_i * pPoly; preserve pPoly.
sca_xi_Mult_pp(short i,const poly pPoly,const ring rRing)523 static inline poly sca_xi_Mult_pp(short i, const poly pPoly, const ring rRing)
524 {
525 assume( rIsSCA(rRing) );
526
527 #ifdef PDEBUG
528 p_Test(pPoly, rRing);
529 #endif
530
531 assume(i <= scaLastAltVar(rRing));
532 assume(scaFirstAltVar(rRing) <= i);
533
534 if( pPoly == NULL )
535 return NULL;
536
537 poly pResult = NULL;
538 poly* ppPrev = &pResult;
539
540 for( poly p = pPoly; p!= NULL; pIter(p) )
541 {
542
543 // terms will be in the same order because of quasi-ordering!
544 poly v = sca_xi_Mult_mm(i, p, rRing);
545
546 #ifdef PDEBUG
547 p_Test(v, rRing);
548 #endif
549
550 if( v != NULL )
551 {
552 *ppPrev = v;
553 ppPrev = &pNext(*ppPrev);
554 }
555 }
556
557
558 #ifdef PDEBUG
559 p_Test(pResult, rRing);
560 #endif
561
562 return(pResult);
563 }
564
565
566 // return new poly = pMonom * pPoly; preserve pPoly and pMonom.
sca_pp_mm_Mult(const poly pPoly,const poly pMonom,const ring rRing)567 static poly sca_pp_mm_Mult(const poly pPoly, const poly pMonom, const ring rRing)
568 {
569 assume( rIsSCA(rRing) );
570
571 #ifdef PDEBUG
572 // PrintS("sca_mm_Mult_pp\n"); // !
573
574 p_Test(pPoly, rRing);
575 p_Test(pMonom, rRing);
576 #endif
577
578 if ((pPoly==NULL) || (pMonom==NULL)) return NULL;
579
580 assume( (pPoly != NULL) && (pMonom !=NULL));
581
582 const int iComponentMonomM = p_GetComp(pMonom, rRing);
583
584 poly pResult = NULL;
585 poly* ppPrev = &pResult;
586
587 for( poly p = pPoly; p!= NULL; pIter(p) )
588 {
589 #ifdef PDEBUG
590 p_Test(p, rRing);
591 #endif
592 const int iComponent = p_GetComp(p, rRing);
593
594 if( iComponentMonomM!=0 )
595 {
596 if( iComponent!=0 ) // TODO: make global if on iComponentMonomM =?= 0
597 {
598 // REPORT_ERROR
599 Werror("sca_pp_mm_Mult: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
600 // what should we do further?!?
601
602 p_Delete( &pResult, rRing); // delete the result
603 return NULL;
604 }
605 #ifdef PDEBUG
606 if(iComponent==0 )
607 {
608 dReportError("sca_pp_mm_Mult: Multiplication in the left module from the right!");
609 // PrintS("mm = "); p_Write(pMonom, rRing);
610 // PrintS("pp = "); p_Write(pPoly, rRing);
611 // assume(iComponent!=0);
612 }
613 #endif
614 }
615
616 // terms will be in the same order because of quasi-ordering!
617 poly v = sca_mm_Mult_mm(pMonom, p, rRing);
618
619 if( v != NULL )
620 {
621 *ppPrev = v;
622 ppPrev = &pNext(*ppPrev); // nice line ;-)
623 }
624 }
625
626 #ifdef PDEBUG
627 p_Test(pResult,rRing);
628 #endif
629
630 return(pResult);
631 }
632
633
634 // return poly = pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
sca_p_mm_Mult(poly pPoly,const poly pMonom,const ring rRing)635 static poly sca_p_mm_Mult(poly pPoly, const poly pMonom, const ring rRing) // !!!!! the MOST used procedure !!!!!
636 {
637 assume( rIsSCA(rRing) );
638
639 #ifdef PDEBUG
640 p_Test(pPoly, rRing);
641 p_Test(pMonom, rRing);
642 #endif
643
644 if( pPoly == NULL )
645 return NULL;
646
647 assume(pMonom!=NULL);
648 //if( pMonom == NULL )
649 //{
650 // // pPoly != NULL =>
651 // p_Delete( &pPoly, rRing );
652 // return NULL;
653 //}
654
655 const int iComponentMonomM = p_GetComp(pMonom, rRing);
656
657 poly p = pPoly; poly* ppPrev = &pPoly;
658
659 loop
660 {
661 #ifdef PDEBUG
662 if( !p_Test(p, rRing) )
663 {
664 PrintS("p is wrong!");
665 p_Write(p,rRing);
666 }
667 #endif
668
669 const int iComponent = p_GetComp(p, rRing);
670
671 if( iComponentMonomM!=0 )
672 {
673 if( iComponent!=0 )
674 {
675 // REPORT_ERROR
676 Werror("sca_p_mm_Mult: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
677 // what should we do further?!?
678
679 p_Delete( &pPoly, rRing); // delete the result
680 return NULL;
681 }
682 #ifdef PDEBUG
683 if(iComponent==0)
684 {
685 dReportError("sca_p_mm_Mult: Multiplication in the left module from the right!");
686 // PrintS("mm = "); p_Write(pMonom, rRing);
687 // PrintS("p = "); p_Write(pPoly, rRing);
688 // assume(iComponent!=0);
689 }
690 #endif
691 }
692
693 // terms will be in the same order because of quasi-ordering!
694 poly v = sca_mm_Mult_m(pMonom, p, rRing);
695
696 if( v != NULL )
697 {
698 ppPrev = &pNext(p);
699
700 // *p is changed if v != NULL ( p == v )
701 pIter(p);
702
703 if( p == NULL )
704 break;
705 }
706 else
707 { // Upps! Zero!!! we must kill this term!
708 p = p_LmDeleteAndNext(p, rRing);
709
710 *ppPrev = p;
711
712 if( p == NULL )
713 break;
714 }
715 }
716
717 #ifdef PDEBUG
718 if( !p_Test(pPoly, rRing) )
719 {
720 PrintS("pPoly is wrong!");
721 p_Write(pPoly, rRing);
722 }
723 #endif
724
725 return(pPoly);
726 }
727
728 //-----------------------------------------------------------------------------------//
729
730 #ifdef PDEBUG
731 #endif
732
733
734
735
736 //-----------------------------------------------------------------------------------//
737
738 // GB computation routines:
739
740 /*4
741 * creates the S-polynomial of p1 and p2
742 * does not destroy p1 and p2
743 */
sca_SPoly(const poly p1,const poly p2,const ring r)744 poly sca_SPoly( const poly p1, const poly p2, const ring r )
745 {
746 assume( rIsSCA(r) );
747
748 const long lCompP1 = p_GetComp(p1,r);
749 const long lCompP2 = p_GetComp(p2,r);
750
751 if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
752 {
753 #ifdef PDEBUG
754 dReportError("sca_SPoly: different non-zero components!\n");
755 #endif
756 return(NULL);
757 }
758
759 poly pL = p_Lcm(p1, p2, r); // pL = lcm( lm(p1), lm(p2) )
760
761 poly m1 = p_One( r);
762 p_ExpVectorDiff(m1, pL, p1, r); // m1 = pL / lm(p1)
763
764 //p_SetComp(m1,0,r);
765 //p_Setm(m1,r);
766 #ifdef PDEBUG
767 p_Test(m1,r);
768 #endif
769
770
771 poly m2 = p_One( r);
772 p_ExpVectorDiff (m2, pL, p2, r); // m2 = pL / lm(p2)
773
774 //p_SetComp(m2,0,r);
775 //p_Setm(m2,r);
776 #ifdef PDEBUG
777 p_Test(m2,r);
778 #endif
779
780 p_Delete(&pL,r);
781
782 number C1 = n_Copy(pGetCoeff(p1),r->cf); // C1 = lc(p1)
783 number C2 = n_Copy(pGetCoeff(p2),r->cf); // C2 = lc(p2)
784
785 number C = n_Gcd(C1,C2,r->cf); // C = gcd(C1, C2)
786
787 if (!n_IsOne(C, r->cf)) // if C != 1
788 {
789 C1=n_Div(C1, C, r->cf); // C1 = C1 / C
790 C2=n_Div(C2, C, r->cf); // C2 = C2 / C
791 }
792
793 n_Delete(&C,r->cf); // destroy the number C
794
795 const int iSignSum = sca_Sign_mm_Mult_mm (m1, p1, r) + sca_Sign_mm_Mult_mm (m2, p2, r);
796 // zero if different signs
797
798 assume( (iSignSum*iSignSum == 0) || (iSignSum*iSignSum == 4) );
799
800 if( iSignSum != 0 ) // the same sign!
801 C2=n_InpNeg (C2, r->cf);
802
803 p_SetCoeff(m1, C2, r); // lc(m1) = C2!!!
804 p_SetCoeff(m2, C1, r); // lc(m2) = C1!!!
805
806 poly tmp1 = nc_mm_Mult_pp (m1, pNext(p1), r); // tmp1 = m1 * tail(p1),
807 p_Delete(&m1,r); // => n_Delete(&C1,r);
808
809 poly tmp2 = nc_mm_Mult_pp (m2, pNext(p2), r); // tmp1 = m2 * tail(p2),
810 p_Delete(&m2,r); // => n_Delete(&C1,r);
811
812 poly spoly = p_Add_q (tmp1, tmp2, r); // spoly = spoly(lt(p1), lt(p2)) + m1 * tail(p1), delete tmp1,2
813
814 if (spoly!=NULL) p_Cleardenom (spoly, r);
815
816 #ifdef PDEBUG
817 p_Test (spoly, r);
818 #endif
819
820 return(spoly);
821 }
822
823
824
825
826 /*2
827 * reduction of p2 with p1
828 * do not destroy p1, but p2
829 * p1 divides p2 -> for use in NF algorithm
830 */
sca_ReduceSpoly(const poly p1,poly p2,const ring r)831 poly sca_ReduceSpoly(const poly p1, poly p2, const ring r)
832 {
833 assume( rIsSCA(r) );
834
835 assume( p1 != NULL );
836
837 const long lCompP1 = p_GetComp (p1, r);
838 const long lCompP2 = p_GetComp (p2, r);
839
840 if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
841 {
842 #ifdef PDEBUG
843 dReportError("sca_ReduceSpoly: different non-zero components!");
844 #endif
845 return(NULL);
846 }
847
848 poly m = p_ISet (1, r);
849 p_ExpVectorDiff (m, p2, p1, r); // m = lm(p2) / lm(p1)
850 //p_Setm(m,r);
851 #ifdef PDEBUG
852 p_Test (m,r);
853 #endif
854
855 number C1 = n_Copy( pGetCoeff(p1), r->cf);
856 number C2 = n_Copy( pGetCoeff(p2), r->cf);
857
858 /* GCD stuff */
859 number C = n_Gcd(C1, C2, r->cf);
860
861 if (!n_IsOne(C, r->cf))
862 {
863 C1 = n_Div(C1, C, r->cf);
864 C2 = n_Div(C2, C, r->cf);
865 }
866 n_Delete(&C,r->cf);
867
868 const int iSign = sca_Sign_mm_Mult_mm( m, p1, r );
869
870 if(iSign == 1)
871 C2 = n_InpNeg(C2,r->cf);
872
873 p_SetCoeff(m, C2, r);
874
875 #ifdef PDEBUG
876 p_Test(m,r);
877 #endif
878
879 p2 = p_LmDeleteAndNext( p2, r );
880
881 p2 = p_Mult_nn(p2, C1, r); // p2 !!!
882 n_Delete(&C1,r->cf);
883
884 poly T = nc_mm_Mult_pp(m, pNext(p1), r);
885 p_Delete(&m, r);
886
887 p2 = p_Add_q(p2, T, r);
888
889 if ( p2!=NULL ) p_Cleardenom(p2,r);
890
891 #ifdef PDEBUG
892 p_Test(p2,r);
893 #endif
894
895 return(p2);
896 }
897
898 // should be used only inside nc_SetupQuotient!
899 // Check whether this our case:
900 // 1. rG is a commutative polynomial ring \otimes anticommutative algebra
901 // 2. factor ideal rGR->qideal contains squares of all alternating variables.
902 //
903 // if yes, make rGR a super-commutative algebra!
904 // NOTE: Factors of SuperCommutative Algebras are supported this way!
905 //
906 // rG == NULL means that there is no separate base G-algebra in this case take rGR == rG
907
908 // special case: bCopy == true (default value: false)
909 // meaning: rGR copies structure from rG
910 // (maybe with some minor changes, which don't change the type!)
sca_SetupQuotient(ring rGR,ring rG,bool bCopy)911 bool sca_SetupQuotient(ring rGR, ring rG, bool bCopy)
912 {
913
914 //////////////////////////////////////////////////////////////////////////
915 // checks...
916 //////////////////////////////////////////////////////////////////////////
917 if( rG == NULL )
918 rG = rGR;
919
920 assume(rGR != NULL);
921 assume(rG != NULL);
922 assume(rIsPluralRing(rG));
923
924 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
925 PrintS("sca_SetupQuotient(rGR, rG, bCopy)");
926
927 {
928 PrintS("\nrG: \n"); rWrite(rG);
929 PrintS("\nrGR: \n"); rWrite(rGR);
930 PrintLn();
931 }
932 #endif
933
934
935 if(bCopy)
936 {
937 if(rIsSCA(rG) && (rG != rGR))
938 return sca_Force(rGR, scaFirstAltVar(rG), scaLastAltVar(rG));
939 else
940 return false;
941 }
942
943 assume(!bCopy);
944
945 const int N = rG->N;
946
947 if(N < 2)
948 return false;
949
950
951 // if( (ncRingType(rG) != nc_skew) || (ncRingType(rG) != nc_comm) )
952 // return false;
953
954 #if OUTPUT
955 PrintS("sca_SetupQuotient: qring?\n");
956 #endif
957
958 if(rGR->qideal == NULL) // there should be a factor!
959 return false;
960
961 #if OUTPUT
962 PrintS("sca_SetupQuotient: qideal!!!\n");
963 #endif
964
965 // if((rG->qideal != NULL) && (rG != rGR) ) // we cannot change from factor to factor at the time, sorry!
966 // return false;
967
968
969 int iAltVarEnd = -1;
970 int iAltVarStart = N+1;
971
972 const nc_struct* NC = rG->GetNC();
973 const ring rBase = rG; //NC->basering;
974 const matrix C = NC->C; // live in rBase!
975 const matrix D = NC->D; // live in rBase!
976
977 #if OUTPUT
978 PrintS("sca_SetupQuotient: AltVars?!\n");
979 #endif
980
981 for(int i = 1; i < N; i++)
982 {
983 for(int j = i + 1; j <= N; j++)
984 {
985 if( MATELEM(D,i,j) != NULL) // !!!???
986 {
987 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
988 Print("Nonzero D[%d, %d]\n", i, j);
989 #endif
990 return false;
991 }
992
993
994 assume(MATELEM(C,i,j) != NULL); // after CallPlural!
995 number c = p_GetCoeff(MATELEM(C,i,j), rBase);
996
997 if( n_IsMOne(c, rBase->cf) ) // !!!???
998 {
999 if( i < iAltVarStart)
1000 iAltVarStart = i;
1001
1002 if( j > iAltVarEnd)
1003 iAltVarEnd = j;
1004 } else
1005 {
1006 if( !n_IsOne(c, rBase->cf) )
1007 {
1008 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1009 Print("Wrong Coeff at: [%d, %d]\n", i, j);
1010 #endif
1011 return false;
1012 }
1013 }
1014 }
1015 }
1016
1017 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1018 Print("AltVars?1: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1019 #endif
1020
1021
1022 if( (iAltVarEnd == -1) || (iAltVarStart == (N+1)) )
1023 return false; // either no alternating varables, or a single one => we are in commutative case!
1024
1025
1026 for(int i = 1; i < N; i++)
1027 {
1028 for(int j = i + 1; j <= N; j++)
1029 {
1030 assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1031 number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1032
1033 if( (iAltVarStart <= i) && (j <= iAltVarEnd) ) // S <= i < j <= E
1034 { // anticommutative part
1035 if( !n_IsMOne(c, rBase->cf) )
1036 {
1037 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1038 Print("Wrong Coeff at: [%d, %d]\n", i, j);
1039 #endif
1040 return false;
1041 }
1042 }
1043 else
1044 { // should commute
1045 if( !n_IsOne(c, rBase->cf) )
1046 {
1047 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1048 Print("Wrong Coeff at: [%d, %d]\n", i, j);
1049 #endif
1050 return false;
1051 }
1052 }
1053 }
1054 }
1055
1056 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1057 Print("AltVars!?: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1058 #endif
1059
1060 assume( 1 <= iAltVarStart );
1061 assume( iAltVarStart < iAltVarEnd );
1062 assume( iAltVarEnd <= N );
1063
1064
1065 // ring rSaveRing = assureCurrentRing(rG);
1066
1067
1068 assume(rGR->qideal != NULL);
1069 assume(rGR->N == rG->N);
1070 // assume(rG->qideal == NULL); // ?
1071
1072 const ideal idQuotient = rGR->qideal;
1073
1074
1075 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1076 PrintS("Analyzing quotient ideal:\n");
1077 idPrint(idQuotient); // in rG!!!
1078 #endif
1079
1080
1081 // check for
1082 // y_{iAltVarStart}^2, y_{iAltVarStart+1}^2, \ldots, y_{iAltVarEnd}^2 (iAltVarEnd > iAltVarStart)
1083 // to be within quotient ideal.
1084
1085 int b = N+1;
1086 int e = -1;
1087
1088 if(rIsSCA(rG))
1089 {
1090 b = si_min(b, scaFirstAltVar(rG));
1091 e = si_max(e, scaLastAltVar(rG));
1092
1093 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1094 Print("AltVars!?: [%d, %d]\n", b, e);
1095 #endif
1096 }
1097
1098 for ( int i = iAltVarStart; (i <= iAltVarEnd); i++ )
1099 if( (i < b) || (i > e) ) // otherwise it's ok since rG is an SCA!
1100 {
1101 poly square = p_One( rG);
1102 p_SetExp(square, i, 2, rG); // square = var(i)^2.
1103 p_Setm(square, rG);
1104
1105 // square = NF( var(i)^2 | Q )
1106 // NOTE: there is no better way to check this in general!
1107 square = nc_NF(idQuotient, NULL, square, 0, 1, rG); // must ran in currRing == rG!
1108
1109 if( square != NULL ) // var(i)^2 is not in Q?
1110 {
1111 p_Delete(&square, rG);
1112 return false;
1113 }
1114 }
1115
1116 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1117 Print("ScaVars!: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1118 #endif
1119
1120
1121 //////////////////////////////////////////////////////////////////////////
1122 // ok... here we go. let's setup it!!!
1123 //////////////////////////////////////////////////////////////////////////
1124 ideal tempQ = id_KillSquares(idQuotient, iAltVarStart, iAltVarEnd, rG); // in rG!!!
1125
1126
1127 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1128 PrintS("Quotient: \n");
1129 iiWriteMatrix((matrix)idQuotient,"__",1, rG, 0);
1130 PrintS("tempSCAQuotient: \n");
1131 iiWriteMatrix((matrix)tempQ,"__",1, rG, 0);
1132 #endif
1133
1134 idSkipZeroes( tempQ );
1135
1136 ncRingType( rGR, nc_exterior );
1137
1138 scaFirstAltVar( rGR, iAltVarStart );
1139 scaLastAltVar( rGR, iAltVarEnd );
1140
1141 if( idIs0(tempQ) )
1142 rGR->GetNC()->SCAQuotient() = NULL;
1143 else
1144 rGR->GetNC()->SCAQuotient() = idrMoveR(tempQ, rG, rGR); // deletes tempQ!
1145
1146 nc_p_ProcsSet(rGR, rGR->p_Procs); // !!!!!!!!!!!!!!!!!
1147
1148
1149 #if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1150 PrintS("SCAQuotient: \n");
1151 if(tempQ != NULL)
1152 iiWriteMatrix((matrix)tempQ,"__",1, rGR, 0);
1153 else
1154 PrintS("(NULL)\n");
1155 #endif
1156
1157 return true;
1158 }
1159
1160
sca_Force(ring rGR,int b,int e)1161 bool sca_Force(ring rGR, int b, int e)
1162 {
1163 assume(rGR != NULL);
1164 assume(rIsPluralRing(rGR));
1165 assume(!rIsSCA(rGR));
1166
1167 const int N = rGR->N;
1168
1169 // ring rSaveRing = currRing;
1170 // if(rSaveRing != rGR)
1171 // rChangeCurrRing(rGR);
1172
1173 const ideal idQuotient = rGR->qideal;
1174
1175 ideal tempQ = idQuotient;
1176
1177 if( b <= N && e >= 1 )
1178 tempQ = id_KillSquares(idQuotient, b, e, rGR);
1179
1180 idSkipZeroes( tempQ );
1181
1182 ncRingType( rGR, nc_exterior );
1183
1184 if( idIs0(tempQ) )
1185 rGR->GetNC()->SCAQuotient() = NULL;
1186 else
1187 rGR->GetNC()->SCAQuotient() = tempQ;
1188
1189
1190 scaFirstAltVar( rGR, b );
1191 scaLastAltVar( rGR, e );
1192
1193
1194 nc_p_ProcsSet(rGR, rGR->p_Procs);
1195
1196 // if(rSaveRing != rGR)
1197 // rChangeCurrRing(rSaveRing);
1198
1199 return true;
1200 }
1201
1202 // return x_i * pPoly; preserve pPoly.
sca_pp_Mult_xi_pp(short i,const poly pPoly,const ring rRing)1203 poly sca_pp_Mult_xi_pp(short i, const poly pPoly, const ring rRing)
1204 {
1205 assume(1 <= i);
1206 assume(i <= rVar(rRing));
1207
1208 if(rIsSCA(rRing))
1209 return sca_xi_Mult_pp(i, pPoly, rRing);
1210
1211
1212
1213 poly xi = p_One( rRing);
1214 p_SetExp(xi, i, 1, rRing);
1215 p_Setm(xi, rRing);
1216
1217 poly pResult = pp_Mult_qq(xi, pPoly, rRing);
1218
1219 p_Delete( &xi, rRing);
1220
1221 return pResult;
1222
1223 }
1224
sca_p_ProcsSet(ring rGR,p_Procs_s * p_Procs)1225 void sca_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
1226 {
1227
1228 // "commutative" procedures:
1229 rGR->p_Procs->p_Mult_mm = sca_p_Mult_mm;
1230 rGR->p_Procs->pp_Mult_mm = sca_pp_Mult_mm;
1231
1232 p_Procs->p_Mult_mm = sca_p_Mult_mm;
1233 p_Procs->pp_Mult_mm = sca_pp_Mult_mm;
1234
1235 // non-commutaitve
1236 p_Procs->p_mm_Mult = sca_p_mm_Mult;
1237 p_Procs->pp_mm_Mult = sca_pp_mm_Mult;
1238
1239 // rGR->GetNC()->p_Procs.SPoly = sca_SPoly;
1240 // rGR->GetNC()->p_Procs.ReduceSPoly = sca_ReduceSpoly;
1241
1242 #if 0
1243
1244 // Multiplication procedures:
1245
1246 p_Procs->p_Mult_mm = sca_p_Mult_mm;
1247 _p_procs->p_Mult_mm = sca_p_Mult_mm;
1248
1249 p_Procs->pp_Mult_mm = sca_pp_Mult_mm;
1250 _p_procs->pp_Mult_mm = sca_pp_Mult_mm;
1251
1252 r->GetNC()->mmMultP() = sca_mm_Mult_p;
1253 r->GetNC()->mmMultPP() = sca_mm_Mult_pp;
1254
1255 /*
1256 // ??????????????????????????????????????? coefficients swell...
1257 r->GetNC()->SPoly() = sca_SPoly;
1258 r->GetNC()->ReduceSPoly() = sca_ReduceSpoly;
1259 */
1260 // r->GetNC()->BucketPolyRed() = gnc_kBucketPolyRed;
1261 // r->GetNC()->BucketPolyRed_Z()= gnc_kBucketPolyRed_Z;
1262 #endif
1263
1264 // local ordering => Mora, otherwise - Buchberger!
1265 if (rHasLocalOrMixedOrdering(rGR))
1266 rGR->GetNC()->p_Procs.GB = cast_A_to_vptr(sca_mora);
1267 else
1268 rGR->GetNC()->p_Procs.GB = cast_A_to_vptr(sca_bba); // sca_gr_bba?
1269 }
1270
1271
1272 // bi-Degree (x, y) of monomial "m"
1273 // x-es and y-s are weighted by wx and wy resp.
1274 // [optional] components have weights by wCx and wCy.
m_GetBiDegree(const poly m,const intvec * wx,const intvec * wy,const intvec * wCx,const intvec * wCy,int & dx,int & dy,const ring r)1275 static inline void m_GetBiDegree(const poly m,
1276 const intvec *wx, const intvec *wy,
1277 const intvec *wCx, const intvec *wCy,
1278 int& dx, int& dy, const ring r)
1279 {
1280 const unsigned int N = r->N;
1281
1282 p_Test(m, r);
1283
1284 assume( wx != NULL );
1285 assume( wy != NULL );
1286
1287 assume( wx->cols() == 1 );
1288 assume( wy->cols() == 1 );
1289
1290 assume( (unsigned int)wx->rows() >= N );
1291 assume( (unsigned int)wy->rows() >= N );
1292
1293 int x = 0;
1294 int y = 0;
1295
1296 for(int i = N; i > 0; i--)
1297 {
1298 const int d = p_GetExp(m, i, r);
1299 x += d * (*wx)[i-1];
1300 y += d * (*wy)[i-1];
1301 }
1302
1303 if( (wCx != NULL) && (wCy != NULL) )
1304 {
1305 const int c = p_GetComp(m, r);
1306
1307 if( wCx->range(c) )
1308 x += (*wCx)[c];
1309
1310 if( wCy->range(c) )
1311 x += (*wCy)[c];
1312 }
1313
1314 dx = x;
1315 dy = y;
1316 }
1317
1318 // returns true if polynom p is bi-homogenous with respect to the given weights
1319 // simultaneously sets bi-Degree
p_IsBiHomogeneous(const poly p,const intvec * wx,const intvec * wy,const intvec * wCx,const intvec * wCy,int & dx,int & dy,const ring r)1320 bool p_IsBiHomogeneous(const poly p,
1321 const intvec *wx, const intvec *wy,
1322 const intvec *wCx, const intvec *wCy,
1323 int &dx, int &dy,
1324 const ring r)
1325 {
1326 if( p == NULL )
1327 {
1328 dx = 0;
1329 dy = 0;
1330 return true;
1331 }
1332
1333 poly q = p;
1334
1335
1336 int ddx, ddy;
1337
1338 m_GetBiDegree( q, wx, wy, wCx, wCy, ddx, ddy, r); // get bi degree of lm(p)
1339
1340 pIter(q);
1341
1342 for(; q != NULL; pIter(q) )
1343 {
1344 int x, y;
1345
1346 m_GetBiDegree( q, wx, wy, wCx, wCy, x, y, r); // get bi degree of q
1347
1348 if ( (x != ddx) || (y != ddy) ) return false;
1349 }
1350
1351 dx = ddx;
1352 dy = ddy;
1353
1354 return true;
1355 }
1356
1357
1358 // returns true if id is bi-homogenous without respect to the given weights
id_IsBiHomogeneous(const ideal id,const intvec * wx,const intvec * wy,const intvec * wCx,const intvec * wCy,const ring r)1359 bool id_IsBiHomogeneous(const ideal id,
1360 const intvec *wx, const intvec *wy,
1361 const intvec *wCx, const intvec *wCy,
1362 const ring r)
1363 {
1364 if (id == NULL) return true; // zero ideal
1365
1366 const int iSize = IDELEMS(id);
1367
1368 if (iSize == 0) return true;
1369
1370 bool b = true;
1371 int x, y;
1372
1373 for(int i = iSize - 1; (i >= 0 ) && b; i--)
1374 b = p_IsBiHomogeneous(id->m[i], wx, wy, wCx, wCy, x, y, r);
1375
1376 return b;
1377 }
1378
1379
1380 // returns an intvector with [nvars(r)] integers [1/0]
1381 // 1 - for commutative variables
1382 // 0 - for anticommutative variables
ivGetSCAXVarWeights(const ring r)1383 intvec *ivGetSCAXVarWeights(const ring r)
1384 {
1385 const unsigned int N = r->N;
1386
1387 const int CommutativeVariable = 0; // bug correction!
1388 const int AntiCommutativeVariable = 0;
1389
1390 intvec* w = new intvec(N, 1, CommutativeVariable);
1391
1392 if(AntiCommutativeVariable != CommutativeVariable)
1393 if( rIsSCA(r) )
1394 {
1395 const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
1396 const unsigned int m_iLastAltVar = scaLastAltVar(r);
1397
1398 for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
1399 {
1400 (*w)[i-1] = AntiCommutativeVariable;
1401 }
1402 }
1403
1404 return w;
1405 }
1406
1407
1408 // returns an intvector with [nvars(r)] integers [1/0]
1409 // 0 - for commutative variables
1410 // 1 - for anticommutative variables
ivGetSCAYVarWeights(const ring r)1411 intvec *ivGetSCAYVarWeights(const ring r)
1412 {
1413 const unsigned int N = r->N;
1414
1415 const int CommutativeVariable = 0;
1416 const int AntiCommutativeVariable = 1;
1417
1418 intvec* w = new intvec(N, 1, CommutativeVariable);
1419
1420 if(AntiCommutativeVariable != CommutativeVariable)
1421 if( rIsSCA(r) )
1422 {
1423 const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
1424 const unsigned int m_iLastAltVar = scaLastAltVar(r);
1425
1426 for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
1427 {
1428 (*w)[i-1] = AntiCommutativeVariable;
1429 }
1430 }
1431 return w;
1432 }
1433
1434
1435
1436
1437 // reduce term lt(m) modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar:
1438 // either create a copy of m or return NULL
m_KillSquares(const poly m,const short iFirstAltVar,const short iLastAltVar,const ring r)1439 static inline poly m_KillSquares(const poly m,
1440 const short iFirstAltVar, const short iLastAltVar,
1441 const ring r)
1442 {
1443 #ifdef PDEBUG
1444 p_Test(m, r);
1445 assume( (iFirstAltVar >= 1) && (iLastAltVar <= rVar(r)) && (iFirstAltVar <= iLastAltVar) );
1446
1447 #if 0
1448 PrintS("m_KillSquares, m = "); // !
1449 p_Write(m, r);
1450 #endif
1451 #endif
1452
1453 assume( m != NULL );
1454
1455 for(short k = iFirstAltVar; k <= iLastAltVar; k++)
1456 if( p_GetExp(m, k, r) > 1 )
1457 return NULL;
1458
1459 return p_Head(m, r);
1460 }
1461
1462
1463 // reduce polynomial p modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
1464 // returns a new poly!
p_KillSquares(const poly p,const short iFirstAltVar,const short iLastAltVar,const ring r)1465 poly p_KillSquares(const poly p,
1466 const short iFirstAltVar, const short iLastAltVar,
1467 const ring r)
1468 {
1469 #ifdef PDEBUG
1470 p_Test(p, r);
1471
1472 assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
1473
1474 #if 0
1475 PrintS("p_KillSquares, p = "); // !
1476 p_Write(p, r);
1477 #endif
1478 #endif
1479
1480
1481 if( p == NULL )
1482 return NULL;
1483
1484 poly pResult = NULL;
1485 poly* ppPrev = &pResult;
1486
1487 for( poly q = p; q!= NULL; pIter(q) )
1488 {
1489 #ifdef PDEBUG
1490 p_Test(q, r);
1491 #endif
1492
1493 // terms will be in the same order because of quasi-ordering!
1494 poly v = m_KillSquares(q, iFirstAltVar, iLastAltVar, r);
1495
1496 if( v != NULL )
1497 {
1498 *ppPrev = v;
1499 ppPrev = &pNext(v);
1500 }
1501
1502 }
1503
1504 #ifdef PDEBUG
1505 p_Test(pResult, r);
1506 #if 0
1507 PrintS("p_KillSquares => "); // !
1508 p_Write(pResult, r);
1509 #endif
1510 #endif
1511
1512 return(pResult);
1513 }
1514
1515
1516
1517
1518 // reduces ideal id modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
1519 // returns the reduced ideal or zero ideal.
id_KillSquares(const ideal id,const short iFirstAltVar,const short iLastAltVar,const ring r,const bool bSkipZeroes)1520 ideal id_KillSquares(const ideal id,
1521 const short iFirstAltVar, const short iLastAltVar,
1522 const ring r, const bool bSkipZeroes)
1523 {
1524 if (id == NULL) return id; // zero ideal
1525
1526 assume( (iFirstAltVar >= 1) && (iLastAltVar <= rVar(r)) && (iFirstAltVar <= iLastAltVar) );
1527
1528 const int iSize = IDELEMS(id);
1529
1530 if (iSize == 0) return id;
1531
1532 ideal temp = idInit(iSize, id->rank);
1533
1534 #if 0
1535 PrintS("<id_KillSquares>\n");
1536 {
1537 PrintS("ideal id: \n");
1538 for (unsigned int i = 0; i < IDELEMS(id); i++)
1539 {
1540 Print("; id[%d] = ", i+1);
1541 p_Write(id->m[i], r);
1542 }
1543 PrintS(";\n");
1544 PrintLn();
1545 }
1546 #endif
1547
1548
1549 for (int j = 0; j < iSize; j++)
1550 temp->m[j] = p_KillSquares(id->m[j], iFirstAltVar, iLastAltVar, r);
1551
1552 if( bSkipZeroes )
1553 idSkipZeroes(temp);
1554
1555 #if 0
1556 PrintS("<id_KillSquares>\n");
1557 {
1558 PrintS("ideal temp: \n");
1559 for (int i = 0; i < IDELEMS(temp); i++)
1560 {
1561 Print("; temp[%d] = ", i+1);
1562 p_Write(temp->m[i], r);
1563 }
1564 PrintS(";\n");
1565 PrintLn();
1566 }
1567 PrintS("</id_KillSquares>\n");
1568 #endif
1569
1570 // temp->rank = idRankFreeModule(temp, r);
1571
1572 return temp;
1573 }
1574
1575
1576
1577
1578 #endif
1579