source: git/kernel/sca.cc @ 71b676

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