source: git/kernel/sca.cc @ 129878

spielwiese
Last change on this file since 129878 was 129878, checked in by Motsak Oleksandr <motsak@…>, 17 years ago
! better SCAQuotient initializing git-svn-id: file:///usr/local/Singular/svn/trunk@9892 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 56.6 KB
Line 
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 *  Version: $Id: sca.cc,v 1.10 2007-02-23 14:41:43 motsak Exp $
10 *******************************************************************/
11
12// #define PDEBUG 2
13#include "mod2.h"
14
15#ifdef HAVE_PLURAL
16// for
17#define PLURAL_INTERNAL_DECLARATIONS
18#include "sca.h"
19
20
21#include "febase.h"
22
23#include "p_polys.h"
24#include "kutil.h"
25#include "ideals.h"
26#include "intvec.h"
27#include "polys.h"
28
29#include "ring.h"
30#include "numbers.h"
31#include "matpol.h"
32#include "kbuckets.h"
33#include "kstd1.h"
34#include "sbuckets.h"
35#include "prCopy.h"
36#include "p_Mult_q.h"
37#include "p_MemAdd.h"
38
39#include "kutil.h"
40#include "kstd1.h"
41
42#include "weight.h"
43
44
45// poly functions defined in p_Procs :
46
47// return pPoly * pMonom; preserve pPoly and pMonom.
48poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing, poly &);
49
50// return pMonom * pPoly; preserve pPoly and pMonom.
51poly sca_mm_Mult_pp(const poly pMonom, const poly pPoly, const ring rRing);
52
53// return pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
54poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing);
55
56// return pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
57poly sca_mm_Mult_p(const poly pMonom, poly pPoly, const ring rRing);
58
59
60// compute the spolynomial of p1 and p2
61poly sca_SPoly(const poly p1, const poly p2, const ring r);
62poly sca_ReduceSpoly(const poly p1, poly p2, const ring r);
63
64// Modified Plural's Buchberger's algorithmus.
65ideal sca_gr_bba(const ideal F, const ideal Q, const intvec *, const intvec *, kStrategy strat);
66
67// Modified modern Sinuglar Buchberger's algorithm.
68ideal sca_bba(const ideal F, const ideal Q, const intvec *w, const intvec *, kStrategy strat);
69
70// Modified modern Sinuglar Mora's algorithm.
71ideal sca_mora(const ideal F, const ideal Q, const intvec *w, const intvec *, kStrategy strat);
72
73
74
75////////////////////////////////////////////////////////////////////////////////////////////////////
76// Super Commutative Algabra extension by Motsak
77////////////////////////////////////////////////////////////////////////////////////////////////////
78
79
80
81
82
83// returns the sign of: lm(pMonomM) * lm(pMonomMM),
84// preserves input, may return +/-1, 0
85inline int sca_Sign_mm_Mult_mm( const poly pMonomM, const poly pMonomMM, const ring rRing )
86{
87#ifdef PDEBUG
88    p_Test(pMonomM,  rRing);
89    p_Test(pMonomMM, rRing);
90#endif
91
92    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
93    const unsigned int iLastAltVar  = scaLastAltVar(rRing);
94
95    unsigned int tpower = 0;
96    unsigned int cpower = 0;
97
98    for( unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
99    {
100      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
101      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
102
103      if( iExpMM != 0 )
104      {
105        if( iExpM != 0 )
106        {
107          return 0; // lm(pMonomM) * lm(pMonomMM) == 0
108        }
109        tpower += cpower; // compute degree of (-1).
110      }
111
112      cpower += iExpM;
113    }
114
115    if( (tpower&1) != 0 ) // degree is odd => negate coeff.
116      return -1;
117
118    return(1);
119}
120
121
122
123
124// returns and changes pMonomM: lt(pMonomM) = lt(pMonomM) * lt(pMonomMM),
125// preserves pMonomMM. may return NULL!
126// if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomM) * lt(pMonomMM)
127// if result == NULL => pMonomM MUST BE deleted manually!
128inline poly sca_m_Mult_mm( poly pMonomM, const poly pMonomMM, const ring rRing )
129{
130#ifdef PDEBUG
131    p_Test(pMonomM,  rRing);
132    p_Test(pMonomMM, rRing);
133#endif
134
135    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
136    const unsigned int iLastAltVar = scaLastAltVar(rRing);
137
138    unsigned int tpower = 0;
139    unsigned int cpower = 0;
140
141    for( unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
142    {
143      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
144      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
145
146      if( iExpMM != 0 )
147      {
148        if( iExpM != 0 ) // result is zero!
149        {
150          return NULL; // we do nothing with pMonomM in this case!
151        }
152
153        tpower += cpower; // compute degree of (-1).
154      }
155
156      cpower += iExpM;
157    }
158
159    p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
160
161    number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
162
163    if( (tpower&1) != 0 ) // degree is odd => negate coeff.
164      nCoeffM = n_Neg(nCoeffM, rRing); // negate nCoeff (will destroy the original number)
165
166    const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
167
168    number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing); // new number!
169
170    p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
171
172#ifdef PDEBUG
173    p_Test(pMonomM, rRing);
174#endif
175
176    return(pMonomM);
177}
178
179// returns and changes pMonomM: lt(pMonomM) = lt(pMonomMM) * lt(pMonomM),
180// preserves pMonomMM. may return NULL!
181// if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomMM) * lt(pMonomM)
182// if result == NULL => pMonomM MUST BE deleted manually!
183inline poly sca_mm_Mult_m( const poly pMonomMM, poly pMonomM, const ring rRing )
184{
185#ifdef PDEBUG
186    p_Test(pMonomM,  rRing);
187    p_Test(pMonomMM, rRing);
188#endif
189
190    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
191    const unsigned int iLastAltVar = scaLastAltVar(rRing);
192
193    unsigned int tpower = 0;
194    unsigned int cpower = 0;
195
196    for( unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
197    {
198      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
199      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
200
201      if( iExpM != 0 )
202      {
203        if( iExpMM != 0 ) // result is zero!
204        {
205          return NULL; // we do nothing with pMonomM in this case!
206        }
207
208        tpower += cpower; // compute degree of (-1).
209      }
210
211      cpower += iExpMM;
212    }
213
214    p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
215
216    number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
217
218    if( (tpower&1) != 0 ) // degree is odd => negate coeff.
219      nCoeffM = n_Neg(nCoeffM, rRing); // negate nCoeff (will destroy the original number), creates new number!
220
221    const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
222
223    number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing); // new number!
224
225    p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
226
227#ifdef PDEBUG
228    p_Test(pMonomM, rRing);
229#endif
230
231    return(pMonomM);
232}
233
234
235
236// returns: result = lt(pMonom1) * lt(pMonom2),
237// preserves pMonom1, pMonom2. may return NULL!
238// if result != NULL => next(result) = NULL, lt(result) = lt(pMonom1) * lt(pMonom2)
239inline poly sca_mm_Mult_mm( poly pMonom1, const poly pMonom2, const ring rRing )
240{
241#ifdef PDEBUG
242    p_Test(pMonom1, rRing);
243    p_Test(pMonom2, rRing);
244#endif
245
246    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
247    const unsigned int iLastAltVar = scaLastAltVar(rRing);
248
249    unsigned int tpower = 0;
250    unsigned int cpower = 0;
251
252    for( unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
253    {
254      const unsigned int iExp1 = p_GetExp(pMonom1, j, rRing);
255      const unsigned int iExp2 = p_GetExp(pMonom2, j, rRing);
256
257      if( iExp2 != 0 )
258      {
259        if( iExp1 != 0 ) // result is zero!
260        {
261          return NULL;
262        }
263
264        tpower += cpower; // compute degree of (-1).
265      }
266
267      cpower += iExp1;
268    }
269
270    poly pResult;
271    omTypeAllocBin(poly, pResult, rRing->PolyBin);
272    p_SetRingOfLm(pResult, rRing);
273
274    p_ExpVectorSum(pResult, pMonom1, pMonom2, rRing); // "exponents" are additive!!!
275
276    pNext(pResult) = NULL;
277
278    const number nCoeff1 = p_GetCoeff(pMonom1, rRing); // no new copy!
279    const number nCoeff2 = p_GetCoeff(pMonom2, rRing); // no new copy!
280
281    number nCoeff = n_Mult(nCoeff1, nCoeff2, rRing); // new number!
282
283    if( (tpower&1) != 0 ) // degree is odd => negate coeff.
284      nCoeff = n_Neg(nCoeff, rRing); // negate nCoeff (will destroy the original number)
285
286    p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
287
288#ifdef PDEBUG
289    p_Test(pResult, rRing);
290#endif
291
292    return(pResult);
293}
294
295// returns: result =  x_i * lt(pMonom),
296// preserves pMonom. may return NULL!
297// if result != NULL => next(result) = NULL, lt(result) = x_i * lt(pMonom)
298inline poly sca_xi_Mult_mm(unsigned int i, const poly pMonom, const ring rRing)
299{
300#ifdef PDEBUG
301    p_Test(pMonom, rRing);
302#endif
303
304    assume( i <= scaLastAltVar(rRing));
305    assume( scaFirstAltVar(rRing) <= i );
306
307    if( p_GetExp(pMonom, i, rRing) != 0 ) // => result is zero!
308      return NULL;
309
310    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
311
312    unsigned int cpower = 0;
313
314    for( unsigned int j = iFirstAltVar; j < i ; j++ )
315      cpower += p_GetExp(pMonom, j, rRing);
316
317    poly pResult = p_LmInit(pMonom, rRing);
318
319    p_SetExp(pResult, i, 1, rRing); // pResult*=X_i &&
320    p_Setm(pResult, rRing);         // addjust degree after previous step!
321
322    number nCoeff = n_Copy(p_GetCoeff(pMonom, rRing), rRing); // new number!
323
324    if( (cpower&1) != 0 ) // degree is odd => negate coeff.
325      nCoeff = n_Neg(nCoeff, rRing); // negate nCoeff (will destroy the original number)
326
327    p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
328
329#ifdef PDEBUG
330    p_Test(pResult, rRing);
331#endif
332
333    return(pResult);
334}
335
336//-----------------------------------------------------------------------------------//
337
338// return poly = pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
339poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing)
340{
341  assume( rIsSCA(rRing) );
342
343#ifdef PDEBUG
344//  Print("sca_p_Mult_mm\n"); // !
345
346  p_Test(pPoly, rRing);
347  p_Test(pMonom, rRing);
348#endif
349
350  if( pPoly == NULL )
351    return NULL;
352
353  if( pMonom == NULL )
354  {
355    // pPoly != NULL =>
356    p_Delete( &pPoly, rRing );
357
358    return NULL;
359  }
360
361  const int iComponentMonomM = p_GetComp(pMonom, rRing);
362
363  poly p = pPoly; poly* ppPrev = &pPoly;
364
365  loop
366  {
367#ifdef PDEBUG
368    p_Test(p, rRing);
369#endif
370    const int iComponent = p_GetComp(p, rRing);
371
372    if( iComponent!=0 )
373    {
374      if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
375      {
376        // REPORT_ERROR
377        Werror("sca_p_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
378        // what should we do further?!?
379
380        p_Delete( &pPoly, rRing); // delete the result AND rest
381        return NULL;
382      }
383#ifdef PDEBUG
384      if(iComponentMonomM==0 )
385      {
386        Print("Multiplication in the left module from the right");
387      }
388#endif
389    }
390
391    // terms will be in the same order because of quasi-ordering!
392    poly v = sca_m_Mult_mm(p, pMonom, rRing);
393
394    if( v != NULL )
395    {
396      ppPrev = &pNext(p); // fixed!
397
398    // *p is changed if v != NULL ( p == v )
399      pIter(p);
400
401      if( p == NULL )
402        break;
403    }
404    else
405    { // Upps! Zero!!! we must kill this term!
406
407      //
408      p = p_LmDeleteAndNext(p, rRing);
409
410      *ppPrev = p;
411
412      if( p == NULL )
413        break;
414    }
415
416
417  }
418
419#ifdef PDEBUG
420  p_Test(pPoly,rRing);
421#endif
422
423  return(pPoly);
424}
425
426// return new poly = pPoly * pMonom; preserve pPoly and pMonom.
427poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing, poly &)
428{
429  assume( rIsSCA(rRing) );
430
431#ifdef PDEBUG
432//  Print("sca_pp_Mult_mm\n"); // !
433
434  p_Test(pPoly, rRing);
435  p_Test(pMonom, rRing);
436#endif
437
438  if( ( pPoly == NULL ) || ( pMonom == NULL ) )
439    return NULL;
440
441  const int iComponentMonomM = p_GetComp(pMonom, rRing);
442
443  poly pResult = NULL;
444  poly* ppPrev = &pResult;
445
446  for( poly p = pPoly; p!= NULL; pIter(p) )
447  {
448#ifdef PDEBUG
449    p_Test(p, rRing);
450#endif
451    const int iComponent = p_GetComp(p, rRing);
452
453    if( iComponent!=0 )
454    {
455      if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
456      {
457        // REPORT_ERROR
458        Werror("sca_pp_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
459        // what should we do further?!?
460
461        p_Delete( &pResult, rRing); // delete the result
462        return NULL;
463      }
464
465#ifdef PDEBUG
466      if(iComponentMonomM==0 )
467      {
468        Print("Multiplication in the left module from the right");
469      }
470#endif
471    }
472
473    // terms will be in the same order because of quasi-ordering!
474    poly v = sca_mm_Mult_mm(p, pMonom, rRing);
475
476    if( v != NULL )
477    {
478      *ppPrev = v;
479      ppPrev = &pNext(v);
480    }
481
482  }
483
484#ifdef PDEBUG
485  p_Test(pResult,rRing);
486#endif
487
488  return(pResult);
489}
490
491//-----------------------------------------------------------------------------------//
492
493// return x_i * pPoly; preserve pPoly.
494inline poly sca_xi_Mult_pp(unsigned int i, const poly pPoly, const ring rRing)
495{
496  assume( rIsSCA(rRing) );
497
498#ifdef PDEBUG
499  p_Test(pPoly, rRing);
500#endif
501
502  assume(i <= scaLastAltVar(rRing));
503  assume(scaFirstAltVar(rRing) <= i);
504
505  if( pPoly == NULL )
506    return NULL;
507
508  poly pResult = NULL;
509  poly* ppPrev = &pResult;
510
511  for( poly p = pPoly; p!= NULL; pIter(p) )
512  {
513
514    // terms will be in the same order because of quasi-ordering!
515    poly v = sca_xi_Mult_mm(i, p, rRing);
516
517#ifdef PDEBUG
518    p_Test(v, rRing);
519#endif
520
521    if( v != NULL )
522    {
523      *ppPrev = v;
524      ppPrev = &pNext(*ppPrev);
525    }
526  }
527
528
529#ifdef PDEBUG
530  p_Test(pResult, rRing);
531#endif
532
533  return(pResult);
534}
535
536
537// return new poly = pMonom * pPoly; preserve pPoly and pMonom.
538poly sca_mm_Mult_pp(const poly pMonom, const poly pPoly, const ring rRing)
539{
540  assume( rIsSCA(rRing) );
541
542#ifdef PDEBUG
543//  Print("sca_mm_Mult_pp\n"); // !
544
545  p_Test(pPoly, rRing);
546  p_Test(pMonom, rRing);
547#endif
548
549  if( ( pPoly == NULL ) || ( pMonom == NULL ) )
550    return NULL;
551
552  const int iComponentMonomM = p_GetComp(pMonom, rRing);
553
554  poly pResult = NULL;
555  poly* ppPrev = &pResult;
556
557  for( poly p = pPoly; p!= NULL; pIter(p) )
558  {
559#ifdef PDEBUG
560    p_Test(p, rRing);
561#endif
562    const int iComponent = p_GetComp(p, rRing);
563
564    if( iComponentMonomM!=0 )
565    {
566      if( iComponent!=0 ) // TODO: make global if on iComponentMonomM =?= 0
567      {
568        // REPORT_ERROR
569        Werror("sca_mm_Mult_pp: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
570        // what should we do further?!?
571
572        p_Delete( &pResult, rRing); // delete the result
573        return NULL;
574      }
575#ifdef PDEBUG
576      if(iComponent==0 )
577      {
578        Print("Multiplication in the left module from the right");
579      }
580#endif
581    }
582
583    // terms will be in the same order because of quasi-ordering!
584    poly v = sca_mm_Mult_mm(pMonom, p, rRing);
585
586    if( v != NULL )
587    {
588      *ppPrev = v;
589      ppPrev = &pNext(*ppPrev); // nice line ;-)
590    }
591  }
592
593#ifdef PDEBUG
594  p_Test(pResult,rRing);
595#endif
596
597  return(pResult);
598}
599
600
601// return poly = pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
602poly sca_mm_Mult_p(const poly pMonom, poly pPoly, const ring rRing) // !!!!! the MOST used procedure !!!!!
603{
604  assume( rIsSCA(rRing) );
605
606#ifdef PDEBUG
607  p_Test(pPoly, rRing);
608  p_Test(pMonom, rRing);
609#endif
610
611  if( pPoly == NULL )
612    return NULL;
613
614  if( pMonom == NULL )
615  {
616    // pPoly != NULL =>
617    p_Delete( &pPoly, rRing );
618    return NULL;
619  }
620
621  const int iComponentMonomM = p_GetComp(pMonom, rRing);
622
623  poly p = pPoly; poly* ppPrev = &pPoly;
624
625  loop
626  {
627#ifdef PDEBUG
628    if( !p_Test(p, rRing) )
629    {
630      Print("p is wrong!");
631      p_Write(p,rRing);
632    }
633#endif
634
635    const int iComponent = p_GetComp(p, rRing);
636
637    if( iComponentMonomM!=0 )
638    {
639      if( iComponent!=0 ) // TODO: make global if on iComponentMonomM =?= 0
640      {
641        // REPORT_ERROR
642        Werror("sca_mm_Mult_p: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
643        // what should we do further?!?
644
645        p_Delete( &pPoly, rRing); // delete the result
646        return NULL;
647      }
648#ifdef PDEBUG
649      if(iComponent==0 )
650      {
651        Print("Multiplication in the left module from the right");
652      }
653#endif
654    }
655
656    // terms will be in the same order because of quasi-ordering!
657    poly v = sca_mm_Mult_m(pMonom, p, rRing);
658
659    if( v != NULL )
660    {
661      ppPrev = &pNext(p);
662
663    // *p is changed if v != NULL ( p == v )
664      pIter(p);
665
666      if( p == NULL )
667        break;
668    }
669    else
670    { // Upps! Zero!!! we must kill this term!
671      p = p_LmDeleteAndNext(p, rRing);
672
673      *ppPrev = p;
674
675      if( p == NULL )
676        break;
677    }
678  }
679
680#ifdef PDEBUG
681    if( !p_Test(pPoly, rRing) )
682    {
683      Print("pPoly is wrong!");
684      p_Write(pPoly, rRing);
685    }
686#endif
687
688  return(pPoly);
689}
690
691//-----------------------------------------------------------------------------------//
692
693#ifdef PDEBUG
694#endif
695
696
697
698
699//-----------------------------------------------------------------------------------//
700
701// GB computation routines:
702
703/*2
704* returns the LCM of the head terms of a and b
705*/
706inline poly p_Lcm(const poly a, const poly b, const long lCompM, const ring r)
707{
708  poly m = p_ISet(1, r);
709
710  const int pVariables = r->N;
711
712  for (int i = pVariables; i; i--)
713  {
714    const int lExpA = p_GetExp (a, i, r);
715    const int lExpB = p_GetExp (b, i, r);
716
717    p_SetExp (m, i, si_max(lExpA, lExpB), r);
718  }
719
720  p_SetComp (m, lCompM, r);
721
722  p_Setm(m,r);
723
724#ifdef PDEBUG
725  p_Test(m,r);
726#endif
727
728  return(m);
729}
730
731/*4
732* creates the S-polynomial of p1 and p2
733* does not destroy p1 and p2
734*/
735poly sca_SPoly( const poly p1, const poly p2, const ring r )
736{
737  assume( rIsSCA(r) );
738
739  const long lCompP1 = p_GetComp(p1,r);
740  const long lCompP2 = p_GetComp(p2,r);
741
742  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
743  {
744#ifdef PDEBUG
745    Print("sca_SPoly: different non-zero components!\n");
746#endif
747    return(NULL);
748  }
749
750  poly pL = p_Lcm(p1, p2, si_max(lCompP1, lCompP2), r);       // pL = lcm( lm(p1), lm(p2) )
751
752  poly m1 = p_ISet(1, r);
753  p_ExpVectorDiff(m1, pL, p1, r);                  // m1 = pL / lm(p1)
754  //p_SetComp(m1,0,r);
755  //p_Setm(m1,r);
756#ifdef PDEBUG
757  p_Test(m1,r);
758#endif
759
760
761  poly m2 = p_ISet(1, r);
762  p_ExpVectorDiff (m2, pL, p2, r);                  // m2 = pL / lm(p2)
763
764  //p_SetComp(m2,0,r);
765  //p_Setm(m2,r);
766#ifdef PDEBUG
767  p_Test(m2,r);
768#endif
769
770  p_Delete(&pL,r);
771
772  number C1  = n_Copy(p_GetCoeff(p1,r),r);      // C1 = lc(p1)
773  number C2  = n_Copy(p_GetCoeff(p2,r),r);      // C2 = lc(p2)
774
775  number C = nGcd(C1,C2,r);                     // C = gcd(C1, C2)
776
777  if (!n_IsOne(C, r))                              // if C != 1
778  {
779    C1=n_Div(C1, C, r);                              // C1 = C1 / C
780    C2=n_Div(C2, C, r);                              // C2 = C2 / C
781  }
782
783  n_Delete(&C,r); // destroy the number C
784
785  const int iSignSum = sca_Sign_mm_Mult_mm (m1, p1, r) + sca_Sign_mm_Mult_mm (m2, p2, r);
786  // zero if different signs
787
788  assume( (iSignSum*iSignSum == 0) || (iSignSum*iSignSum == 4) );
789
790  if( iSignSum != 0 ) // the same sign!
791    C2=n_Neg (C2, r);
792
793  p_SetCoeff(m1, C2, r);                           // lc(m1) = C2!!!
794  p_SetCoeff(m2, C1, r);                           // lc(m2) = C1!!!
795
796  poly tmp1 = nc_mm_Mult_pp (m1, pNext(p1), r);    // tmp1 = m1 * tail(p1),
797  p_Delete(&m1,r);  //  => n_Delete(&C1,r);
798
799  poly tmp2 = nc_mm_Mult_pp (m2, pNext(p2), r);    // tmp1 = m2 * tail(p2),
800  p_Delete(&m2,r);  //  => n_Delete(&C1,r);
801
802  poly spoly = p_Add_q (tmp1, tmp2, r); // spoly = spoly(lt(p1), lt(p2)) + m1 * tail(p1), delete tmp1,2
803
804  if (spoly!=NULL) pCleardenom (spoly); // r?
805//  if (spoly!=NULL) pContent (spoly); // r?
806
807#ifdef PDEBUG
808  p_Test (spoly, r);
809#endif
810
811  return(spoly);
812}
813
814
815
816
817/*2
818* reduction of p2 with p1
819* do not destroy p1, but p2
820* p1 divides p2 -> for use in NF algorithm
821*/
822poly sca_ReduceSpoly(const poly p1, poly p2, const ring r)
823{
824  assume( rIsSCA(r) );
825
826  assume( p1 != NULL );
827
828  const long lCompP1 = p_GetComp (p1, r);
829  const long lCompP2 = p_GetComp (p2, r);
830
831  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
832  {
833#ifdef PDEBUG
834    Print("sca_ReduceSpoly: different non-zero components!");
835#endif
836    return(NULL);
837  }
838
839  poly m = p_ISet (1, r);
840  p_ExpVectorDiff (m, p2, p1, r);                      // m = lm(p2) / lm(p1)
841  //p_Setm(m,r);
842#ifdef PDEBUG
843  p_Test (m,r);
844#endif
845
846  number C1 = n_Copy( p_GetCoeff(p1, r), r);
847  number C2 = n_Copy( p_GetCoeff(p2, r), r);
848
849  /* GCD stuff */
850  number C = nGcd(C1, C2, r);
851
852  if (!n_IsOne(C, r))
853  {
854    C1 = n_Div(C1, C, r);
855    C2 = n_Div(C2, C, r);
856  }
857  n_Delete(&C,r);
858
859  const int iSign = sca_Sign_mm_Mult_mm( m, p1, r );
860
861  if(iSign == 1)
862    C2 = n_Neg(C2,r);
863
864  p_SetCoeff(m, C2, r);
865
866#ifdef PDEBUG
867  p_Test(m,r);
868#endif
869
870  p2 = p_LmDeleteAndNext( p2, r );
871
872  p2 = p_Mult_nn(p2, C1, r); // p2 !!!
873  n_Delete(&C1,r);
874
875  poly T = nc_mm_Mult_pp(m, pNext(p1), r);
876  p_Delete(&m, r);
877
878  p2 = p_Add_q(p2, T, r);
879
880  if ( p2!=NULL ) pContent(p2); // r?
881
882#ifdef PDEBUG
883  p_Test(p2,r);
884#endif
885
886  return(p2);
887}
888
889
890void addLObject(LObject& h, kStrategy& strat)
891{
892  if(h.IsNull()) return;
893
894  strat->initEcart(&h);
895  h.sev=0; // pGetShortExpVector(h.p);
896
897  // add h into S and L
898  int pos=posInS(strat, strat->sl, h.p, h.ecart);
899
900  if ( (pos <= strat->sl) && (pComparePolys(h.p, strat->S[pos])) )
901  {
902    if (TEST_OPT_PROT)
903      PrintS("d\n");
904  }
905  else
906  {
907    if (TEST_OPT_INTSTRATEGY)
908    {
909      pCleardenom(h.p);
910    }
911    else
912    {
913      pNorm(h.p);
914      pContent(h.p);
915    }
916
917    if ((strat->syzComp==0)||(!strat->homog))
918    {
919      h.p = redtailBba(h.p,pos-1,strat);
920
921      if (TEST_OPT_INTSTRATEGY)
922      {
923//        pCleardenom(h.p);
924        pContent(h.p);
925      }
926      else
927      {
928        pNorm(h.p);
929      }
930    }
931
932    if(h.IsNull()) return;
933
934    /* statistic */
935    if (TEST_OPT_PROT)
936    {
937      PrintS("s\n");
938    }
939
940    if (TEST_OPT_DEBUG)
941    {
942      PrintS("new s:");
943      wrp(h.p);
944      PrintLn();
945    }
946
947    enterpairs(h.p, strat->sl, h.ecart, 0, strat);
948
949    pos=0;
950
951    if (strat->sl!=-1) pos = posInS(strat, strat->sl, h.p, h.ecart);
952
953    strat->enterS(h, pos, strat, -1);
954
955    if (h.lcm!=NULL) pLmFree(h.lcm);
956  }
957
958
959}
960
961
962ideal sca_gr_bba(const ideal F, const ideal Q, const intvec *, const intvec *, kStrategy strat)
963{
964  assume(rIsSCA(currRing));
965
966  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
967  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
968 
969  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
970  const ideal tempQ = currRing->nc->SCAQuotient();
971
972  bool bIdHomog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
973
974  assume( !bIdHomog || strat->homog ); //  bIdHomog =====[implies]>>>>> strat->homog
975
976  strat->homog = strat->homog && bIdHomog;
977
978#ifdef PDEBUG
979  assume( strat->homog == bIdHomog );
980#endif /*PDEBUG*/
981
982
983#if 0
984   PrintS("<sca::bba>\n");
985  {
986    Print("ideal F: \n");
987    for (int i = 0; i < tempF->idelems(); i++)
988    {
989      Print("; F[%d] = ", i+1);
990      p_Write(tempF->m[i], currRing);
991    }
992    Print(";\n");
993    PrintLn();
994  }
995#endif
996
997
998#ifdef PDEBUG
999//   PrintS("sca_gr_bba\n");
1000#endif /*PDEBUG*/
1001
1002
1003
1004  int srmax, lrmax;
1005  int olddeg, reduc;
1006  int red_result = 1;
1007  int hilbeledeg = 1, hilbcount = 0, minimcnt = 0;
1008
1009  initBuchMoraCrit(strat); // set Gebauer, honey, sugarCrit
1010
1011  nc_gr_initBba(tempF,strat); // set enterS, red, initEcart, initEcartPair
1012
1013  initBuchMoraPos(strat);
1014
1015
1016  // ?? set spSpolyShort, reduce ???
1017
1018  initBuchMora(tempF, tempQ, strat); // currRing->nc->SCAQuotient() instead of Q == squares!!!!!!!
1019
1020  strat->posInT=posInT110; // !!!
1021
1022  srmax = strat->sl;
1023  reduc = olddeg = lrmax = 0;
1024
1025
1026  /* compute------------------------------------------------------- */
1027  for(; strat->Ll >= 0;
1028#ifdef KDEBUG
1029    strat->P.lcm = NULL,
1030#endif
1031    kTest(strat)
1032    )
1033  {
1034    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
1035
1036    if (TEST_OPT_DEBUG) messageSets(strat);
1037
1038    if (strat->Ll== 0) strat->interpt=TRUE;
1039
1040    if (TEST_OPT_DEGBOUND
1041    && ((strat->honey
1042    && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1043       || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
1044    {
1045      /*
1046      *stops computation if
1047      * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
1048      *a predefined number Kstd1_deg
1049      */
1050      while (strat->Ll >= 0) deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1051      break;
1052    }
1053
1054    /* picks the last element from the lazyset L */
1055    strat->P = strat->L[strat->Ll];
1056    strat->Ll--;
1057
1058    //kTest(strat);
1059
1060    assume(pNext(strat->P.p) != strat->tail);
1061
1062//     if( pNext(strat->P.p) == strat->tail )
1063//     {
1064//       // deletes the int spoly and computes SPoly
1065//       pLmFree(strat->P.p); // ???
1066//       strat->P.p = sca_SPoly(strat->P.p1, strat->P.p2, currRing);
1067//     }
1068
1069    if(strat->P.IsNull()) continue;
1070
1071//     poly save = NULL;
1072//
1073//     if(pNext(strat->P.p) != NULL)
1074//       save = p_Copy(strat->P.p, currRing);
1075
1076    strat->initEcart(&strat->P); // remove it?
1077
1078    if (TEST_OPT_PROT)
1079      message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(), &olddeg,&reduc,strat, red_result);
1080
1081    // reduction of the element chosen from L wrt S
1082    strat->red(&strat->P,strat);
1083
1084    if(strat->P.IsNull()) continue;
1085
1086    addLObject(strat->P, strat);
1087
1088    if (strat->sl > srmax) srmax = strat->sl;
1089
1090    const poly save = strat->P.p;
1091
1092#ifdef PDEBUG
1093      p_Test(save, currRing);
1094#endif
1095    assume( save != NULL );
1096   
1097    // SCA Specials:
1098
1099    {
1100      const poly pNext = pNext(save);
1101     
1102      if( pNext != NULL )
1103      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
1104      if( p_GetExp(save, i, currRing) != 0 )
1105      {
1106        assume(p_GetExp(save, i, currRing) == 1);
1107       
1108        const poly tt = sca_pp_Mult_xi_pp(i, pNext, currRing);
1109
1110#ifdef PDEBUG
1111        p_Test(tt, currRing);
1112#endif
1113
1114        if( tt == NULL) continue;
1115
1116        LObject h(tt); // h = x_i * P
1117
1118        if (TEST_OPT_INTSTRATEGY)
1119        {
1120//           h.pCleardenom(); // also does a pContent
1121          pContent(h.p);
1122        }
1123        else
1124        {
1125          h.pNorm();
1126        }
1127
1128        strat->initEcart(&h);
1129
1130
1131//         if (pOrdSgn==-1)
1132//         {
1133//           cancelunit(&h);  // tries to cancel a unit
1134//           deleteHC(&h, strat);
1135//         }
1136
1137//         if(h.IsNull()) continue;
1138
1139//         if (TEST_OPT_PROT)
1140//           message((strat->honey ? h.ecart : 0) + h.pFDeg(), &olddeg, &reduc, strat, red_result);
1141
1142//         strat->red(&h, strat); // wrt S
1143//         if(h.IsNull()) continue;
1144
1145//         poly save = p_Copy(h.p, currRing);
1146
1147        int pos;
1148
1149        if (strat->Ll==-1)
1150          pos =0;
1151        else
1152          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1153
1154        h.sev = pGetShortExpVector(h.p);
1155        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
1156
1157  //       h.p = save;
1158  //       addLObject(h, strat);
1159
1160//         if (strat->sl > srmax) srmax = strat->sl;
1161      }
1162
1163      // p_Delete( &save, currRing );
1164    }
1165
1166
1167  } // for(;;)
1168
1169
1170  if (TEST_OPT_DEBUG) messageSets(strat);
1171
1172  /* release temp data-------------------------------- */
1173  exitBuchMora(strat);
1174
1175  if (TEST_OPT_WEIGHTM)
1176  {
1177    pFDeg=pFDegOld;
1178    pLDeg=pLDegOld;
1179    if (ecartWeights)
1180    {
1181      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(int));
1182      ecartWeights=NULL;
1183    }
1184  }
1185
1186  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
1187
1188  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
1189 
1190  id_Delete(&tempF, currRing);
1191
1192
1193  /* complete reduction of the standard basis--------- */
1194  if (TEST_OPT_REDSB){
1195//     completeReduce(strat); // ???
1196
1197    ideal I = strat->Shdl;
1198    ideal erg = kInterRed(I,tempQ);
1199    assume(I!=erg);
1200    id_Delete(&I, currRing);
1201    strat->Shdl = erg;
1202  }
1203
1204
1205#if 0
1206   PrintS("</sca::bba>\n");
1207#endif
1208
1209  return (strat->Shdl);
1210}
1211
1212
1213// should be used only inside nc_SetupQuotient!
1214// Check whether this our case:
1215//  1. rG is  a commutative polynomial ring \otimes anticommutative algebra
1216//  2. factor ideal rGR->qideal contains squares of all alternating variables.
1217//
1218// if yes, make rGR a super-commutative algebra!
1219// NOTE: Factors of SuperCommutative Algebras are supported this way!
1220bool sca_SetupQuotient(ring rGR, const ring rG)
1221{
1222//   return false; // test Plural
1223
1224  //////////////////////////////////////////////////////////////////////////
1225  // checks...
1226  //////////////////////////////////////////////////////////////////////////
1227  assume(rGR != NULL);
1228  assume(rG  != NULL);
1229
1230  assume(rIsPluralRing(rG));
1231
1232  const int N = rG->N;
1233
1234  if(N < 2)
1235    return false;
1236
1237  if(ncRingType(rG) != nc_skew)
1238    return false;
1239
1240  if(rGR->qideal == NULL) // there will be a factor!
1241    return false;
1242
1243  if(rG->qideal != NULL) // we cannot change from factor to factor!
1244    return false;
1245
1246
1247  int iAltVarEnd = -1;
1248  int iAltVarStart   = N+1;
1249
1250  const ring rBase = rG->nc->basering;
1251  const matrix C   = rG->nc->C; // live in rBase!
1252
1253  for(int i = 1; i < N; i++)
1254  {
1255    for(int j = i + 1; j <= N; j++)
1256    {
1257      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1258      number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1259     
1260      if( n_IsMOne(c, rBase) )
1261      {       
1262        if( i < iAltVarStart) 
1263          iAltVarStart = i;
1264         
1265        if( j > iAltVarEnd)
1266          iAltVarEnd = j;
1267      } else
1268      {
1269        if( !n_IsOne(c, rBase) )
1270        {
1271#ifdef PDEBUG
1272//           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1273#endif
1274          return false;
1275        }
1276      }
1277    }
1278  }
1279
1280  if( (iAltVarEnd == -1) || (iAltVarStart == (N+1)) )
1281    return false; // either no alternating varables, or a single one => we are in commutative case!
1282 
1283  for(int i = 1; i < N; i++)
1284  {
1285    for(int j = i + 1; j <= N; j++)
1286    {
1287      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1288      number c = p_GetCoeff(MATELEM(C,i,j), rBase);     
1289
1290      if( (iAltVarStart <= i) && (j <= iAltVarEnd) ) // S <= i < j <= E
1291      { // anticommutative part
1292        if( !n_IsMOne(c, rBase) )
1293        {
1294#ifdef PDEBUG
1295//           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1296#endif
1297          return false;
1298        }
1299      } else
1300      { // should commute
1301        if( !n_IsOne(c, rBase) )
1302        {
1303#ifdef PDEBUG
1304//           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1305#endif
1306          return false;
1307        }
1308      }
1309    }
1310  }
1311
1312  assume( 1            <= iAltVarStart );
1313  assume( iAltVarStart < iAltVarEnd   );
1314  assume( iAltVarEnd   <= N            );
1315
1316  bool bWeChangeRing = false;
1317  // for sanity
1318  ring rSaveRing = currRing;
1319
1320  if(rSaveRing != rG)
1321  {
1322    rChangeCurrRing(rG);
1323    bWeChangeRing = true;
1324  }
1325
1326  assume(rGR->qideal != NULL);
1327  assume(rG->qideal == NULL);
1328 
1329  const ideal idQuotient = rGR->qideal;
1330
1331  // check for
1332  // y_{iAltVarStart}^2, y_{iAltVarStart+1}^2, \ldots, y_{iAltVarEnd}^2  (iAltVarEnd > iAltVarStart)
1333  // to be within quotient ideal.
1334
1335  bool bSCA = true;
1336
1337  for ( int i = iAltVarStart; (i <= iAltVarEnd) && bSCA; i++ )
1338  {
1339    poly square = p_ISet(1, rSaveRing);
1340    p_SetExp(square, i, 2, rSaveRing); // square = var(i)^2.
1341    p_Setm(square, rSaveRing);
1342
1343    // square = NF( var(i)^2 | Q )
1344    // NOTE: rSaveRing == currRing now!
1345    // NOTE: there is no better way to check this in general!
1346    square = kNF(idQuotient, NULL, square, 0, 0); 
1347   
1348    if( square != NULL ) // var(i)^2 is not in Q?
1349    {
1350      p_Delete(&square, rSaveRing);
1351      bSCA = false;     
1352    }   
1353  }
1354 
1355
1356  if (bWeChangeRing)
1357  {
1358    rChangeCurrRing(rSaveRing);
1359  }
1360
1361  if(!bSCA) return false;
1362 
1363
1364
1365#ifdef PDEBUG
1366//  Print("AltVars: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1367#endif
1368
1369
1370  //////////////////////////////////////////////////////////////////////////
1371  // ok... here we go. let's setup it!!!
1372  //////////////////////////////////////////////////////////////////////////
1373  ideal tempQ = id_KillSquares(idQuotient, iAltVarStart, iAltVarEnd, rG); // in rG!!!
1374
1375  idSkipZeroes( tempQ );
1376 
1377  if( idIs0(tempQ) )
1378    rGR->nc->SCAQuotient() = NULL;
1379  else
1380    rGR->nc->SCAQuotient() = idrMoveR(tempQ, rG, rGR); // deletes tempQ!
1381
1382  ncRingType( rGR, nc_exterior );
1383
1384  scaFirstAltVar( rGR, iAltVarStart );
1385  scaLastAltVar( rGR, iAltVarEnd );
1386
1387
1388
1389  sca_p_ProcsSet(rGR, rGR->p_Procs);
1390
1391 
1392  return true;
1393}
1394
1395// return x_i * pPoly; preserve pPoly.
1396poly sca_pp_Mult_xi_pp(unsigned int i, const poly pPoly, const ring rRing)
1397{
1398  assume(1 <= i);
1399  assume(i <= rRing->N);
1400
1401  if(rIsSCA(rRing))
1402    return sca_xi_Mult_pp(i, pPoly, rRing);
1403
1404
1405
1406  poly xi =  p_ISet(1, rRing);
1407  p_SetExp(xi, i, 1, rRing);
1408  p_Setm(xi, rRing);
1409
1410  poly pResult = pp_Mult_qq(xi, pPoly, rRing);
1411
1412  p_Delete( &xi, rRing);
1413
1414  return pResult;
1415
1416}
1417
1418
1419///////////////////////////////////////////////////////////////////////////////////////
1420//************* SCA BBA *************************************************************//
1421///////////////////////////////////////////////////////////////////////////////////////
1422
1423// Under development!!!
1424ideal sca_bba (const ideal F, const ideal Q, const intvec *w, const intvec *hilb, kStrategy strat)
1425{
1426  assume(rIsSCA(currRing));
1427
1428  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
1429  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
1430
1431  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
1432
1433  bool bIdHomog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
1434
1435  assume( !bIdHomog || strat->homog ); //  bIdHomog =====[implies]>>>>> strat->homog
1436
1437  strat->homog = strat->homog && bIdHomog;
1438
1439#ifdef PDEBUG
1440  assume( strat->homog == bIdHomog );
1441#endif /*PDEBUG*/
1442
1443
1444//    PrintS("<sca>\n");
1445
1446  om_Opts.MinTrack = 5; // ???
1447
1448  int   srmax, lrmax, red_result = 1;
1449  int   olddeg, reduc;
1450  int hilbeledeg=1, hilbcount=0, minimcnt=0;
1451
1452  BOOLEAN withT = FALSE;
1453
1454  initBuchMoraCrit(strat); // sets Gebauer, honey, sugarCrit // sca - ok???
1455  initBuchMoraPos(strat); // sets strat->posInL, strat->posInT // check!! (Plural's: )
1456
1457//   initHilbCrit(F, Q, &hilb, strat);
1458
1459//  nc_gr_initBba(F,strat);
1460  initBba(tempF, strat); // set enterS, red, initEcart, initEcartPair
1461
1462  /*set enterS, spSpolyShort, reduce, red, initEcart, initEcartPair*/
1463  // ?? set spSpolyShort, reduce ???
1464  initBuchMora(tempF, NULL, strat); // Q = squares!!!
1465
1466//   if (strat->minim>0) strat->M = idInit(IDELEMS(F),F->rank);
1467
1468  srmax = strat->sl;
1469  reduc = olddeg = lrmax = 0;
1470
1471#define NO_BUCKETS
1472
1473#ifndef NO_BUCKETS
1474  if (!TEST_OPT_NOT_BUCKETS)
1475    strat->use_buckets = 1;
1476#endif
1477
1478  // redtailBBa against T for inhomogenous input
1479  if (!K_TEST_OPT_OLDSTD)
1480    withT = ! strat->homog;
1481
1482  // strat->posInT = posInT_pLength;
1483  kTest_TS(strat);
1484
1485#undef HAVE_TAIL_RING
1486
1487#ifdef HAVE_TAIL_RING
1488  kStratInitChangeTailRing(strat);
1489#endif
1490
1491  ///////////////////////////////////////////////////////////////
1492  // SCA:
1493
1494  /* compute------------------------------------------------------- */
1495  while (strat->Ll >= 0)
1496  {
1497    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
1498
1499#ifdef KDEBUG
1500//     loop_count++;
1501    if (TEST_OPT_DEBUG) messageSets(strat);
1502#endif
1503
1504    if (strat->Ll== 0) strat->interpt=TRUE;
1505
1506    if (TEST_OPT_DEGBOUND
1507        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1508            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
1509    {
1510      /*
1511       *stops computation if
1512       * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
1513       *a predefined number Kstd1_deg
1514       */
1515      while ((strat->Ll >= 0)
1516  && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
1517        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1518            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))
1519  )
1520        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1521      if (strat->Ll<0) break;
1522      else strat->noClearS=TRUE;
1523    }
1524
1525    /* picks the last element from the lazyset L */
1526    strat->P = strat->L[strat->Ll];
1527    strat->Ll--;
1528
1529
1530    assume(pNext(strat->P.p) != strat->tail);
1531
1532/*    if (pNext(strat->P.p) == strat->tail)
1533    {
1534      // deletes the short spoly
1535      pLmFree(strat->P.p);
1536      strat->P.p = NULL;
1537      poly m1 = NULL, m2 = NULL;
1538
1539      // check that spoly creation is ok
1540      while (strat->tailRing != currRing &&
1541             !kCheckSpolyCreation(&(strat->P), strat, m1, m2))
1542      {
1543        assume(m1 == NULL && m2 == NULL);
1544        // if not, change to a ring where exponents are at least
1545        // large enough
1546        kStratChangeTailRing(strat);
1547      }
1548
1549#ifdef PDEBUG
1550      Print("ksCreateSpoly!#?");
1551#endif
1552
1553      // create the real one
1554      ksCreateSpoly(&(strat->P), NULL, strat->use_buckets,
1555                    strat->tailRing, m1, m2, strat->R); //?????????
1556    }
1557    else */
1558
1559    if (strat->P.p1 == NULL)
1560    {
1561//       if (strat->minim > 0)
1562//         strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
1563
1564
1565      // for input polys, prepare reduction
1566      strat->P.PrepareRed(strat->use_buckets);
1567    }
1568
1569    if (TEST_OPT_PROT)
1570      message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
1571              &olddeg,&reduc,strat, red_result);
1572
1573    /* reduction of the element choosen from L */
1574    red_result = strat->red(&strat->P,strat);
1575
1576
1577    // reduction to non-zero new poly
1578    if (red_result == 1)
1579    {
1580      /* statistic */
1581      if (TEST_OPT_PROT) PrintS("s");
1582
1583      // get the polynomial (canonicalize bucket, make sure P.p is set)
1584      strat->P.GetP(strat->lmBin);
1585
1586      int pos = posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
1587
1588      // reduce the tail and normalize poly
1589      if (TEST_OPT_INTSTRATEGY)
1590      {
1591        strat->P.pCleardenom();
1592        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1593        {
1594          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT); // !!!
1595          strat->P.pCleardenom();
1596
1597        }
1598      }
1599      else
1600      {
1601
1602        strat->P.pNorm();
1603        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1604          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
1605      }
1606
1607
1608#ifdef KDEBUG
1609      if (TEST_OPT_DEBUG){PrintS(" ns:");strat->P.wrp();PrintLn();}
1610#endif
1611
1612//       // min_std stuff
1613//       if ((strat->P.p1==NULL) && (strat->minim>0))
1614//       {
1615//         if (strat->minim==1)
1616//         {
1617//           strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
1618//           p_Delete(&strat->P.p2, currRing, strat->tailRing);
1619//         }
1620//         else
1621//         {
1622//           strat->M->m[minimcnt]=strat->P.p2;
1623//           strat->P.p2=NULL;
1624//         }
1625//         if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
1626//           pNext(strat->M->m[minimcnt])
1627//             = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
1628//                                            strat->tailRing, currRing,
1629//                                            currRing->PolyBin);
1630//         minimcnt++;
1631//       }
1632
1633      // enter into S, L, and T
1634      if(withT)
1635        enterT(strat->P, strat);
1636
1637      // L
1638      enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
1639
1640      // posInS only depends on the leading term
1641      strat->enterS(strat->P, pos, strat, strat->tl);
1642
1643//       if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
1644
1645//      Print("[%d]",hilbeledeg);
1646      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
1647
1648      if (strat->sl>srmax) srmax = strat->sl;
1649
1650      // //////////////////////////////////////////////////////////
1651      // SCA:
1652      const poly pSave = strat->P.p;
1653      const poly pNext = pNext(pSave);
1654
1655//       if(0)
1656      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
1657      if( p_GetExp(pSave, i, currRing) != 0 )
1658      {
1659        assume(p_GetExp(pSave, i, currRing) == 1);
1660       
1661        const poly pNew = sca_pp_Mult_xi_pp(i, pNext, currRing);
1662
1663#ifdef PDEBUG
1664        p_Test(pNew, currRing);
1665#endif
1666
1667        if( pNew == NULL) continue;
1668
1669        LObject h(pNew); // h = x_i * strat->P
1670
1671        if (TEST_OPT_INTSTRATEGY)
1672        {
1673//           h.pCleardenom(); // also does a pContent
1674          pContent(h.p);
1675        }
1676        else
1677        {
1678          h.pNorm();
1679        }
1680
1681        strat->initEcart(&h);
1682        h.sev = pGetShortExpVector(h.p);
1683
1684        int pos;
1685
1686        if (strat->Ll==-1)
1687          pos =0;
1688        else
1689          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1690
1691        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
1692
1693
1694/*
1695        h.sev = pGetShortExpVector(h.p);
1696        strat->initEcart(&h);
1697
1698        h.PrepareRed(strat->use_buckets);
1699
1700        // reduction of the element choosen from L(?)
1701        red_result = strat->red(&h,strat);
1702
1703        // reduction to non-zero new poly
1704        if (red_result != 1) continue;
1705
1706
1707        int pos = posInS(strat,strat->sl,h.p,h.ecart);
1708
1709        // reduce the tail and normalize poly
1710        if (TEST_OPT_INTSTRATEGY)
1711        {
1712          h.pCleardenom();
1713          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1714          {
1715            h.p = redtailBba(&(h),pos-1,strat, withT); // !!!
1716            h.pCleardenom();
1717          }
1718        }
1719        else
1720        {
1721
1722          h.pNorm();
1723          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1724            h.p = redtailBba(&(h),pos-1,strat, withT);
1725        }
1726
1727
1728  #ifdef KDEBUG
1729        if (TEST_OPT_DEBUG){PrintS(" N:");h.wrp();PrintLn();}
1730  #endif
1731
1732//         h.PrepareRed(strat->use_buckets); // ???
1733
1734        h.sev = pGetShortExpVector(h.p);
1735        strat->initEcart(&h);
1736
1737        if (strat->Ll==-1)
1738          pos = 0;
1739        else
1740          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1741
1742         enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);*/
1743
1744      } // for all x_i \in Ann(lm(P))
1745
1746
1747
1748
1749
1750    } // if red(P) != NULL
1751
1752//     else if (strat->P.p1 == NULL && strat->minim > 0)
1753//     {
1754//       p_Delete(&strat->P.p2, currRing, strat->tailRing);
1755//     }
1756
1757// #ifdef KDEBUG
1758    memset(&(strat->P), 0, sizeof(strat->P));
1759// #endif
1760
1761    kTest_TS(strat);
1762
1763//     Print("\n$\n");
1764
1765  }
1766
1767#ifdef KDEBUG
1768  if (TEST_OPT_DEBUG) messageSets(strat);
1769#endif
1770
1771  /* complete reduction of the standard basis--------- */
1772  if (TEST_OPT_REDSB) completeReduce(strat);
1773
1774  /* release temp data-------------------------------- */
1775  id_Delete(&tempF, currRing);
1776
1777  exitBuchMora(strat);
1778
1779  if (TEST_OPT_WEIGHTM)
1780  {
1781    pRestoreDegProcs(pFDegOld, pLDegOld);
1782    if (ecartWeights)
1783    {
1784      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
1785      ecartWeights=NULL;
1786    }
1787  }
1788
1789  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
1790
1791//   if (Q!=NULL) updateResult(strat->Shdl,Q,strat);
1792
1793//    PrintS("</sca>\n");
1794
1795  return (strat->Shdl);
1796}
1797
1798
1799// //////////////////////////////////////////////////////////////////////////////
1800// sca mora...
1801
1802// returns TRUE if mora should use buckets, false otherwise
1803static BOOLEAN kMoraUseBucket(kStrategy strat)
1804{
1805#ifdef MORA_USE_BUCKETS
1806  if (TEST_OPT_NOT_BUCKETS)
1807    return FALSE;
1808  if (strat->red == redFirst)
1809  {
1810#ifdef NO_LDEG
1811    if (!strat->syzComp)
1812      return TRUE;
1813#else
1814    if ((strat->homog || strat->honey) && !strat->syzComp)
1815      return TRUE;
1816#endif
1817  }
1818  else
1819  {
1820    assume(strat->red == redEcart);
1821    if (strat->honey && !strat->syzComp)
1822      return TRUE;
1823  }
1824#endif
1825  return FALSE;
1826}
1827
1828
1829#ifdef HAVE_ASSUME
1830static int sca_mora_count = 0;
1831static int sca_mora_loop_count;
1832#endif
1833
1834// ideal sca_mora (ideal F, ideal Q, intvec *w, intvec *, kStrategy strat)
1835ideal sca_mora(const ideal F, const ideal Q, const intvec *w, const intvec *, kStrategy strat)
1836{
1837  assume(rIsSCA(currRing));
1838
1839  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
1840  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
1841 
1842  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
1843 
1844  const ideal tempQ = currRing->nc->SCAQuotient();
1845
1846  bool bIdHomog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
1847
1848  assume( !bIdHomog || strat->homog ); //  bIdHomog =====[implies]>>>>> strat->homog
1849
1850  strat->homog = strat->homog && bIdHomog;
1851
1852#ifdef PDEBUG
1853  assume( strat->homog == bIdHomog );
1854#endif /*PDEBUG*/
1855
1856#ifdef HAVE_ASSUME
1857  sca_mora_count++;
1858  sca_mora_loop_count = 0;
1859#endif
1860
1861#ifdef KDEBUG
1862  om_Opts.MinTrack = 5;
1863#endif
1864  int srmax;
1865  int lrmax = 0;
1866  int olddeg = 0;
1867  int reduc = 0;
1868  int red_result = 1;
1869  int hilbeledeg=1,hilbcount=0;
1870
1871  strat->update = TRUE;
1872  /*- setting global variables ------------------- -*/
1873  initBuchMoraCrit(strat);
1874//   initHilbCrit(F,NULL,&hilb,strat); // no Q!
1875  initMora(tempF, strat);
1876  initBuchMoraPos(strat);
1877  /*Shdl=*/initBuchMora(tempF, tempQ, strat); // temp Q, F!
1878//   if (TEST_OPT_FASTHC) missingAxis(&strat->lastAxis,strat);
1879  /*updateS in initBuchMora has Hecketest
1880  * and could have put strat->kHEdgdeFound FALSE*/
1881#if 0
1882  if (ppNoether!=NULL)
1883  {
1884    strat->kHEdgeFound = TRUE;
1885  }
1886  if (strat->kHEdgeFound && strat->update)
1887  {
1888    firstUpdate(strat);
1889    updateLHC(strat);
1890    reorderL(strat);
1891  }
1892  if (TEST_OPT_FASTHC && (strat->lastAxis) && strat->posInLOldFlag)
1893  {
1894    strat->posInLOld = strat->posInL;
1895    strat->posInLOldFlag = FALSE;
1896    strat->posInL = posInL10;
1897    updateL(strat);
1898    reorderL(strat);
1899  }
1900#endif
1901
1902  srmax = strat->sl;
1903  kTest_TS(strat);
1904
1905  strat->use_buckets = kMoraUseBucket(strat);
1906  /*- compute-------------------------------------------*/
1907
1908#undef HAVE_TAIL_RING
1909
1910#ifdef HAVE_TAIL_RING
1911//  if (strat->homog && strat->red == redFirst)
1912//     kStratInitChangeTailRing(strat);
1913#endif
1914
1915
1916  while (strat->Ll >= 0)
1917  {
1918#ifdef HAVE_ASSUME
1919    sca_mora_loop_count++;
1920#endif
1921    if (lrmax< strat->Ll) lrmax=strat->Ll; /*stat*/
1922    //test_int_std(strat->kIdeal);
1923    if (TEST_OPT_DEBUG) messageSets(strat);
1924    if (TEST_OPT_DEGBOUND
1925    && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg))
1926    {
1927      /*
1928      * stops computation if
1929      * - 24 (degBound)
1930      *   && upper degree is bigger than Kstd1_deg
1931      */
1932      while ((strat->Ll >= 0)
1933        && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
1934        && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg)
1935      )
1936      {
1937        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1938        //if (TEST_OPT_PROT)
1939        //{
1940        //   PrintS("D"); mflush();
1941        //}
1942      }
1943      if (strat->Ll<0) break;
1944      else strat->noClearS=TRUE;
1945    }
1946    strat->P = strat->L[strat->Ll];/*- picks the last element from the lazyset L -*/
1947    if (strat->Ll==0) strat->interpt=TRUE;
1948    strat->Ll--;
1949
1950    // create the real Spoly
1951    assume(pNext(strat->P.p) != strat->tail);
1952
1953    if (strat->P.p1 == NULL)
1954    {
1955      // for input polys, prepare reduction (buckets !)
1956      strat->P.SetLength(strat->length_pLength);
1957      strat->P.PrepareRed(strat->use_buckets);
1958    }
1959
1960    if (!strat->P.IsNull())
1961    {
1962      // might be NULL from noether !!!
1963      if (TEST_OPT_PROT)
1964        message(strat->P.ecart+strat->P.GetpFDeg(),&olddeg,&reduc,strat, red_result);
1965      // reduce
1966      red_result = strat->red(&strat->P,strat);
1967    }
1968
1969    if (! strat->P.IsNull())
1970    {
1971      strat->P.GetP();
1972      // statistics
1973      if (TEST_OPT_PROT) PrintS("s");
1974      // normalization
1975      if (!TEST_OPT_INTSTRATEGY)
1976        strat->P.pNorm();
1977      // tailreduction
1978      strat->P.p = redtail(&(strat->P),strat->sl,strat);
1979      // set ecart -- might have changed because of tail reductions
1980      if ((!strat->noTailReduction) && (!strat->honey))
1981        strat->initEcart(&strat->P);
1982      // cancel unit
1983      cancelunit(&strat->P);
1984      // for char 0, clear denominators
1985      if (TEST_OPT_INTSTRATEGY)
1986        strat->P.pCleardenom();
1987
1988      // put in T
1989      enterT(strat->P,strat);
1990      // build new pairs
1991      enterpairs(strat->P.p,strat->sl,strat->P.ecart,0,strat, strat->tl);
1992      // put in S
1993      strat->enterS(strat->P,
1994                    posInS(strat,strat->sl,strat->P.p, strat->P.ecart),
1995                    strat, strat->tl);
1996
1997
1998      // clear strat->P
1999      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
2000      strat->P.lcm=NULL;
2001
2002      if (strat->sl>srmax) srmax = strat->sl; /*stat.*/
2003      if (strat->Ll>lrmax) lrmax = strat->Ll;
2004
2005
2006
2007      // //////////////////////////////////////////////////////////
2008      // SCA:
2009      const poly pSave = strat->P.p;
2010      const poly pNext = pNext(pSave);
2011
2012      if(pNext != NULL)
2013      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
2014      if( p_GetExp(pSave, i, currRing) != 0 )
2015      {
2016
2017        assume(p_GetExp(pSave, i, currRing) == 1);
2018
2019        const poly pNew = sca_pp_Mult_xi_pp(i, pNext, currRing);
2020
2021#ifdef PDEBUG
2022        p_Test(pNew, currRing);
2023#endif
2024
2025        if( pNew == NULL) continue;
2026
2027        LObject h(pNew); // h = x_i * strat->P
2028
2029        if (TEST_OPT_INTSTRATEGY)
2030           h.pCleardenom(); // also does a pContent
2031        else
2032          h.pNorm();
2033
2034        strat->initEcart(&h);
2035        h.sev = pGetShortExpVector(h.p);
2036
2037        int pos;
2038
2039        if (strat->Ll==-1)
2040          pos = 0;
2041        else
2042          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
2043
2044        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
2045
2046        if (strat->Ll>lrmax) lrmax = strat->Ll;
2047      }
2048
2049#ifdef KDEBUG
2050      // make sure kTest_TS does not complain about strat->P
2051      memset(&strat->P,0,sizeof(strat->P));
2052#endif
2053    }
2054#if 0
2055    if (strat->kHEdgeFound)
2056    {
2057      if ((BTEST1(27))
2058      || ((TEST_OPT_MULTBOUND) && (scMult0Int((strat->Shdl)) < mu)))
2059      {
2060        // obachman: is this still used ???
2061        /*
2062        * stops computation if strat->kHEdgeFound and
2063        * - 27 (finiteDeterminacyTest)
2064        * or
2065        * - 23
2066        *   (multBound)
2067        *   && multiplicity of the ideal is smaller then a predefined number mu
2068        */
2069        while (strat->Ll >= 0) deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
2070      }
2071    }
2072#endif
2073    kTest_TS(strat);
2074  }
2075  /*- complete reduction of the standard basis------------------------ -*/
2076  if (TEST_OPT_REDSB) completeReduce(strat);
2077  /*- release temp data------------------------------- -*/
2078  exitBuchMora(strat);
2079  /*- polynomials used for HECKE: HC, noether -*/
2080  if (BTEST1(27))
2081  {
2082    if (strat->kHEdge!=NULL)
2083      Kstd1_mu=pFDeg(strat->kHEdge,currRing);
2084    else
2085      Kstd1_mu=-1;
2086  }
2087  pDelete(&strat->kHEdge);
2088  strat->update = TRUE; //???
2089  strat->lastAxis = 0; //???
2090  pDelete(&strat->kNoether);
2091  omFreeSize((ADDRESS)strat->NotUsedAxis,(pVariables+1)*sizeof(BOOLEAN));
2092  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
2093  if (TEST_OPT_WEIGHTM)
2094  {
2095    pRestoreDegProcs(pFDegOld, pLDegOld);
2096    if (ecartWeights)
2097    {
2098      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
2099      ecartWeights=NULL;
2100    }
2101  }
2102  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
2103  idTest(strat->Shdl);
2104 
2105  id_Delete( &tempF, currRing);
2106 
2107  return (strat->Shdl);
2108}
2109
2110
2111
2112
2113
2114
2115void sca_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
2116{
2117
2118  // "commutative" procedures:
2119  rGR->p_Procs->p_Mult_mm     = sca_p_Mult_mm;
2120  rGR->p_Procs->pp_Mult_mm    = sca_pp_Mult_mm;
2121
2122  p_Procs->p_Mult_mm          = sca_p_Mult_mm;
2123  p_Procs->pp_Mult_mm         = sca_pp_Mult_mm;
2124
2125  // non-commutaitve
2126  rGR->nc->p_Procs.mm_Mult_p   = sca_mm_Mult_p;
2127  rGR->nc->p_Procs.mm_Mult_pp  = sca_mm_Mult_pp;
2128
2129
2130  if (pOrdSgn==-1)
2131  {
2132#ifdef PDEBUG
2133//           Print("Local case => GB == mora!\n");
2134#endif
2135    rGR->nc->p_Procs.GB          = sca_mora; // local ordering => Mora, otherwise - Buchberger!
2136  }
2137  else
2138  {
2139#ifdef PDEBUG
2140//           Print("Global case => GB == bba!\n");
2141#endif
2142    rGR->nc->p_Procs.GB          = sca_gr_bba; // sca_bba?
2143  }
2144
2145
2146//   rGR->nc->p_Procs.GlobalGB    = sca_gr_bba;
2147//   rGR->nc->p_Procs.LocalGB     = sca_mora;
2148
2149
2150//   rGR->nc->p_Procs.SPoly         = sca_SPoly;
2151//   rGR->nc->p_Procs.ReduceSPoly   = sca_ReduceSpoly;
2152
2153#if 0
2154
2155        // Multiplication procedures:
2156
2157        p_Procs->p_Mult_mm   = sca_p_Mult_mm;
2158        _p_procs->p_Mult_mm  = sca_p_Mult_mm;
2159
2160        p_Procs->pp_Mult_mm  = sca_pp_Mult_mm;
2161        _p_procs->pp_Mult_mm = sca_pp_Mult_mm;
2162
2163        r->nc->mmMultP()     = sca_mm_Mult_p;
2164        r->nc->mmMultPP()    = sca_mm_Mult_pp;
2165
2166        r->nc->GB()            = sca_gr_bba;
2167/*
2168        // ??????????????????????????????????????? coefficients swell...
2169        r->nc->SPoly()         = sca_SPoly;
2170        r->nc->ReduceSPoly()   = sca_ReduceSpoly;
2171*/
2172//         r->nc->BucketPolyRed() = gnc_kBucketPolyRed;
2173//         r->nc->BucketPolyRed_Z()= gnc_kBucketPolyRed_Z;
2174
2175#endif
2176}
2177
2178
2179// bi-Degree (x, y) of monomial "m"
2180// x-es and y-s are weighted by wx and wy resp.
2181// [optional] components have weights by wCx and wCy.
2182inline void m_GetBiDegree(const poly m, 
2183  const intvec *wx, const intvec *wy, 
2184  const intvec *wCx, const intvec *wCy, 
2185  int& dx, int& dy, const ring r)
2186{
2187  const unsigned int N  = r->N;
2188 
2189  p_Test(m, r);
2190 
2191  assume( wx != NULL );
2192  assume( wy != NULL );
2193 
2194  assume( wx->cols() == 1 );
2195  assume( wy->cols() == 1 );
2196
2197  assume( wx->rows() >= N );
2198  assume( wy->rows() >= N );
2199
2200  int x = 0;
2201  int y = 0;
2202 
2203  for(int i = N; i > 0; i--)
2204  {
2205    const int d = p_GetExp(m, i, r);
2206    x += d * (*wx)[i-1];
2207    y += d * (*wy)[i-1];
2208  }
2209 
2210  if( (wCx != NULL) && (wCy != NULL) )
2211  {
2212    const int c = p_GetComp(m, r);
2213
2214    if( wCx->range(c) )
2215      x += (*wCx)[c];
2216
2217    if( wCy->range(c) )
2218      x += (*wCy)[c];     
2219  }
2220 
2221  dx = x;
2222  dy = y;
2223}
2224
2225// returns true if polynom p is bi-homogenous with respect to the given weights
2226// simultaneously sets bi-Degree
2227bool p_IsBiHomogeneous(const poly p, 
2228  const intvec *wx, const intvec *wy, 
2229  const intvec *wCx, const intvec *wCy, 
2230  int &dx, int &dy,
2231  const ring r)
2232{
2233  if( p == NULL ) 
2234  {
2235    dx = 0;
2236    dy = 0;
2237    return true;
2238  }
2239
2240  poly q = p;
2241
2242
2243  int ddx, ddy;
2244
2245  m_GetBiDegree( q, wx, wy, wCx, wCy, ddx, ddy, r); // get bi degree of lm(p)
2246
2247  pIter(q);
2248
2249  for(; q != NULL; pIter(q) )
2250  {
2251    int x, y;   
2252   
2253    m_GetBiDegree( q, wx, wy, wCx, wCy, x, y, r); // get bi degree of q
2254   
2255    if ( (x != ddx) || (y != ddy) ) return false;
2256  }
2257 
2258  dx = ddx;
2259  dy = ddy;
2260
2261  return true;
2262}
2263
2264
2265// returns true if id is bi-homogenous without respect to the given weights
2266bool id_IsBiHomogeneous(const ideal id, 
2267  const intvec *wx, const intvec *wy, 
2268  const intvec *wCx, const intvec *wCy, 
2269  const ring r)
2270{
2271  if (id == NULL) return true; // zero ideal
2272
2273  const int iSize = IDELEMS(id);
2274
2275  if (iSize == 0) return true;
2276
2277  bool b = true;
2278  int x, y;
2279
2280  for(int i = iSize - 1; (i >= 0 ) && b; i--)
2281    b = p_IsBiHomogeneous(id->m[i], wx, wy, wCx, wCy, x, y, r);
2282
2283  return b;
2284}
2285
2286
2287// returns an intvector with [nvars(r)] integers [1/0]
2288// 1 - for commutative variables
2289// 0 - for anticommutative variables
2290intvec *ivGetSCAXVarWeights(const ring r)
2291{
2292  const unsigned int N  = r->N;
2293
2294  const int CommutativeVariable = 1;
2295  const int AntiCommutativeVariable = 0;
2296 
2297  intvec* w = new intvec(N, 1, CommutativeVariable);
2298 
2299  if( rIsSCA(r) )
2300  {
2301    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2302    const unsigned int m_iLastAltVar  = scaLastAltVar(r); 
2303
2304    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2305    {
2306      (*w)[i-1] = AntiCommutativeVariable;
2307    } 
2308  }
2309  return w;
2310}
2311
2312
2313// returns an intvector with [nvars(r)] integers [1/0]
2314// 0 - for commutative variables
2315// 1 - for anticommutative variables
2316intvec *ivGetSCAYVarWeights(const ring r)
2317{
2318  const unsigned int N  = r->N;
2319
2320  const int CommutativeVariable = 0;
2321  const int AntiCommutativeVariable = 1;
2322 
2323  intvec* w = new intvec(N, 1, CommutativeVariable);
2324 
2325  if( rIsSCA(r) )
2326  {
2327    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2328    const unsigned int m_iLastAltVar  = scaLastAltVar(r); 
2329
2330    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2331    {
2332      (*w)[i-1] = AntiCommutativeVariable;
2333    } 
2334  }
2335  return w;
2336}
2337
2338
2339
2340
2341// reduce term lt(m) modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar:
2342// either create a copy of m or return NULL
2343inline poly m_KillSquares(const poly m, 
2344  const unsigned int iFirstAltVar, const unsigned int iLastAltVar, 
2345  const ring r)
2346{ 
2347#ifdef PDEBUG
2348  p_Test(m, r);
2349
2350#if 0
2351  Print("m_KillSquares, m = "); // !
2352  p_Write(m, r);
2353#endif
2354#endif
2355
2356  assume( m != NULL );
2357
2358  for(int k = iFirstAltVar; k <= iLastAltVar; k++)
2359    if( p_GetExp(m, k, r) > 1 )
2360      return NULL; 
2361 
2362  return p_Head(m, r);
2363}
2364
2365
2366// reduce polynomial p modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2367poly p_KillSquares(const poly p, 
2368  const unsigned int iFirstAltVar, const unsigned int iLastAltVar, 
2369  const ring r)
2370{ 
2371#ifdef PDEBUG
2372  p_Test(p, r);
2373
2374#if 0
2375  Print("p_KillSquares, p = "); // !
2376  p_Write(p, r);
2377#endif
2378#endif
2379
2380
2381  if( p == NULL )
2382    return NULL;
2383
2384  poly pResult = NULL;
2385  poly* ppPrev = &pResult;
2386
2387  for( poly q = p; q!= NULL; pIter(q) )
2388  {
2389#ifdef PDEBUG
2390    p_Test(q, r);
2391#endif
2392
2393    // terms will be in the same order because of quasi-ordering!   
2394    poly v = m_KillSquares(q, iFirstAltVar, iLastAltVar, r);
2395
2396    if( v != NULL )
2397    {
2398      *ppPrev = v;
2399      ppPrev = &pNext(v);
2400    }
2401
2402  }
2403
2404#ifdef PDEBUG
2405  p_Test(pResult, r);
2406#if 0
2407  Print("p_KillSquares => "); // !
2408  p_Write(pResult, r);
2409#endif
2410#endif
2411
2412  return(pResult);
2413}
2414 
2415
2416
2417
2418// reduces ideal id modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2419// returns the reduced ideal or zero ideal.
2420ideal id_KillSquares(const ideal id, 
2421  const unsigned int iFirstAltVar, const unsigned int iLastAltVar, 
2422  const ring r)
2423{
2424  if (id == NULL) return id; // zero ideal
2425
2426  const int iSize = id->idelems();
2427
2428  if (iSize == 0) return id;
2429 
2430  ideal temp = idInit(iSize, 1);
2431 
2432#if 0
2433   PrintS("<id_KillSquares>\n");
2434  {
2435    Print("ideal id: \n");
2436    for (int i = 0; i < id->idelems(); i++)
2437    {
2438      Print("; id[%d] = ", i+1);
2439      p_Write(id->m[i], r);
2440    }
2441    Print(";\n");
2442    PrintLn();
2443  }
2444#endif
2445 
2446
2447  for (int j = 0; j < iSize; j++)
2448    temp->m[j] = p_KillSquares(id->m[j], iFirstAltVar, iLastAltVar, r);
2449
2450  idSkipZeroes(temp);
2451
2452#if 0
2453   PrintS("<id_KillSquares>\n");
2454  {
2455    Print("ideal temp: \n");
2456    for (int i = 0; i < temp->idelems(); i++)
2457    {
2458      Print("; temp[%d] = ", i+1);
2459      p_Write(temp->m[i], r);
2460    }
2461    Print(";\n");
2462    PrintLn();
2463  }
2464   PrintS("</id_KillSquares>\n");
2465#endif
2466
2467  return temp;
2468}
2469
2470
2471
2472
2473#endif
Note: See TracBrowser for help on using the repository browser.