source: git/kernel/sca.cc @ 9108d3

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