source: git/kernel/sca.cc @ b1c0a9

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