source: git/kernel/sca.cc @ 096c99

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