source: git/kernel/sca.cc @ a610ee

spielwiese
Last change on this file since a610ee was a610ee, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: Print vs. PrintS git-svn-id: file:///usr/local/Singular/svn/trunk@11452 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 61.6 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.32 2009-02-23 13:50:52 Singular 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_Test(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_Test(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_Test(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_Test(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 = 1;
1480  int e = N; 
1481
1482  if(rIsSCA(rG))
1483  {
1484    b = si_max(b, scaFirstAltVar(rG));
1485    e = si_min(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 - 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  /* compute------------------------------------------------------- */
1753  while (strat->Ll >= 0)
1754  {
1755    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
1756
1757#ifdef KDEBUG
1758//     loop_count++;
1759    if (TEST_OPT_DEBUG) messageSets(strat);
1760#endif
1761
1762    if (strat->Ll== 0) strat->interpt=TRUE;
1763
1764    if (TEST_OPT_DEGBOUND
1765        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1766            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
1767    {
1768      /*
1769       *stops computation if
1770       * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
1771       *a predefined number Kstd1_deg
1772       */
1773      while ((strat->Ll >= 0)
1774  && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
1775        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1776            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))
1777  )
1778        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1779      if (strat->Ll<0) break;
1780      else strat->noClearS=TRUE;
1781    }
1782
1783    /* picks the last element from the lazyset L */
1784    strat->P = strat->L[strat->Ll];
1785    strat->Ll--;
1786
1787
1788//    assume(pNext(strat->P.p) != strat->tail);
1789
1790    if(strat->P.IsNull()) continue;
1791
1792    if (pNext(strat->P.p) == strat->tail)
1793    {
1794      // deletes the short spoly
1795      pLmFree(strat->P.p);
1796
1797      strat->P.p = nc_CreateSpoly(strat->P.p1, strat->P.p2, currRing);
1798      if (strat->P.p!=NULL) strat->initEcart(&strat->P);
1799    }//    else
1800
1801
1802    if(strat->P.IsNull()) continue;
1803
1804    if (strat->P.p1 == NULL)
1805    {
1806//       if (strat->minim > 0)
1807//         strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
1808
1809
1810      // for input polys, prepare reduction
1811        strat->P.PrepareRed(strat->use_buckets);
1812    }
1813
1814    if (TEST_OPT_PROT)
1815      message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
1816              &olddeg,&reduc,strat, red_result);
1817
1818    /* reduction of the element choosen from L */
1819    red_result = strat->red(&strat->P,strat);
1820
1821
1822    // reduction to non-zero new poly
1823    if (red_result == 1)
1824    {
1825      /* statistic */
1826      if (TEST_OPT_PROT) PrintS("s");
1827
1828      // get the polynomial (canonicalize bucket, make sure P.p is set)
1829      strat->P.GetP(strat->lmBin);
1830
1831      int pos = posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
1832
1833      // reduce the tail and normalize poly
1834      if (TEST_OPT_INTSTRATEGY)
1835      {
1836        strat->P.pCleardenom();
1837        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1838        {
1839          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT); // !!!
1840          strat->P.pCleardenom();
1841        }
1842      }
1843      else
1844      {
1845        strat->P.pNorm();
1846        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1847          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
1848      }
1849      strat->P.is_normalized=nIsOne(pGetCoeff(strat->P.p));
1850
1851#ifdef KDEBUG
1852      if (TEST_OPT_DEBUG){PrintS(" ns:");p_wrp(strat->P.p,currRing);PrintLn();}
1853#endif
1854
1855//       // min_std stuff
1856//       if ((strat->P.p1==NULL) && (strat->minim>0))
1857//       {
1858//         if (strat->minim==1)
1859//         {
1860//           strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
1861//           p_Delete(&strat->P.p2, currRing, strat->tailRing);
1862//         }
1863//         else
1864//         {
1865//           strat->M->m[minimcnt]=strat->P.p2;
1866//           strat->P.p2=NULL;
1867//         }
1868//         if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
1869//           pNext(strat->M->m[minimcnt])
1870//             = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
1871//                                            strat->tailRing, currRing,
1872//                                            currRing->PolyBin);
1873//         minimcnt++;
1874//       }
1875
1876      // enter into S, L, and T
1877      //if(withT)
1878        enterT(strat->P, strat);
1879
1880      // L
1881      enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
1882
1883      // posInS only depends on the leading term
1884      strat->enterS(strat->P, pos, strat, strat->tl);
1885
1886//       if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
1887
1888//      Print("[%d]",hilbeledeg);
1889      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
1890
1891      if (strat->sl>srmax) srmax = strat->sl;
1892
1893      // //////////////////////////////////////////////////////////
1894      // SCA:
1895      const poly pSave = strat->P.p;
1896      const poly p_next = pNext(pSave);
1897
1898//       if(0)
1899      if( p_next != NULL )
1900      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
1901      if( p_GetExp(pSave, i, currRing) != 0 )
1902      {
1903        assume(p_GetExp(pSave, i, currRing) == 1);
1904        const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
1905
1906#ifdef PDEBUG
1907        p_Test(p_new, currRing);
1908#endif
1909
1910        if( p_new == NULL) continue;
1911
1912        LObject h(p_new); // h = x_i * strat->P
1913
1914        if (TEST_OPT_INTSTRATEGY)
1915        {
1916          pContent(h.p);
1917        }
1918        else
1919        {
1920          h.pNorm();
1921        }
1922
1923        strat->initEcart(&h);
1924        h.sev = pGetShortExpVector(h.p);
1925
1926        int pos;
1927        if (strat->Ll==-1)
1928          pos =0;
1929        else
1930          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1931
1932        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
1933/*
1934        h.sev = pGetShortExpVector(h.p);
1935        strat->initEcart(&h);
1936
1937        h.PrepareRed(strat->use_buckets);
1938
1939        // reduction of the element choosen from L(?)
1940        red_result = strat->red(&h,strat);
1941
1942        // reduction to non-zero new poly
1943        if (red_result != 1) continue;
1944
1945
1946        int pos = posInS(strat,strat->sl,h.p,h.ecart);
1947
1948        // reduce the tail and normalize poly
1949        if (TEST_OPT_INTSTRATEGY)
1950        {
1951          h.pCleardenom();
1952          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1953          {
1954            h.p = redtailBba(&(h),pos-1,strat, withT); // !!!
1955            h.pCleardenom();
1956          }
1957        }
1958        else
1959        {
1960          h.pNorm();
1961          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1962            h.p = redtailBba(&(h),pos-1,strat, withT);
1963        }
1964
1965#ifdef KDEBUG
1966        if (TEST_OPT_DEBUG){PrintS(" N:");h.wrp();PrintLn();}
1967#endif
1968
1969//        h.PrepareRed(strat->use_buckets); // ???
1970
1971        h.sev = pGetShortExpVector(h.p);
1972        strat->initEcart(&h);
1973
1974        if (strat->Ll==-1)
1975          pos = 0;
1976        else
1977          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1978
1979         enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);*/
1980
1981      } // for all x_i \in Ann(lm(P))
1982    } // if red(P) != NULL
1983
1984//     else if (strat->P.p1 == NULL && strat->minim > 0)
1985//     {
1986//       p_Delete(&strat->P.p2, currRing, strat->tailRing);
1987//     }
1988
1989
1990   
1991#ifdef KDEBUG
1992//    memset(&(strat->P), 0, sizeof(strat->P));
1993#endif
1994
1995    kTest_TS(strat); // even of T is not used!
1996
1997//     Print("\n$\n");
1998
1999  }
2000
2001#ifdef KDEBUG
2002  if (TEST_OPT_DEBUG) messageSets(strat);
2003#endif
2004
2005  /* complete reduction of the standard basis--------- */
2006
2007  if (TEST_OPT_REDSB){
2008    completeReduce(strat);
2009  }
2010
2011 
2012  /* release temp data-------------------------------- */
2013
2014  exitBuchMora(strat); // cleanT!
2015
2016  id_Delete(&tempF, currRing);
2017 
2018  if (TEST_OPT_WEIGHTM)
2019  {
2020    pRestoreDegProcs(pFDegOld, pLDegOld);
2021    if (ecartWeights)
2022    {
2023      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
2024      ecartWeights=NULL;
2025    }
2026  }
2027
2028  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
2029
2030
2031 
2032  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
2033
2034
2035  if (TEST_OPT_REDSB) // ???
2036  {
2037    // must be at the very end (after exitBuchMora) as it changes the S set!!!
2038    ideal I = strat->Shdl;
2039    ideal erg = kInterRedOld(I,tempQ);
2040    assume(I!=erg);
2041    id_Delete(&I, currRing);
2042    strat->Shdl = erg;
2043  }
2044 
2045 
2046#if MYTEST
2047  PrintS("\n\n</sca_bba>\n\n");
2048#endif
2049 
2050  return (strat->Shdl);
2051}
2052
2053
2054// //////////////////////////////////////////////////////////////////////////////
2055// sca mora...
2056
2057// returns TRUE if mora should use buckets, false otherwise
2058static BOOLEAN kMoraUseBucket(kStrategy strat)
2059{
2060#ifdef MORA_USE_BUCKETS
2061  if (TEST_OPT_NOT_BUCKETS)
2062    return FALSE;
2063  if (strat->red == redFirst)
2064  {
2065#ifdef NO_LDEG
2066    if (!strat->syzComp)
2067      return TRUE;
2068#else
2069    if ((strat->homog || strat->honey) && !strat->syzComp)
2070      return TRUE;
2071#endif
2072  }
2073  else
2074  {
2075    assume(strat->red == redEcart);
2076    if (strat->honey && !strat->syzComp)
2077      return TRUE;
2078  }
2079#endif
2080  return FALSE;
2081}
2082
2083
2084#ifdef HAVE_ASSUME
2085static int sca_mora_count = 0;
2086static int sca_mora_loop_count;
2087#endif
2088
2089// ideal sca_mora (ideal F, ideal Q, intvec *w, intvec *, kStrategy strat)
2090ideal sca_mora(const ideal F, const ideal Q, const intvec *w, const intvec *, kStrategy strat)
2091{
2092  assume(rIsSCA(currRing));
2093
2094  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
2095  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
2096
2097  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
2098
2099  ideal tempQ = Q;
2100
2101  if(Q == currQuotient)
2102    tempQ = SCAQuotient(currRing);
2103
2104  bool bIdHomog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
2105
2106  assume( !bIdHomog || strat->homog ); //  bIdHomog =====[implies]>>>>> strat->homog
2107
2108  strat->homog = strat->homog && bIdHomog;
2109
2110#ifdef PDEBUG
2111  assume( strat->homog == bIdHomog );
2112#endif /*PDEBUG*/
2113
2114#ifdef HAVE_ASSUME
2115  sca_mora_count++;
2116  sca_mora_loop_count = 0;
2117#endif
2118
2119#ifdef KDEBUG
2120  om_Opts.MinTrack = 5;
2121#endif
2122
2123  int srmax;
2124  int lrmax = 0;
2125  int olddeg = 0;
2126  int reduc = 0;
2127  int red_result = 1;
2128//  int hilbeledeg=1;
2129  int hilbcount=0;
2130
2131  strat->update = TRUE;
2132  /*- setting global variables ------------------- -*/
2133  initBuchMoraCrit(strat);
2134//   initHilbCrit(F,NULL,&hilb,strat); // no Q!
2135  initMora(tempF, strat);
2136  initBuchMoraPos(strat);
2137  /*Shdl=*/initBuchMora(tempF, tempQ, strat); // temp Q, F!
2138//   if (TEST_OPT_FASTHC) missingAxis(&strat->lastAxis,strat);
2139  /*updateS in initBuchMora has Hecketest
2140  * and could have put strat->kHEdgdeFound FALSE*/
2141#if 0
2142  if (ppNoether!=NULL)
2143  {
2144    strat->kHEdgeFound = TRUE;
2145  }
2146  if (strat->kHEdgeFound && strat->update)
2147  {
2148    firstUpdate(strat);
2149    updateLHC(strat);
2150    reorderL(strat);
2151  }
2152  if (TEST_OPT_FASTHC && (strat->lastAxis) && strat->posInLOldFlag)
2153  {
2154    strat->posInLOld = strat->posInL;
2155    strat->posInLOldFlag = FALSE;
2156    strat->posInL = posInL10;
2157    updateL(strat);
2158    reorderL(strat);
2159  }
2160#endif
2161
2162  srmax = strat->sl;
2163  kTest_TS(strat);
2164
2165  strat->use_buckets = kMoraUseBucket(strat);
2166  /*- compute-------------------------------------------*/
2167
2168#undef HAVE_TAIL_RING
2169
2170#ifdef HAVE_TAIL_RING
2171//  if (strat->homog && strat->red == redFirst)
2172//     kStratInitChangeTailRing(strat);
2173#endif
2174
2175
2176  while (strat->Ll >= 0)
2177  {
2178#ifdef HAVE_ASSUME
2179    sca_mora_loop_count++;
2180#endif
2181    if (lrmax< strat->Ll) lrmax=strat->Ll; /*stat*/
2182    //test_int_std(strat->kIdeal);
2183#ifdef KDEBUG
2184    if (TEST_OPT_DEBUG) messageSets(strat);
2185#endif
2186    if (TEST_OPT_DEGBOUND
2187    && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg))
2188    {
2189      /*
2190      * stops computation if
2191      * - 24 (degBound)
2192      *   && upper degree is bigger than Kstd1_deg
2193      */
2194      while ((strat->Ll >= 0)
2195        && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
2196        && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg)
2197      )
2198      {
2199        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
2200        //if (TEST_OPT_PROT)
2201        //{
2202        //   PrintS("D"); mflush();
2203        //}
2204      }
2205      if (strat->Ll<0) break;
2206      else strat->noClearS=TRUE;
2207    }
2208    strat->P = strat->L[strat->Ll];/*- picks the last element from the lazyset L -*/
2209    if (strat->Ll==0) strat->interpt=TRUE;
2210    strat->Ll--;
2211
2212    // create the real Spoly
2213//    assume(pNext(strat->P.p) != strat->tail);
2214
2215    if(strat->P.IsNull()) continue;
2216
2217
2218    if( pNext(strat->P.p) == strat->tail )
2219    {
2220      // deletes the int spoly and computes SPoly
2221      pLmFree(strat->P.p); // ???
2222      strat->P.p = nc_CreateSpoly(strat->P.p1, strat->P.p2, currRing);
2223    }
2224
2225
2226
2227    if (strat->P.p1 == NULL)
2228    {
2229      // for input polys, prepare reduction (buckets !)
2230      strat->P.SetLength(strat->length_pLength);
2231      strat->P.PrepareRed(strat->use_buckets);
2232    }
2233
2234    if (!strat->P.IsNull())
2235    {
2236      // might be NULL from noether !!!
2237      if (TEST_OPT_PROT)
2238        message(strat->P.ecart+strat->P.GetpFDeg(),&olddeg,&reduc,strat, red_result);
2239      // reduce
2240      red_result = strat->red(&strat->P,strat);
2241    }
2242
2243    if (! strat->P.IsNull())
2244    {
2245      strat->P.GetP();
2246      // statistics
2247      if (TEST_OPT_PROT) PrintS("s");
2248      // normalization
2249      if (!TEST_OPT_INTSTRATEGY)
2250        strat->P.pNorm();
2251      // tailreduction
2252      strat->P.p = redtail(&(strat->P),strat->sl,strat);
2253      // set ecart -- might have changed because of tail reductions
2254      if ((!strat->noTailReduction) && (!strat->honey))
2255        strat->initEcart(&strat->P);
2256      // cancel unit
2257      cancelunit(&strat->P);
2258      // for char 0, clear denominators
2259      if (TEST_OPT_INTSTRATEGY)
2260        strat->P.pCleardenom();
2261
2262      // put in T
2263      enterT(strat->P,strat);
2264      // build new pairs
2265      enterpairs(strat->P.p,strat->sl,strat->P.ecart,0,strat, strat->tl);
2266      // put in S
2267      strat->enterS(strat->P,
2268                    posInS(strat,strat->sl,strat->P.p, strat->P.ecart),
2269                    strat, strat->tl);
2270
2271
2272      // clear strat->P
2273      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
2274      strat->P.lcm=NULL;
2275
2276      if (strat->sl>srmax) srmax = strat->sl; /*stat.*/
2277      if (strat->Ll>lrmax) lrmax = strat->Ll;
2278
2279
2280
2281      // //////////////////////////////////////////////////////////
2282      // SCA:
2283      const poly pSave = strat->P.p;
2284      const poly p_next = pNext(pSave);
2285
2286      if(p_next != NULL)
2287      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
2288      if( p_GetExp(pSave, i, currRing) != 0 )
2289      {
2290
2291        assume(p_GetExp(pSave, i, currRing) == 1);
2292
2293        const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
2294
2295#ifdef PDEBUG
2296        p_Test(p_new, currRing);
2297#endif
2298
2299        if( p_new == NULL) continue;
2300
2301        LObject h(p_new); // h = x_i * strat->P
2302
2303        if (TEST_OPT_INTSTRATEGY)
2304           h.pCleardenom(); // also does a pContent
2305        else
2306          h.pNorm();
2307
2308        strat->initEcart(&h);
2309        h.sev = pGetShortExpVector(h.p);
2310
2311        int pos;
2312
2313        if (strat->Ll==-1)
2314          pos = 0;
2315        else
2316          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
2317
2318        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
2319
2320        if (strat->Ll>lrmax) lrmax = strat->Ll;
2321      }
2322
2323#ifdef KDEBUG
2324      // make sure kTest_TS does not complain about strat->P
2325      memset(&strat->P,0,sizeof(strat->P));
2326#endif
2327    }
2328#if 0
2329    if (strat->kHEdgeFound)
2330    {
2331      if ((BTEST1(27))
2332      || ((TEST_OPT_MULTBOUND) && (scMult0Int((strat->Shdl)) < mu)))
2333      {
2334        // obachman: is this still used ???
2335        /*
2336        * stops computation if strat->kHEdgeFound and
2337        * - 27 (finiteDeterminacyTest)
2338        * or
2339        * - 23
2340        *   (multBound)
2341        *   && multiplicity of the ideal is smaller then a predefined number mu
2342        */
2343        while (strat->Ll >= 0) deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
2344      }
2345    }
2346#endif
2347    kTest_TS(strat);
2348  }
2349  /*- complete reduction of the standard basis------------------------ -*/
2350  if (TEST_OPT_REDSB) completeReduce(strat);
2351  /*- release temp data------------------------------- -*/
2352  exitBuchMora(strat);
2353  /*- polynomials used for HECKE: HC, noether -*/
2354  if (BTEST1(27))
2355  {
2356    if (strat->kHEdge!=NULL)
2357      Kstd1_mu=pFDeg(strat->kHEdge,currRing);
2358    else
2359      Kstd1_mu=-1;
2360  }
2361  pDelete(&strat->kHEdge);
2362  strat->update = TRUE; //???
2363  strat->lastAxis = 0; //???
2364  pDelete(&strat->kNoether);
2365  omFreeSize((ADDRESS)strat->NotUsedAxis,(pVariables+1)*sizeof(BOOLEAN));
2366  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
2367  if (TEST_OPT_WEIGHTM)
2368  {
2369    pRestoreDegProcs(pFDegOld, pLDegOld);
2370    if (ecartWeights)
2371    {
2372      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
2373      ecartWeights=NULL;
2374    }
2375  }
2376  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
2377  idTest(strat->Shdl);
2378
2379  id_Delete( &tempF, currRing);
2380
2381  return (strat->Shdl);
2382}
2383
2384
2385
2386
2387
2388
2389void sca_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
2390{
2391
2392  // "commutative" procedures:
2393  rGR->p_Procs->p_Mult_mm     = sca_p_Mult_mm;
2394  rGR->p_Procs->pp_Mult_mm    = sca_pp_Mult_mm;
2395
2396  p_Procs->p_Mult_mm          = sca_p_Mult_mm;
2397  p_Procs->pp_Mult_mm         = sca_pp_Mult_mm;
2398
2399  // non-commutaitve
2400  rGR->GetNC()->p_Procs.mm_Mult_p   = sca_mm_Mult_p;
2401  rGR->GetNC()->p_Procs.mm_Mult_pp  = sca_mm_Mult_pp;
2402
2403
2404  if (rHasLocalOrMixedOrdering(rGR))
2405  {
2406#ifdef PDEBUG
2407//           Print("Local case => GB == mora!\n");
2408#endif
2409    rGR->GetNC()->p_Procs.GB          = sca_mora; // local ordering => Mora, otherwise - Buchberger!
2410  }
2411  else
2412  {
2413#ifdef PDEBUG
2414//           Print("Global case => GB == bba!\n");
2415#endif
2416    rGR->GetNC()->p_Procs.GB          = sca_bba; // sca_gr_bba; // sca_bba? // sca_bba;
2417  }
2418
2419
2420//   rGR->GetNC()->p_Procs.GlobalGB    = sca_gr_bba;
2421//   rGR->GetNC()->p_Procs.LocalGB     = sca_mora;
2422
2423
2424//   rGR->GetNC()->p_Procs.SPoly         = sca_SPoly;
2425//   rGR->GetNC()->p_Procs.ReduceSPoly   = sca_ReduceSpoly;
2426
2427#if 0
2428
2429        // Multiplication procedures:
2430
2431        p_Procs->p_Mult_mm   = sca_p_Mult_mm;
2432        _p_procs->p_Mult_mm  = sca_p_Mult_mm;
2433
2434        p_Procs->pp_Mult_mm  = sca_pp_Mult_mm;
2435        _p_procs->pp_Mult_mm = sca_pp_Mult_mm;
2436
2437        r->GetNC()->mmMultP()     = sca_mm_Mult_p;
2438        r->GetNC()->mmMultPP()    = sca_mm_Mult_pp;
2439
2440        r->GetNC()->GB()            = sca_gr_bba;
2441/*
2442        // ??????????????????????????????????????? coefficients swell...
2443        r->GetNC()->SPoly()         = sca_SPoly;
2444        r->GetNC()->ReduceSPoly()   = sca_ReduceSpoly;
2445*/
2446//         r->GetNC()->BucketPolyRed() = gnc_kBucketPolyRed;
2447//         r->GetNC()->BucketPolyRed_Z()= gnc_kBucketPolyRed_Z;
2448
2449#endif
2450}
2451
2452
2453// bi-Degree (x, y) of monomial "m"
2454// x-es and y-s are weighted by wx and wy resp.
2455// [optional] components have weights by wCx and wCy.
2456inline void m_GetBiDegree(const poly m,
2457  const intvec *wx, const intvec *wy,
2458  const intvec *wCx, const intvec *wCy,
2459  int& dx, int& dy, const ring r)
2460{
2461  const unsigned int N  = r->N;
2462
2463  p_Test(m, r);
2464
2465  assume( wx != NULL );
2466  assume( wy != NULL );
2467
2468  assume( wx->cols() == 1 );
2469  assume( wy->cols() == 1 );
2470
2471  assume( (unsigned int)wx->rows() >= N );
2472  assume( (unsigned int)wy->rows() >= N );
2473
2474  int x = 0;
2475  int y = 0;
2476
2477  for(int i = N; i > 0; i--)
2478  {
2479    const int d = p_GetExp(m, i, r);
2480    x += d * (*wx)[i-1];
2481    y += d * (*wy)[i-1];
2482  }
2483
2484  if( (wCx != NULL) && (wCy != NULL) )
2485  {
2486    const int c = p_GetComp(m, r);
2487
2488    if( wCx->range(c) )
2489      x += (*wCx)[c];
2490
2491    if( wCy->range(c) )
2492      x += (*wCy)[c];
2493  }
2494
2495  dx = x;
2496  dy = y;
2497}
2498
2499// returns true if polynom p is bi-homogenous with respect to the given weights
2500// simultaneously sets bi-Degree
2501bool p_IsBiHomogeneous(const poly p,
2502  const intvec *wx, const intvec *wy,
2503  const intvec *wCx, const intvec *wCy,
2504  int &dx, int &dy,
2505  const ring r)
2506{
2507  if( p == NULL )
2508  {
2509    dx = 0;
2510    dy = 0;
2511    return true;
2512  }
2513
2514  poly q = p;
2515
2516
2517  int ddx, ddy;
2518
2519  m_GetBiDegree( q, wx, wy, wCx, wCy, ddx, ddy, r); // get bi degree of lm(p)
2520
2521  pIter(q);
2522
2523  for(; q != NULL; pIter(q) )
2524  {
2525    int x, y;
2526
2527    m_GetBiDegree( q, wx, wy, wCx, wCy, x, y, r); // get bi degree of q
2528
2529    if ( (x != ddx) || (y != ddy) ) return false;
2530  }
2531
2532  dx = ddx;
2533  dy = ddy;
2534
2535  return true;
2536}
2537
2538
2539// returns true if id is bi-homogenous without respect to the given weights
2540bool id_IsBiHomogeneous(const ideal id,
2541  const intvec *wx, const intvec *wy,
2542  const intvec *wCx, const intvec *wCy,
2543  const ring r)
2544{
2545  if (id == NULL) return true; // zero ideal
2546
2547  const int iSize = IDELEMS(id);
2548
2549  if (iSize == 0) return true;
2550
2551  bool b = true;
2552  int x, y;
2553
2554  for(int i = iSize - 1; (i >= 0 ) && b; i--)
2555    b = p_IsBiHomogeneous(id->m[i], wx, wy, wCx, wCy, x, y, r);
2556
2557  return b;
2558}
2559
2560
2561// returns an intvector with [nvars(r)] integers [1/0]
2562// 1 - for commutative variables
2563// 0 - for anticommutative variables
2564intvec *ivGetSCAXVarWeights(const ring r)
2565{
2566  const unsigned int N  = r->N;
2567
2568  const int CommutativeVariable = 0; // bug correction!
2569  const int AntiCommutativeVariable = 0;
2570
2571  intvec* w = new intvec(N, 1, CommutativeVariable);
2572
2573  if(AntiCommutativeVariable != CommutativeVariable)
2574  if( rIsSCA(r) )
2575  {
2576    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2577    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
2578
2579    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2580    {
2581      (*w)[i-1] = AntiCommutativeVariable;
2582    }
2583  }
2584
2585  return w;
2586}
2587
2588
2589// returns an intvector with [nvars(r)] integers [1/0]
2590// 0 - for commutative variables
2591// 1 - for anticommutative variables
2592intvec *ivGetSCAYVarWeights(const ring r)
2593{
2594  const unsigned int N  = r->N;
2595
2596  const int CommutativeVariable = 0;
2597  const int AntiCommutativeVariable = 1;
2598
2599  intvec* w = new intvec(N, 1, CommutativeVariable);
2600
2601  if(AntiCommutativeVariable != CommutativeVariable)
2602  if( rIsSCA(r) )
2603  {
2604    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2605    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
2606
2607    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2608    {
2609      (*w)[i-1] = AntiCommutativeVariable;
2610    }
2611  }
2612  return w;
2613}
2614
2615
2616
2617
2618// reduce term lt(m) modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar:
2619// either create a copy of m or return NULL
2620inline poly m_KillSquares(const poly m,
2621  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2622  const ring r)
2623{
2624#ifdef PDEBUG
2625  p_Test(m, r);
2626  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2627
2628#if 0
2629  Print("m_KillSquares, m = "); // !
2630  p_Write(m, r);
2631#endif
2632#endif
2633
2634  assume( m != NULL );
2635
2636  for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
2637    if( p_GetExp(m, k, r) > 1 )
2638      return NULL;
2639
2640  return p_Head(m, r);
2641}
2642
2643
2644// reduce polynomial p modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2645poly p_KillSquares(const poly p,
2646  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2647  const ring r)
2648{
2649#ifdef PDEBUG
2650  p_Test(p, r);
2651
2652  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2653
2654#if 0
2655  Print("p_KillSquares, p = "); // !
2656  p_Write(p, r);
2657#endif
2658#endif
2659
2660
2661  if( p == NULL )
2662    return NULL;
2663
2664  poly pResult = NULL;
2665  poly* ppPrev = &pResult;
2666
2667  for( poly q = p; q!= NULL; pIter(q) )
2668  {
2669#ifdef PDEBUG
2670    p_Test(q, r);
2671#endif
2672
2673    // terms will be in the same order because of quasi-ordering!
2674    poly v = m_KillSquares(q, iFirstAltVar, iLastAltVar, r);
2675
2676    if( v != NULL )
2677    {
2678      *ppPrev = v;
2679      ppPrev = &pNext(v);
2680    }
2681
2682  }
2683
2684#ifdef PDEBUG
2685  p_Test(pResult, r);
2686#if 0
2687  Print("p_KillSquares => "); // !
2688  p_Write(pResult, r);
2689#endif
2690#endif
2691
2692  return(pResult);
2693}
2694
2695
2696
2697
2698// reduces ideal id modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2699// returns the reduced ideal or zero ideal.
2700ideal id_KillSquares(const ideal id,
2701  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2702  const ring r)
2703{
2704  if (id == NULL) return id; // zero ideal
2705
2706  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2707
2708  const int iSize = IDELEMS(id);
2709
2710  if (iSize == 0) return id;
2711
2712  ideal temp = idInit(iSize, id->rank);
2713
2714#if 0
2715   PrintS("<id_KillSquares>\n");
2716  {
2717    PrintS("ideal id: \n");
2718    for (int i = 0; i < IDELEMS(id); i++)
2719    {
2720      Print("; id[%d] = ", i+1);
2721      p_Write(id->m[i], r);
2722    }
2723    PrintS(";\n");
2724    PrintLn();
2725  }
2726#endif
2727
2728
2729  for (int j = 0; j < iSize; j++)
2730    temp->m[j] = p_KillSquares(id->m[j], iFirstAltVar, iLastAltVar, r);
2731
2732  idSkipZeroes(temp);
2733
2734#if 0
2735   PrintS("<id_KillSquares>\n");
2736  {
2737    PrintS("ideal temp: \n");
2738    for (int i = 0; i < IDELEMS(temp); i++)
2739    {
2740      Print("; temp[%d] = ", i+1);
2741      p_Write(temp->m[i], r);
2742    }
2743    PrintS(";\n");
2744    PrintLn();
2745  }
2746   PrintS("</id_KillSquares>\n");
2747#endif
2748
2749//  temp->rank = idRankFreeModule(temp, r);
2750
2751  return temp;
2752}
2753
2754
2755
2756
2757#endif
Note: See TracBrowser for help on using the repository browser.