source: git/kernel/sca.cc @ f115489

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