source: git/kernel/sca.cc @ 651f6f

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