source: git/kernel/sca.cc @ 098f98f

spielwiese
Last change on this file since 098f98f was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 64.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    sca.cc
6 *  Purpose: supercommutative kernel procedures
7 *  Author:  motsak (Oleksandr Motsak)
8 *  Created: 2006/12/18
9 *  Version: $Id$
10 *******************************************************************/
11
12// set it here if needed.
13#define OUTPUT 0
14#define MYTEST 0
15
16#if MYTEST
17#define OM_CHECK 4
18#define OM_TRACK 5
19#endif
20
21// #define PDEBUG 2
22#include "mod2.h"
23
24#ifdef HAVE_PLURAL
25// for
26#define PLURAL_INTERNAL_DECLARATIONS
27#include "sca.h"
28#include "gring.h"
29
30
31#include "febase.h"
32
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 Algebra extension by Oleksandr
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    register unsigned int tpower = 0;
115    register unsigned int cpower = 0;
116
117    for( register 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 ) // TODO: think about eliminating there if-s...
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    register unsigned int tpower = 0;
166    register unsigned int cpower = 0;
167
168    for( register 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    register unsigned int tpower = 0;
230    register unsigned int cpower = 0;
231
232    for( register 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    register unsigned int tpower = 0;
295    register unsigned int cpower = 0;
296
297    for( register 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    register unsigned int cpower = 0;
367
368    for( register 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 ((defined(PDEBUG) && OUTPUT) || 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 nc_struct* NC = rG->GetNC();
1363  const ring rBase = NC->basering;
1364  const matrix C   = NC->C; // live in rBase!
1365  const matrix D   = NC->D; // live in rBase!
1366
1367#if OUTPUT
1368  PrintS("sca_SetupQuotient: AltVars?!\n");
1369#endif
1370
1371  for(int i = 1; i < N; i++)
1372  {
1373    for(int j = i + 1; j <= N; j++)
1374    {
1375      if( MATELEM(D,i,j) != NULL) // !!!???
1376      {
1377#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1378        Print("Nonzero D[%d, %d]\n", i, j);
1379#endif
1380        return false;
1381      }
1382     
1383     
1384      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1385      number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1386
1387      if( n_IsMOne(c, rBase) ) // !!!???
1388      {
1389        if( i < iAltVarStart)
1390          iAltVarStart = i;
1391
1392        if( j > iAltVarEnd)
1393          iAltVarEnd = j;
1394      } else
1395      {
1396        if( !n_IsOne(c, rBase) )
1397        {
1398#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1399          Print("Wrong Coeff at: [%d, %d]\n", i, j);
1400#endif
1401          return false;
1402        }
1403      }
1404    }
1405  }
1406
1407#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1408  Print("AltVars?1: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1409#endif
1410
1411
1412  if( (iAltVarEnd == -1) || (iAltVarStart == (N+1)) )
1413    return false; // either no alternating varables, or a single one => we are in commutative case!
1414
1415
1416  for(int i = 1; i < N; i++)
1417  {
1418    for(int j = i + 1; j <= N; j++)
1419    {
1420      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1421      number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1422
1423      if( (iAltVarStart <= i) && (j <= iAltVarEnd) ) // S <= i < j <= E
1424      { // anticommutative part
1425        if( !n_IsMOne(c, rBase) )
1426        {
1427#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1428           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1429#endif
1430          return false;
1431        }
1432      } else
1433      { // should commute
1434        if( !n_IsOne(c, rBase) )
1435        {
1436#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1437           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1438#endif
1439          return false;
1440        }
1441      }
1442    }
1443  }
1444
1445#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1446  Print("AltVars!?: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1447#endif
1448
1449  assume( 1            <= iAltVarStart );
1450  assume( iAltVarStart < iAltVarEnd   );
1451  assume( iAltVarEnd   <= N            );
1452
1453
1454  ring rSaveRing = assureCurrentRing(rG);
1455
1456
1457  assume(rGR->qideal != NULL);
1458  assume(rGR->N == rG->N);
1459//  assume(rG->qideal == NULL); // ?
1460
1461  const ideal idQuotient = rGR->qideal;
1462
1463
1464#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1465  PrintS("Analyzing quotient ideal:\n");
1466  idPrint(idQuotient); // in rG!!!
1467#endif
1468
1469
1470  // check for
1471  // y_{iAltVarStart}^2, y_{iAltVarStart+1}^2, \ldots, y_{iAltVarEnd}^2  (iAltVarEnd > iAltVarStart)
1472  // to be within quotient ideal.
1473
1474  bool bSCA = true;
1475
1476  int b = N+1;
1477  int e = -1; 
1478
1479  if(rIsSCA(rG))
1480  {
1481    b = si_min(b, scaFirstAltVar(rG));
1482    e = si_max(e, scaLastAltVar(rG));
1483
1484#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1485    Print("AltVars!?: [%d, %d]\n", b, e);
1486#endif
1487   
1488  }
1489
1490 
1491  for ( int i = iAltVarStart; (i <= iAltVarEnd) && bSCA; i++ )
1492    if( (i < b) || (i > e) ) // otherwise it's ok since rG is an SCA!
1493  {
1494    poly square = p_One( rG);
1495    p_SetExp(square, i, 2, rG); // square = var(i)^2.
1496    p_Setm(square, rG);
1497
1498    // square = NF( var(i)^2 | Q )
1499    // NOTE: rSaveRing == currRing now!
1500    // NOTE: there is no better way to check this in general!
1501    square = kNF(idQuotient, NULL, square, 0, 1); // must ran in currRing == rG!
1502
1503    if( square != NULL ) // var(i)^2 is not in Q?
1504    {
1505      p_Delete(&square, rG);
1506      bSCA = false;
1507      break;
1508    }
1509  }
1510
1511  assureCurrentRing(rSaveRing);
1512 
1513  if(!bSCA) return false;
1514
1515
1516#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1517  Print("ScaVars!: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1518#endif
1519
1520
1521  //////////////////////////////////////////////////////////////////////////
1522  // ok... here we go. let's setup it!!!
1523  //////////////////////////////////////////////////////////////////////////
1524  ideal tempQ = id_KillSquares(idQuotient, iAltVarStart, iAltVarEnd, rG); // in rG!!!
1525
1526
1527#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1528  PrintS("Quotient: \n");
1529  iiWriteMatrix((matrix)idQuotient,"__",1);
1530  PrintS("tempSCAQuotient: \n");
1531  iiWriteMatrix((matrix)tempQ,"__",1);
1532#endif
1533 
1534  idSkipZeroes( tempQ );
1535
1536  ncRingType( rGR, nc_exterior );
1537
1538  scaFirstAltVar( rGR, iAltVarStart );
1539  scaLastAltVar( rGR, iAltVarEnd );
1540 
1541  if( idIs0(tempQ) )
1542    rGR->GetNC()->SCAQuotient() = NULL;
1543  else
1544    rGR->GetNC()->SCAQuotient() = idrMoveR(tempQ, rG, rGR); // deletes tempQ!
1545
1546  nc_p_ProcsSet(rGR, rGR->p_Procs); // !!!!!!!!!!!!!!!!!
1547
1548
1549#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1550  PrintS("SCAQuotient: \n");
1551  if(tempQ != NULL)
1552    iiWriteMatrix((matrix)tempQ,"__",1);
1553  else
1554    PrintS("(NULL)\n");
1555#endif
1556 
1557  return true;
1558}
1559
1560
1561bool sca_Force(ring rGR, int b, int e)
1562{
1563  assume(rGR != NULL);
1564  assume(rIsPluralRing(rGR));
1565  assume(!rIsSCA(rGR));
1566
1567  const int N = rGR->N;
1568
1569  ring rSaveRing = currRing;
1570
1571  if(rSaveRing != rGR)
1572    rChangeCurrRing(rGR);
1573
1574  const ideal idQuotient = rGR->qideal;
1575
1576  ideal tempQ = idQuotient;
1577
1578  if( b <= N && e >= 1 )
1579    tempQ = id_KillSquares(idQuotient, b, e, rGR);
1580
1581  idSkipZeroes( tempQ );
1582
1583  ncRingType( rGR, nc_exterior );
1584
1585  if( idIs0(tempQ) )
1586    rGR->GetNC()->SCAQuotient() = NULL;
1587  else
1588    rGR->GetNC()->SCAQuotient() = tempQ;
1589
1590 
1591  scaFirstAltVar( rGR, b );
1592  scaLastAltVar( rGR, e );
1593
1594
1595  nc_p_ProcsSet(rGR, rGR->p_Procs);
1596
1597  if(rSaveRing != rGR)
1598    rChangeCurrRing(rSaveRing);
1599
1600  return true;
1601}
1602
1603// return x_i * pPoly; preserve pPoly.
1604poly sca_pp_Mult_xi_pp(unsigned int i, const poly pPoly, const ring rRing)
1605{
1606  assume(1 <= i);
1607  assume(i <= (unsigned int)rRing->N);
1608
1609  if(rIsSCA(rRing))
1610    return sca_xi_Mult_pp(i, pPoly, rRing);
1611
1612
1613
1614  poly xi =  p_One( rRing);
1615  p_SetExp(xi, i, 1, rRing);
1616  p_Setm(xi, rRing);
1617
1618  poly pResult = pp_Mult_qq(xi, pPoly, rRing);
1619
1620  p_Delete( &xi, rRing);
1621
1622  return pResult;
1623
1624}
1625
1626
1627///////////////////////////////////////////////////////////////////////////////////////
1628//************* SCA BBA *************************************************************//
1629///////////////////////////////////////////////////////////////////////////////////////
1630
1631// Under development!!!
1632ideal sca_bba (const ideal F, const ideal Q, const intvec *w, const intvec * /*hilb*/, kStrategy strat)
1633{
1634#if MYTEST
1635  PrintS("\n\n<sca_bba>\n\n");
1636#endif
1637
1638  assume(rIsSCA(currRing));
1639
1640#ifndef NDEBUG
1641  idTest(F);
1642  idTest(Q);
1643#endif
1644
1645#if MYTEST
1646  PrintS("\ncurrRing: \n");
1647  rWrite(currRing);
1648#ifdef RDEBUG
1649//  rDebugPrint(currRing);
1650#endif
1651
1652  PrintS("\n\nF: \n");
1653  idPrint(F);
1654  PrintS("\n\nQ: \n");
1655  idPrint(Q);
1656
1657  PrintLn();
1658#endif
1659 
1660
1661  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
1662  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
1663
1664  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
1665
1666  ideal tempQ = Q;
1667
1668  if(Q == currQuotient)
1669    tempQ = SCAQuotient(currRing);
1670
1671  // Q or tempQ will not be used below :(((
1672
1673
1674#if MYTEST
1675
1676  PrintS("tempF: \n");
1677  idPrint(tempF);
1678  PrintS("tempQ: \n");
1679  idPrint(tempQ);
1680#endif
1681 
1682  strat->z2homog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
1683   // redo no_prod_crit:
1684  const BOOLEAN bIsSCA  = rIsSCA(currRing) && strat->z2homog; // for Z_2 prod-crit
1685  strat->no_prod_crit   = ! bIsSCA;
1686
1687//  strat->homog = strat->homog && strat->z2homog; // ?
1688
1689
1690 
1691#ifdef KDEBUG
1692  om_Opts.MinTrack = 5;
1693#endif
1694 
1695  int   srmax, lrmax, red_result = 1;
1696  int   olddeg, reduc;
1697
1698//  int hilbeledeg = 1, minimcnt = 0;
1699  int hilbcount = 0;
1700
1701  BOOLEAN withT = FALSE;
1702
1703  initBuchMoraCrit(strat); // sets Gebauer, honey, sugarCrit // sca - ok???
1704  initBuchMoraPos(strat); // sets strat->posInL, strat->posInT // check!! (Plural's: )
1705
1706//   initHilbCrit(F, Q, &hilb, strat);
1707
1708//  nc_gr_initBba(F,strat);
1709  initBba(tempF, strat); // set enterS, red, initEcart, initEcartPair
1710
1711  /*set enterS, spSpolyShort, reduce, red, initEcart, initEcartPair*/
1712  // ?? set spSpolyShort, reduce ???
1713  initBuchMora(tempF, tempQ, strat); // tempQ = Q without squares!!!
1714
1715//   if (strat->minim>0) strat->M = idInit(IDELEMS(F),F->rank);
1716
1717  srmax = strat->sl;
1718  reduc = olddeg = lrmax = 0;
1719
1720#define NO_BUCKETS
1721
1722#ifndef NO_BUCKETS
1723  if (!TEST_OPT_NOT_BUCKETS)
1724    strat->use_buckets = 1;
1725#endif
1726
1727  // redtailBBa against T for inhomogenous input
1728  if (!K_TEST_OPT_OLDSTD)
1729    withT = ! strat->homog;
1730
1731  // strat->posInT = posInT_pLength;
1732  kTest_TS(strat);
1733
1734#undef HAVE_TAIL_RING
1735
1736#ifdef HAVE_TAIL_RING
1737  kStratInitChangeTailRing(strat);
1738#endif
1739
1740  ///////////////////////////////////////////////////////////////
1741  // SCA:
1742
1743  //  due to std( SB, p).
1744  // Note that after initBuchMora :: initSSpecial all these additional
1745  // elements are in S and T (and some pairs are in L, which also has no initiall
1746  // elements!!!)
1747  if(TEST_OPT_SB_1)
1748  {
1749    // For all additional elements...
1750    for (int iNewElement = strat->newIdeal; iNewElement < IDELEMS(tempF); iNewElement++)
1751    {
1752      const poly pSave = tempF->m[iNewElement];
1753
1754      if( pSave != NULL )
1755      {
1756//        tempF->m[iNewElement] = NULL;
1757
1758        const poly p_next = pNext(pSave);
1759
1760        if(p_next != NULL)
1761          for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
1762            if( p_GetExp(pSave, i, currRing) != 0 )
1763            {
1764              assume(p_GetExp(pSave, i, currRing) == 1);
1765
1766              const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
1767
1768#ifdef PDEBUG
1769              p_Test(p_new, currRing);
1770#endif
1771
1772              if( p_new == NULL) continue;
1773
1774              LObject h(p_new); // h = x_i * strat->P
1775
1776              if (TEST_OPT_INTSTRATEGY)
1777                h.pCleardenom(); // also does a pContent
1778              else
1779                h.pNorm();
1780
1781              strat->initEcart(&h);
1782              h.sev = pGetShortExpVector(h.p);
1783
1784              int pos = 0;
1785
1786              if (strat->Ll != -1)
1787                pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1788
1789              enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
1790            }         
1791
1792      }       
1793
1794    }
1795  } 
1796
1797 
1798 
1799
1800  /* compute------------------------------------------------------- */
1801  while (strat->Ll >= 0)
1802  {
1803    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
1804
1805#ifdef KDEBUG
1806//     loop_count++;
1807    if (TEST_OPT_DEBUG) messageSets(strat);
1808#endif
1809
1810    if (strat->Ll== 0) strat->interpt=TRUE;
1811
1812    if (TEST_OPT_DEGBOUND
1813        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1814            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
1815    {
1816      /*
1817       *stops computation if
1818       * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
1819       *a predefined number Kstd1_deg
1820       */
1821      while ((strat->Ll >= 0)
1822  && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
1823        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1824            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))
1825  )
1826        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1827      if (strat->Ll<0) break;
1828      else strat->noClearS=TRUE;
1829    }
1830
1831    /* picks the last element from the lazyset L */
1832    strat->P = strat->L[strat->Ll];
1833    strat->Ll--;
1834
1835
1836//    assume(pNext(strat->P.p) != strat->tail);
1837
1838    if(strat->P.IsNull()) continue;
1839
1840    if (pNext(strat->P.p) == strat->tail)
1841    {
1842      // deletes the short spoly
1843      pLmFree(strat->P.p);
1844
1845      strat->P.p = nc_CreateSpoly(strat->P.p1, strat->P.p2, currRing);
1846      if (strat->P.p!=NULL) strat->initEcart(&strat->P);
1847    }//    else
1848
1849
1850    if(strat->P.IsNull()) continue;
1851
1852    if (strat->P.p1 == NULL)
1853    {
1854//       if (strat->minim > 0)
1855//         strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
1856
1857
1858      // for input polys, prepare reduction
1859        strat->P.PrepareRed(strat->use_buckets);
1860    }
1861
1862    if (TEST_OPT_PROT)
1863      message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
1864              &olddeg,&reduc,strat, red_result);
1865
1866    /* reduction of the element choosen from L */
1867    red_result = strat->red(&strat->P,strat);
1868
1869
1870    // reduction to non-zero new poly
1871    if (red_result == 1)
1872    {
1873      /* statistic */
1874      if (TEST_OPT_PROT) PrintS("s");
1875
1876      // get the polynomial (canonicalize bucket, make sure P.p is set)
1877      strat->P.GetP(strat->lmBin);
1878
1879      int pos = posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
1880
1881      // reduce the tail and normalize poly
1882      if (TEST_OPT_INTSTRATEGY)
1883      {
1884        strat->P.pCleardenom();
1885        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1886        {
1887          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT); // !!!
1888          strat->P.pCleardenom();
1889        }
1890      }
1891      else
1892      {
1893        strat->P.pNorm();
1894        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1895          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
1896      }
1897      strat->P.is_normalized=nIsOne(pGetCoeff(strat->P.p));
1898
1899#ifdef KDEBUG
1900      if (TEST_OPT_DEBUG){PrintS(" ns:");p_wrp(strat->P.p,currRing);PrintLn();}
1901#endif
1902
1903//       // min_std stuff
1904//       if ((strat->P.p1==NULL) && (strat->minim>0))
1905//       {
1906//         if (strat->minim==1)
1907//         {
1908//           strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
1909//           p_Delete(&strat->P.p2, currRing, strat->tailRing);
1910//         }
1911//         else
1912//         {
1913//           strat->M->m[minimcnt]=strat->P.p2;
1914//           strat->P.p2=NULL;
1915//         }
1916//         if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
1917//           pNext(strat->M->m[minimcnt])
1918//             = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
1919//                                            strat->tailRing, currRing,
1920//                                            currRing->PolyBin);
1921//         minimcnt++;
1922//       }
1923
1924      // enter into S, L, and T
1925      //if(withT)
1926        enterT(strat->P, strat);
1927
1928      // L
1929      enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
1930
1931      // posInS only depends on the leading term
1932      strat->enterS(strat->P, pos, strat, strat->tl);
1933
1934//       if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
1935
1936//      Print("[%d]",hilbeledeg);
1937      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
1938
1939      if (strat->sl>srmax) srmax = strat->sl;
1940
1941      // //////////////////////////////////////////////////////////
1942      // SCA:
1943      const poly pSave = strat->P.p;
1944      const poly p_next = pNext(pSave);
1945
1946//       if(0)
1947      if( p_next != NULL )
1948      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
1949      if( p_GetExp(pSave, i, currRing) != 0 )
1950      {
1951        assume(p_GetExp(pSave, i, currRing) == 1);
1952        const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
1953
1954#ifdef PDEBUG
1955        p_Test(p_new, currRing);
1956#endif
1957
1958        if( p_new == NULL) continue;
1959
1960        LObject h(p_new); // h = x_i * strat->P
1961
1962        if (TEST_OPT_INTSTRATEGY)
1963        {
1964//          pContent(h.p);
1965          h.pCleardenom(); // also does a pContent
1966        }
1967        else
1968        {
1969          h.pNorm();
1970        }
1971
1972        strat->initEcart(&h);
1973        h.sev = pGetShortExpVector(h.p);
1974
1975        int pos = 0;
1976
1977        if (strat->Ll != -1)
1978          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
1979
1980        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
1981/*
1982        h.sev = pGetShortExpVector(h.p);
1983        strat->initEcart(&h);
1984
1985        h.PrepareRed(strat->use_buckets);
1986
1987        // reduction of the element choosen from L(?)
1988        red_result = strat->red(&h,strat);
1989
1990        // reduction to non-zero new poly
1991        if (red_result != 1) continue;
1992
1993
1994        int pos = posInS(strat,strat->sl,h.p,h.ecart);
1995
1996        // reduce the tail and normalize poly
1997        if (TEST_OPT_INTSTRATEGY)
1998        {
1999          h.pCleardenom();
2000          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
2001          {
2002            h.p = redtailBba(&(h),pos-1,strat, withT); // !!!
2003            h.pCleardenom();
2004          }
2005        }
2006        else
2007        {
2008          h.pNorm();
2009          if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
2010            h.p = redtailBba(&(h),pos-1,strat, withT);
2011        }
2012
2013#ifdef KDEBUG
2014        if (TEST_OPT_DEBUG){PrintS(" N:");h.wrp();PrintLn();}
2015#endif
2016
2017//        h.PrepareRed(strat->use_buckets); // ???
2018
2019        h.sev = pGetShortExpVector(h.p);
2020        strat->initEcart(&h);
2021
2022        if (strat->Ll==-1)
2023          pos = 0;
2024        else
2025          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
2026
2027         enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);*/
2028
2029      } // for all x_i \in Ann(lm(P))
2030    } // if red(P) != NULL
2031
2032//     else if (strat->P.p1 == NULL && strat->minim > 0)
2033//     {
2034//       p_Delete(&strat->P.p2, currRing, strat->tailRing);
2035//     }
2036
2037
2038   
2039#ifdef KDEBUG
2040//    memset(&(strat->P), 0, sizeof(strat->P));
2041#endif
2042
2043    kTest_TS(strat); // even of T is not used!
2044
2045//     Print("\n$\n");
2046
2047  }
2048
2049#ifdef KDEBUG
2050  if (TEST_OPT_DEBUG) messageSets(strat);
2051#endif
2052
2053  /* complete reduction of the standard basis--------- */
2054
2055  if (TEST_OPT_REDSB){
2056    completeReduce(strat);
2057  }
2058
2059 
2060  /* release temp data-------------------------------- */
2061
2062  exitBuchMora(strat); // cleanT!
2063
2064  id_Delete(&tempF, currRing);
2065 
2066  if (TEST_OPT_WEIGHTM)
2067  {
2068    pRestoreDegProcs(pFDegOld, pLDegOld);
2069    if (ecartWeights)
2070    {
2071      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
2072      ecartWeights=NULL;
2073    }
2074  }
2075
2076  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
2077
2078
2079 
2080  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
2081
2082
2083  if (TEST_OPT_REDSB) // ???
2084  {
2085    // must be at the very end (after exitBuchMora) as it changes the S set!!!
2086    ideal I = strat->Shdl;
2087    ideal erg = kInterRedOld(I,tempQ);
2088    assume(I!=erg);
2089    id_Delete(&I, currRing);
2090    strat->Shdl = erg;
2091  }
2092 
2093 
2094#if MYTEST
2095  PrintS("\n\n</sca_bba>\n\n");
2096#endif
2097 
2098  return (strat->Shdl);
2099}
2100
2101
2102// //////////////////////////////////////////////////////////////////////////////
2103// sca mora...
2104
2105// returns TRUE if mora should use buckets, false otherwise
2106static BOOLEAN kMoraUseBucket(kStrategy strat)
2107{
2108#ifdef MORA_USE_BUCKETS
2109  if (TEST_OPT_NOT_BUCKETS)
2110    return FALSE;
2111  if (strat->red == redFirst)
2112  {
2113#ifdef NO_LDEG
2114    if (!strat->syzComp)
2115      return TRUE;
2116#else
2117    if ((strat->homog || strat->honey) && !strat->syzComp)
2118      return TRUE;
2119#endif
2120  }
2121  else
2122  {
2123    assume(strat->red == redEcart);
2124    if (strat->honey && !strat->syzComp)
2125      return TRUE;
2126  }
2127#endif
2128  return FALSE;
2129}
2130
2131
2132#ifdef HAVE_ASSUME
2133static int sca_mora_count = 0;
2134static int sca_mora_loop_count;
2135#endif
2136
2137// ideal sca_mora (ideal F, ideal Q, intvec *w, intvec *, kStrategy strat)
2138ideal sca_mora(const ideal F, const ideal Q, const intvec *w, const intvec *, kStrategy strat)
2139{
2140  assume(rIsSCA(currRing));
2141
2142  const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
2143  const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
2144
2145  ideal tempF = id_KillSquares(F, m_iFirstAltVar, m_iLastAltVar, currRing);
2146
2147  ideal tempQ = Q;
2148
2149  if(Q == currQuotient)
2150    tempQ = SCAQuotient(currRing);
2151
2152  bool bIdHomog = id_IsSCAHomogeneous(tempF, NULL, NULL, currRing); // wCx == wCy == NULL!
2153
2154  assume( !bIdHomog || strat->homog ); //  bIdHomog =====[implies]>>>>> strat->homog
2155
2156  strat->homog = strat->homog && bIdHomog;
2157
2158#ifdef PDEBUG
2159  assume( strat->homog == bIdHomog );
2160#endif /*PDEBUG*/
2161
2162#ifdef HAVE_ASSUME
2163  sca_mora_count++;
2164  sca_mora_loop_count = 0;
2165#endif
2166
2167#ifdef KDEBUG
2168  om_Opts.MinTrack = 5;
2169#endif
2170
2171
2172  strat->update = TRUE;
2173  /*- setting global variables ------------------- -*/
2174  initBuchMoraCrit(strat);
2175//   initHilbCrit(F,NULL,&hilb,strat); // no Q!
2176  initMora(tempF, strat);
2177  initBuchMoraPos(strat);
2178  /*Shdl=*/initBuchMora(tempF, tempQ, strat); // temp Q, F!
2179//   if (TEST_OPT_FASTHC) missingAxis(&strat->lastAxis,strat);
2180  /*updateS in initBuchMora has Hecketest
2181  * and could have put strat->kHEdgdeFound FALSE*/
2182#if 0
2183  if (ppNoether!=NULL)
2184  {
2185    strat->kHEdgeFound = TRUE;
2186  }
2187  if (strat->kHEdgeFound && strat->update)
2188  {
2189    firstUpdate(strat);
2190    updateLHC(strat);
2191    reorderL(strat);
2192  }
2193  if (TEST_OPT_FASTHC && (strat->lastAxis) && strat->posInLOldFlag)
2194  {
2195    strat->posInLOld = strat->posInL;
2196    strat->posInLOldFlag = FALSE;
2197    strat->posInL = posInL10;
2198    updateL(strat);
2199    reorderL(strat);
2200  }
2201#endif
2202  strat->use_buckets = kMoraUseBucket(strat);
2203
2204  kTest_TS(strat);
2205
2206
2207  int srmax = strat->sl;
2208  int lrmax = strat->Ll;
2209  int olddeg = 0;
2210  int reduc = 0;
2211  int red_result = 1;
2212//  int hilbeledeg=1;
2213  int hilbcount=0;
2214
2215 
2216  /*- compute-------------------------------------------*/
2217
2218#undef HAVE_TAIL_RING
2219
2220#ifdef HAVE_TAIL_RING
2221//  if (strat->homog && strat->red == redFirst)
2222//     kStratInitChangeTailRing(strat);
2223#endif
2224
2225
2226
2227
2228 
2229//  due to std( SB, p)
2230  if(TEST_OPT_SB_1)
2231  {
2232    for (int iNewElement = strat->newIdeal; iNewElement < IDELEMS(tempF); iNewElement++)
2233    {
2234
2235      const poly pSave = tempF->m[iNewElement];
2236
2237      if( pSave != NULL )
2238      {
2239//        tempF->m[iNewElement] = NULL;
2240
2241        const poly p_next = pNext(pSave);
2242
2243        if(p_next != NULL)
2244          for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
2245            if( p_GetExp(pSave, i, currRing) != 0 )
2246            {
2247
2248              assume(p_GetExp(pSave, i, currRing) == 1);
2249
2250              const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
2251
2252#ifdef PDEBUG
2253              p_Test(p_new, currRing);
2254#endif
2255
2256              if( p_new == NULL) continue;
2257
2258              LObject h(p_new); // h = x_i * strat->P
2259
2260              if (TEST_OPT_INTSTRATEGY)
2261                h.pCleardenom(); // also does a pContent
2262              else
2263                h.pNorm();
2264
2265              strat->initEcart(&h);
2266              h.sev = pGetShortExpVector(h.p);
2267
2268              int pos = 0;
2269
2270              if (strat->Ll != -1)
2271                pos = strat->posInL(strat->L,strat->Ll,&h,strat);
2272
2273              enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
2274
2275              if (strat->Ll>lrmax) lrmax = strat->Ll;
2276            }         
2277     
2278      }     
2279
2280    }
2281  } 
2282
2283
2284
2285
2286  while (strat->Ll >= 0)
2287  {
2288#ifdef HAVE_ASSUME
2289    sca_mora_loop_count++;
2290#endif
2291    if (lrmax< strat->Ll) lrmax=strat->Ll; /*stat*/
2292    //test_int_std(strat->kIdeal);
2293#ifdef KDEBUG
2294    if (TEST_OPT_DEBUG) messageSets(strat);
2295#endif
2296    if (TEST_OPT_DEGBOUND
2297    && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg))
2298    {
2299      /*
2300      * stops computation if
2301      * - 24 (degBound)
2302      *   && upper degree is bigger than Kstd1_deg
2303      */
2304      while ((strat->Ll >= 0)
2305        && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
2306        && (strat->L[strat->Ll].ecart+strat->L[strat->Ll].GetpFDeg()> Kstd1_deg)
2307      )
2308      {
2309        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
2310        //if (TEST_OPT_PROT)
2311        //{
2312        //   PrintS("D"); mflush();
2313        //}
2314      }
2315      if (strat->Ll<0) break;
2316      else strat->noClearS=TRUE;
2317    }
2318    strat->P = strat->L[strat->Ll];/*- picks the last element from the lazyset L -*/
2319    if (strat->Ll==0) strat->interpt=TRUE;
2320    strat->Ll--;
2321
2322    // create the real Spoly
2323//    assume(pNext(strat->P.p) != strat->tail);
2324
2325    if(strat->P.IsNull()) continue;
2326
2327
2328    if( pNext(strat->P.p) == strat->tail )
2329    {
2330      // deletes the int spoly and computes SPoly
2331      pLmFree(strat->P.p); // ???
2332      strat->P.p = nc_CreateSpoly(strat->P.p1, strat->P.p2, currRing);
2333    }
2334
2335
2336
2337    if (strat->P.p1 == NULL)
2338    {
2339      // for input polys, prepare reduction (buckets !)
2340      strat->P.SetLength(strat->length_pLength);
2341      strat->P.PrepareRed(strat->use_buckets);
2342    }
2343
2344    if (!strat->P.IsNull())
2345    {
2346      // might be NULL from noether !!!
2347      if (TEST_OPT_PROT)
2348        message(strat->P.ecart+strat->P.GetpFDeg(),&olddeg,&reduc,strat, red_result);
2349      // reduce
2350      red_result = strat->red(&strat->P,strat);
2351    }
2352
2353    if (! strat->P.IsNull())
2354    {
2355      strat->P.GetP();
2356      // statistics
2357      if (TEST_OPT_PROT) PrintS("s");
2358      // normalization
2359      if (!TEST_OPT_INTSTRATEGY)
2360        strat->P.pNorm();
2361      // tailreduction
2362      strat->P.p = redtail(&(strat->P),strat->sl,strat);
2363      // set ecart -- might have changed because of tail reductions
2364      if ((!strat->noTailReduction) && (!strat->honey))
2365        strat->initEcart(&strat->P);
2366      // cancel unit
2367      cancelunit(&strat->P);
2368      // for char 0, clear denominators
2369      if (TEST_OPT_INTSTRATEGY)
2370        strat->P.pCleardenom();
2371
2372      // put in T
2373      enterT(strat->P,strat);
2374      // build new pairs
2375      enterpairs(strat->P.p,strat->sl,strat->P.ecart,0,strat, strat->tl);
2376      // put in S
2377      strat->enterS(strat->P,
2378                    posInS(strat,strat->sl,strat->P.p, strat->P.ecart),
2379                    strat, strat->tl);
2380
2381
2382      // clear strat->P
2383      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
2384      strat->P.lcm=NULL;
2385
2386      if (strat->sl>srmax) srmax = strat->sl; /*stat.*/
2387      if (strat->Ll>lrmax) lrmax = strat->Ll;
2388
2389
2390
2391      // //////////////////////////////////////////////////////////
2392      // SCA:
2393      const poly pSave = strat->P.p;
2394      const poly p_next = pNext(pSave);
2395
2396      if(p_next != NULL)
2397      for( unsigned int i = m_iFirstAltVar; i <= m_iLastAltVar; i++ )
2398      if( p_GetExp(pSave, i, currRing) != 0 )
2399      {
2400
2401        assume(p_GetExp(pSave, i, currRing) == 1);
2402
2403        const poly p_new = sca_pp_Mult_xi_pp(i, p_next, currRing);
2404
2405#ifdef PDEBUG
2406        p_Test(p_new, currRing);
2407#endif
2408
2409        if( p_new == NULL) continue;
2410
2411        LObject h(p_new); // h = x_i * strat->P
2412
2413        if (TEST_OPT_INTSTRATEGY)
2414           h.pCleardenom(); // also does a pContent
2415        else
2416          h.pNorm();
2417
2418        strat->initEcart(&h);
2419        h.sev = pGetShortExpVector(h.p);
2420
2421        int pos = 0;
2422
2423        if (strat->Ll != -1)
2424          pos = strat->posInL(strat->L,strat->Ll,&h,strat);
2425
2426        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
2427
2428        if (strat->Ll>lrmax) lrmax = strat->Ll;
2429      }
2430
2431#ifdef KDEBUG
2432      // make sure kTest_TS does not complain about strat->P
2433      memset(&strat->P,0,sizeof(strat->P));
2434#endif
2435    }
2436#if 0
2437    if (strat->kHEdgeFound)
2438    {
2439      if ((TEST_OPT_FINDET)
2440      || ((TEST_OPT_MULTBOUND) && (scMult0Int((strat->Shdl)) < mu)))
2441      {
2442        // obachman: is this still used ???
2443        /*
2444        * stops computation if strat->kHEdgeFound and
2445        * - 27 (finiteDeterminacyTest)
2446        * or
2447        * - 23
2448        *   (multBound)
2449        *   && multiplicity of the ideal is smaller then a predefined number mu
2450        */
2451        while (strat->Ll >= 0) deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
2452      }
2453    }
2454#endif
2455    kTest_TS(strat);
2456  }
2457  /*- complete reduction of the standard basis------------------------ -*/
2458  if (TEST_OPT_REDSB) completeReduce(strat);
2459  /*- release temp data------------------------------- -*/
2460  exitBuchMora(strat);
2461  /*- polynomials used for HECKE: HC, noether -*/
2462  if (TEST_OPT_FINDET)
2463  {
2464    if (strat->kHEdge!=NULL)
2465      Kstd1_mu=pFDeg(strat->kHEdge,currRing);
2466    else
2467      Kstd1_mu=-1;
2468  }
2469  pDelete(&strat->kHEdge);
2470  strat->update = TRUE; //???
2471  strat->lastAxis = 0; //???
2472  pDelete(&strat->kNoether);
2473  omFreeSize((ADDRESS)strat->NotUsedAxis,(pVariables+1)*sizeof(BOOLEAN));
2474  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
2475  if (TEST_OPT_WEIGHTM)
2476  {
2477    pRestoreDegProcs(pFDegOld, pLDegOld);
2478    if (ecartWeights)
2479    {
2480      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
2481      ecartWeights=NULL;
2482    }
2483  }
2484  if (tempQ!=NULL) updateResult(strat->Shdl,tempQ,strat);
2485  idTest(strat->Shdl);
2486
2487  id_Delete( &tempF, currRing);
2488
2489  return (strat->Shdl);
2490}
2491
2492
2493
2494
2495
2496
2497void sca_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
2498{
2499
2500  // "commutative" procedures:
2501  rGR->p_Procs->p_Mult_mm     = sca_p_Mult_mm;
2502  rGR->p_Procs->pp_Mult_mm    = sca_pp_Mult_mm;
2503
2504  p_Procs->p_Mult_mm          = sca_p_Mult_mm;
2505  p_Procs->pp_Mult_mm         = sca_pp_Mult_mm;
2506
2507  // non-commutaitve
2508  rGR->GetNC()->p_Procs.mm_Mult_p   = sca_mm_Mult_p;
2509  rGR->GetNC()->p_Procs.mm_Mult_pp  = sca_mm_Mult_pp;
2510
2511
2512  if (rHasLocalOrMixedOrdering(rGR))
2513  {
2514#ifdef PDEBUG
2515//           Print("Local case => GB == mora!\n");
2516#endif
2517    rGR->GetNC()->p_Procs.GB          = sca_mora; // local ordering => Mora, otherwise - Buchberger!
2518  }
2519  else
2520  {
2521#ifdef PDEBUG
2522//           Print("Global case => GB == bba!\n");
2523#endif
2524    rGR->GetNC()->p_Procs.GB          = sca_bba; // sca_gr_bba; // sca_bba? // sca_bba;
2525  }
2526
2527
2528//   rGR->GetNC()->p_Procs.GlobalGB    = sca_gr_bba;
2529//   rGR->GetNC()->p_Procs.LocalGB     = sca_mora;
2530
2531
2532//   rGR->GetNC()->p_Procs.SPoly         = sca_SPoly;
2533//   rGR->GetNC()->p_Procs.ReduceSPoly   = sca_ReduceSpoly;
2534
2535#if 0
2536
2537        // Multiplication procedures:
2538
2539        p_Procs->p_Mult_mm   = sca_p_Mult_mm;
2540        _p_procs->p_Mult_mm  = sca_p_Mult_mm;
2541
2542        p_Procs->pp_Mult_mm  = sca_pp_Mult_mm;
2543        _p_procs->pp_Mult_mm = sca_pp_Mult_mm;
2544
2545        r->GetNC()->mmMultP()     = sca_mm_Mult_p;
2546        r->GetNC()->mmMultPP()    = sca_mm_Mult_pp;
2547
2548        r->GetNC()->GB()            = sca_gr_bba;
2549/*
2550        // ??????????????????????????????????????? coefficients swell...
2551        r->GetNC()->SPoly()         = sca_SPoly;
2552        r->GetNC()->ReduceSPoly()   = sca_ReduceSpoly;
2553*/
2554//         r->GetNC()->BucketPolyRed() = gnc_kBucketPolyRed;
2555//         r->GetNC()->BucketPolyRed_Z()= gnc_kBucketPolyRed_Z;
2556
2557#endif
2558}
2559
2560
2561// bi-Degree (x, y) of monomial "m"
2562// x-es and y-s are weighted by wx and wy resp.
2563// [optional] components have weights by wCx and wCy.
2564inline void m_GetBiDegree(const poly m,
2565  const intvec *wx, const intvec *wy,
2566  const intvec *wCx, const intvec *wCy,
2567  int& dx, int& dy, const ring r)
2568{
2569  const unsigned int N  = r->N;
2570
2571  p_Test(m, r);
2572
2573  assume( wx != NULL );
2574  assume( wy != NULL );
2575
2576  assume( wx->cols() == 1 );
2577  assume( wy->cols() == 1 );
2578
2579  assume( (unsigned int)wx->rows() >= N );
2580  assume( (unsigned int)wy->rows() >= N );
2581
2582  int x = 0;
2583  int y = 0;
2584
2585  for(int i = N; i > 0; i--)
2586  {
2587    const int d = p_GetExp(m, i, r);
2588    x += d * (*wx)[i-1];
2589    y += d * (*wy)[i-1];
2590  }
2591
2592  if( (wCx != NULL) && (wCy != NULL) )
2593  {
2594    const int c = p_GetComp(m, r);
2595
2596    if( wCx->range(c) )
2597      x += (*wCx)[c];
2598
2599    if( wCy->range(c) )
2600      x += (*wCy)[c];
2601  }
2602
2603  dx = x;
2604  dy = y;
2605}
2606
2607// returns true if polynom p is bi-homogenous with respect to the given weights
2608// simultaneously sets bi-Degree
2609bool p_IsBiHomogeneous(const poly p,
2610  const intvec *wx, const intvec *wy,
2611  const intvec *wCx, const intvec *wCy,
2612  int &dx, int &dy,
2613  const ring r)
2614{
2615  if( p == NULL )
2616  {
2617    dx = 0;
2618    dy = 0;
2619    return true;
2620  }
2621
2622  poly q = p;
2623
2624
2625  int ddx, ddy;
2626
2627  m_GetBiDegree( q, wx, wy, wCx, wCy, ddx, ddy, r); // get bi degree of lm(p)
2628
2629  pIter(q);
2630
2631  for(; q != NULL; pIter(q) )
2632  {
2633    int x, y;
2634
2635    m_GetBiDegree( q, wx, wy, wCx, wCy, x, y, r); // get bi degree of q
2636
2637    if ( (x != ddx) || (y != ddy) ) return false;
2638  }
2639
2640  dx = ddx;
2641  dy = ddy;
2642
2643  return true;
2644}
2645
2646
2647// returns true if id is bi-homogenous without respect to the given weights
2648bool id_IsBiHomogeneous(const ideal id,
2649  const intvec *wx, const intvec *wy,
2650  const intvec *wCx, const intvec *wCy,
2651  const ring r)
2652{
2653  if (id == NULL) return true; // zero ideal
2654
2655  const int iSize = IDELEMS(id);
2656
2657  if (iSize == 0) return true;
2658
2659  bool b = true;
2660  int x, y;
2661
2662  for(int i = iSize - 1; (i >= 0 ) && b; i--)
2663    b = p_IsBiHomogeneous(id->m[i], wx, wy, wCx, wCy, x, y, r);
2664
2665  return b;
2666}
2667
2668
2669// returns an intvector with [nvars(r)] integers [1/0]
2670// 1 - for commutative variables
2671// 0 - for anticommutative variables
2672intvec *ivGetSCAXVarWeights(const ring r)
2673{
2674  const unsigned int N  = r->N;
2675
2676  const int CommutativeVariable = 0; // bug correction!
2677  const int AntiCommutativeVariable = 0;
2678
2679  intvec* w = new intvec(N, 1, CommutativeVariable);
2680
2681  if(AntiCommutativeVariable != CommutativeVariable)
2682  if( rIsSCA(r) )
2683  {
2684    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2685    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
2686
2687    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2688    {
2689      (*w)[i-1] = AntiCommutativeVariable;
2690    }
2691  }
2692
2693  return w;
2694}
2695
2696
2697// returns an intvector with [nvars(r)] integers [1/0]
2698// 0 - for commutative variables
2699// 1 - for anticommutative variables
2700intvec *ivGetSCAYVarWeights(const ring r)
2701{
2702  const unsigned int N  = r->N;
2703
2704  const int CommutativeVariable = 0;
2705  const int AntiCommutativeVariable = 1;
2706
2707  intvec* w = new intvec(N, 1, CommutativeVariable);
2708
2709  if(AntiCommutativeVariable != CommutativeVariable)
2710  if( rIsSCA(r) )
2711  {
2712    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
2713    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
2714
2715    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
2716    {
2717      (*w)[i-1] = AntiCommutativeVariable;
2718    }
2719  }
2720  return w;
2721}
2722
2723
2724
2725
2726// reduce term lt(m) modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar:
2727// either create a copy of m or return NULL
2728inline poly m_KillSquares(const poly m,
2729  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2730  const ring r)
2731{
2732#ifdef PDEBUG
2733  p_Test(m, r);
2734  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2735
2736#if 0
2737  Print("m_KillSquares, m = "); // !
2738  p_Write(m, r);
2739#endif
2740#endif
2741
2742  assume( m != NULL );
2743
2744  for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
2745    if( p_GetExp(m, k, r) > 1 )
2746      return NULL;
2747
2748  return p_Head(m, r);
2749}
2750
2751
2752// reduce polynomial p modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2753// returns a new poly!
2754poly p_KillSquares(const poly p,
2755  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2756  const ring r)
2757{
2758#ifdef PDEBUG
2759  p_Test(p, r);
2760
2761  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2762
2763#if 0
2764  Print("p_KillSquares, p = "); // !
2765  p_Write(p, r);
2766#endif
2767#endif
2768
2769
2770  if( p == NULL )
2771    return NULL;
2772
2773  poly pResult = NULL;
2774  poly* ppPrev = &pResult;
2775
2776  for( poly q = p; q!= NULL; pIter(q) )
2777  {
2778#ifdef PDEBUG
2779    p_Test(q, r);
2780#endif
2781
2782    // terms will be in the same order because of quasi-ordering!
2783    poly v = m_KillSquares(q, iFirstAltVar, iLastAltVar, r);
2784
2785    if( v != NULL )
2786    {
2787      *ppPrev = v;
2788      ppPrev = &pNext(v);
2789    }
2790
2791  }
2792
2793#ifdef PDEBUG
2794  p_Test(pResult, r);
2795#if 0
2796  Print("p_KillSquares => "); // !
2797  p_Write(pResult, r);
2798#endif
2799#endif
2800
2801  return(pResult);
2802}
2803
2804
2805
2806
2807// reduces ideal id modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
2808// returns the reduced ideal or zero ideal.
2809ideal id_KillSquares(const ideal id,
2810  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
2811  const ring r, const bool bSkipZeroes)
2812{
2813  if (id == NULL) return id; // zero ideal
2814
2815  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
2816
2817  const int iSize = IDELEMS(id);
2818
2819  if (iSize == 0) return id;
2820
2821  ideal temp = idInit(iSize, id->rank);
2822
2823#if 0
2824   PrintS("<id_KillSquares>\n");
2825  {
2826    PrintS("ideal id: \n");
2827    for (int i = 0; i < IDELEMS(id); i++)
2828    {
2829      Print("; id[%d] = ", i+1);
2830      p_Write(id->m[i], r);
2831    }
2832    PrintS(";\n");
2833    PrintLn();
2834  }
2835#endif
2836
2837
2838  for (int j = 0; j < iSize; j++)
2839    temp->m[j] = p_KillSquares(id->m[j], iFirstAltVar, iLastAltVar, r);
2840
2841  if( bSkipZeroes )
2842    idSkipZeroes(temp);
2843
2844#if 0
2845   PrintS("<id_KillSquares>\n");
2846  {
2847    PrintS("ideal temp: \n");
2848    for (int i = 0; i < IDELEMS(temp); i++)
2849    {
2850      Print("; temp[%d] = ", i+1);
2851      p_Write(temp->m[i], r);
2852    }
2853    PrintS(";\n");
2854    PrintLn();
2855  }
2856   PrintS("</id_KillSquares>\n");
2857#endif
2858
2859//  temp->rank = idRankFreeModule(temp, r);
2860
2861  return temp;
2862}
2863
2864
2865
2866
2867#endif
Note: See TracBrowser for help on using the repository browser.