source: git/kernel/sca.cc @ 0a8ee5

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