source: git/kernel/GBEngine/kspoly.cc @ e009ef

fieker-DuValspielwiese
Last change on this file since e009ef was e40da9f, checked in by Adi Popescu <adi_popescum@…>, 8 years ago
trying to fix Christian sba algorithm
  • Property mode set to 100644
File size: 27.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5*  ABSTRACT -  Routines for Spoly creation and reductions
6*/
7
8// #define PDEBUG 2
9
10
11
12#include <kernel/mod2.h>
13#include <misc/options.h>
14#include <kernel/GBEngine/kutil.h>
15#include <coeffs/numbers.h>
16#include <polys/monomials/p_polys.h>
17#include <polys/templates/p_Procs.h>
18#include <polys/nc/nc.h>
19#ifdef KDEBUG
20#endif
21#ifdef HAVE_RINGS
22#include <kernel/polys.h>
23#endif
24
25//#define ADIDEBUG 0
26
27#ifdef KDEBUG
28int red_count = 0;
29int create_count = 0;
30// define this if reductions are reported on TEST_OPT_DEBUG
31#define TEST_OPT_DEBUG_RED
32#endif
33
34/***************************************************************
35 *
36 * Reduces PR with PW
37 * Assumes PR != NULL, PW != NULL, Lm(PW) divides Lm(PR)
38 *
39 ***************************************************************/
40int ksReducePoly(LObject* PR,
41                 TObject* PW,
42                 poly spNoether,
43                 number *coef,
44                 kStrategy strat)
45{
46#ifdef KDEBUG
47  red_count++;
48#ifdef TEST_OPT_DEBUG_RED
49  if (TEST_OPT_DEBUG)
50  {
51    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
52    PW->wrp();
53    //printf("\necart(PR)-ecart(PW): %i\n",PR->ecart-PW->ecart);
54    //pWrite(PR->p);
55  }
56#endif
57#endif
58  int ret = 0;
59  ring tailRing = PR->tailRing;
60  kTest_L(PR);
61  kTest_T(PW);
62
63  poly p1 = PR->GetLmTailRing();   // p2 | p1
64  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
65  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
66  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
67  p_CheckPolyRing(p1, tailRing);
68  p_CheckPolyRing(p2, tailRing);
69
70  pAssume1(p2 != NULL && p1 != NULL &&
71           p_DivisibleBy(p2,  p1, tailRing));
72
73  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
74           (p_GetComp(p2, tailRing) == 0 &&
75            p_MaxComp(pNext(p2),tailRing) == 0));
76
77#ifdef HAVE_PLURAL
78  if (rIsPluralRing(currRing))
79  {
80    // for the time being: we know currRing==strat->tailRing
81    // no exp-bound checking needed
82    // (only needed if exp-bound(tailring)<exp-b(currRing))
83    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
84    else
85    {
86      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
87      assume(_p != NULL);
88      nc_PolyPolyRed(_p, p2,coef, currRing);
89      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
90      PR->pLength=0; // usually not used, GetpLength re-computes it if needed
91    }
92    return 0;
93  }
94#endif
95
96  if (t2==NULL)           // Divisor is just one term, therefore it will
97  {                       // just cancel the leading term
98    PR->LmDeleteAndIter();
99    if (coef != NULL) *coef = n_Init(1, tailRing);
100    return 0;
101  }
102
103  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
104
105  if (tailRing != currRing)
106  {
107    // check that reduction does not violate exp bound
108    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
109    {
110      // undo changes of lm
111      p_ExpVectorAdd(lm, p2, tailRing);
112      if (strat == NULL) return 2;
113      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
114      tailRing = strat->tailRing;
115      p1 = PR->GetLmTailRing();
116      p2 = PW->GetLmTailRing();
117      t2 = pNext(p2);
118      lm = p1;
119      p_ExpVectorSub(lm, p2, tailRing);
120      ret = 1;
121    }
122  }
123
124  // take care of coef buisness
125  if (! n_IsOne(pGetCoeff(p2), tailRing))
126  {
127    number bn = pGetCoeff(lm);
128    number an = pGetCoeff(p2);
129    int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
130    p_SetCoeff(lm, bn, tailRing);
131    if ((ct == 0) || (ct == 2))
132      PR->Tail_Mult_nn(an);
133    if (coef != NULL) *coef = an;
134    else n_Delete(&an, tailRing);
135  }
136  else
137  {
138    if (coef != NULL) *coef = n_Init(1, tailRing);
139  }
140
141
142  // and finally,
143  PR->Tail_Minus_mm_Mult_qq(lm, t2, pLength(t2) /*PW->GetpLength() - 1*/, spNoether);
144  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
145  PR->LmDeleteAndIter();
146
147  // the following is commented out: shrinking
148#ifdef HAVE_SHIFTBBA_NONEXISTENT
149  if ( (currRing->isLPring) && (!strat->homog) )
150  {
151    // assume? h->p in currRing
152    PR->GetP();
153    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
154    PR->Clear(); // does the right things
155    PR->p = qq;
156    PR->t_p = NULL;
157    PR->SetShortExpVector();
158  }
159#endif
160
161#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
162  if (TEST_OPT_DEBUG)
163  {
164    Print(" to: "); PR->wrp(); Print("\n");
165    //printf("\nt^%i ", PR->ecart);pWrite(pHead(PR->p));
166  }
167#endif
168  return ret;
169}
170
171/***************************************************************
172 *
173 * Reduces PR with PW
174 * Assumes PR != NULL, PW != NULL, Lm(PW) divides Lm(PR)
175 *
176 ***************************************************************/
177 
178int ksReducePolySig(LObject* PR,
179                 TObject* PW,
180                 long /*idx*/,
181                 poly spNoether,
182                 number *coef,
183                 kStrategy strat)
184{
185#ifdef KDEBUG
186  red_count++;
187#ifdef TEST_OPT_DEBUG_RED
188  if (TEST_OPT_DEBUG)
189  {
190    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
191    PW->wrp();
192  }
193#endif
194#endif
195  int ret = 0;
196  ring tailRing = PR->tailRing;
197  kTest_L(PR);
198  kTest_T(PW);
199
200  // signature-based stuff:
201  // checking for sig-safeness first
202  // NOTE: This has to be done in the current ring
203  //
204  /**********************************************
205   *
206   * TODO:
207   * --------------------------------------------
208   * if strat->sbaOrder == 1
209   * Since we are subdividing lower index and
210   * current index reductions it is enough to
211   * look at the polynomial part of the signature
212   * for a check. This should speed-up checking
213   * a lot!
214   * if !strat->sbaOrder == 0
215   * We are not subdividing lower and current index
216   * due to the fact that we are using the induced
217   * Schreyer order
218   *
219   * nevertheless, this different behaviour is
220   * taken care of by is_sigsafe
221   * => one reduction procedure can be used for
222   * both, the incremental and the non-incremental
223   * attempt!
224   * --------------------------------------------
225   *
226   *********************************************/
227  //printf("COMPARE IDX: %ld -- %ld\n",idx,strat->currIdx);
228  if (!PW->is_sigsafe)
229  {
230    poly sigMult = pCopy(PW->sig);   // copy signature of reducer
231//#if 1
232#ifdef DEBUGF5
233    printf("IN KSREDUCEPOLYSIG: \n");
234    pWrite(pHead(f1));
235    pWrite(pHead(f2));
236    pWrite(sigMult);
237    printf("--------------\n");
238#endif
239    p_ExpVectorAddSub(sigMult,PR->GetLmCurrRing(),PW->GetLmCurrRing(),currRing);
240//#if 1
241#ifdef DEBUGF5
242    printf("------------------- IN KSREDUCEPOLYSIG: --------------------\n");
243    pWrite(pHead(f1));
244    pWrite(pHead(f2));
245    pWrite(sigMult);
246    pWrite(PR->sig);
247    printf("--------------\n");
248#endif
249    int sigSafe = p_LmCmp(PR->sig,sigMult,currRing);
250    // now we can delete the copied polynomial data used for checking for
251    // sig-safeness of the reduction step
252//#if 1
253#ifdef DEBUGF5
254    printf("%d -- %d sig\n",sigSafe,PW->is_sigsafe);
255
256#endif
257    //pDelete(&f1);
258    pDelete(&sigMult);
259    // go on with the computations only if the signature of p2 is greater than the
260    // signature of fm*p1
261    if(sigSafe != 1)
262    {
263      PR->is_redundant = TRUE;
264      return 3;
265    }
266    //PW->is_sigsafe  = TRUE;
267  }
268  PR->is_redundant = FALSE;
269  poly p1 = PR->GetLmTailRing();   // p2 | p1
270  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
271  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
272  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
273  p_CheckPolyRing(p1, tailRing);
274  p_CheckPolyRing(p2, tailRing);
275
276  pAssume1(p2 != NULL && p1 != NULL &&
277      p_DivisibleBy(p2,  p1, tailRing));
278
279  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
280      (p_GetComp(p2, tailRing) == 0 &&
281       p_MaxComp(pNext(p2),tailRing) == 0));
282
283#ifdef HAVE_PLURAL
284  if (rIsPluralRing(currRing))
285  {
286    // for the time being: we know currRing==strat->tailRing
287    // no exp-bound checking needed
288    // (only needed if exp-bound(tailring)<exp-b(currRing))
289    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
290    else
291    {
292      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
293      assume(_p != NULL);
294      nc_PolyPolyRed(_p, p2, coef, currRing);
295      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
296      PR->pLength=0; // usaully not used, GetpLength re-comoutes it if needed
297    }
298    return 0;
299  }
300#endif
301
302  if (t2==NULL)           // Divisor is just one term, therefore it will
303  {                       // just cancel the leading term
304    PR->LmDeleteAndIter();
305    if (coef != NULL) *coef = n_Init(1, tailRing);
306    return 0;
307  }
308
309  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
310
311  if (tailRing != currRing)
312  {
313    // check that reduction does not violate exp bound
314    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
315    {
316      // undo changes of lm
317      p_ExpVectorAdd(lm, p2, tailRing);
318      if (strat == NULL) return 2;
319      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
320      tailRing = strat->tailRing;
321      p1 = PR->GetLmTailRing();
322      p2 = PW->GetLmTailRing();
323      t2 = pNext(p2);
324      lm = p1;
325      p_ExpVectorSub(lm, p2, tailRing);
326      ret = 1;
327    }
328  }
329
330  // take care of coef buisness
331  if (! n_IsOne(pGetCoeff(p2), tailRing))
332  {
333    number bn = pGetCoeff(lm);
334    number an = pGetCoeff(p2);
335    int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
336    p_SetCoeff(lm, bn, tailRing);
337    if ((ct == 0) || (ct == 2))
338      PR->Tail_Mult_nn(an);
339    if (coef != NULL) *coef = an;
340    else n_Delete(&an, tailRing);
341  }
342  else
343  {
344    if (coef != NULL) *coef = n_Init(1, tailRing);
345  }
346
347
348  // and finally,
349  PR->Tail_Minus_mm_Mult_qq(lm, t2, PW->GetpLength() - 1, spNoether);
350  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
351  PR->LmDeleteAndIter();
352
353  // the following is commented out: shrinking
354#ifdef HAVE_SHIFTBBA_NONEXISTENT
355  if ( (currRing->isLPring) && (!strat->homog) )
356  {
357    // assume? h->p in currRing
358    PR->GetP();
359    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
360    PR->Clear(); // does the right things
361    PR->p = qq;
362    PR->t_p = NULL;
363    PR->SetShortExpVector();
364  }
365#endif
366
367#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
368  if (TEST_OPT_DEBUG)
369  {
370    Print(" to: "); PR->wrp(); Print("\n");
371  }
372#endif
373  return ret;
374}
375 
376int ksReducePolySigRing(LObject* PR,
377                 TObject* PW,
378                 long /*idx*/,
379                 poly spNoether,
380                 number *coef,
381                 kStrategy strat)
382{
383#ifdef ADIDEBUG
384printf("\nksReducePolySig\n");
385pWrite(PR->p);pWrite(PR->sig);
386pWrite(PW->p);pWrite(PW->sig);
387#endif
388#ifdef KDEBUG
389  red_count++;
390#ifdef TEST_OPT_DEBUG_RED
391  if (TEST_OPT_DEBUG)
392  {
393    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
394    PW->wrp();
395  }
396#endif
397#endif
398  int ret = 0;
399  ring tailRing = PR->tailRing;
400  kTest_L(PR);
401  kTest_T(PW);
402
403  // signature-based stuff:
404  // checking for sig-safeness first
405  // NOTE: This has to be done in the current ring
406  //
407  /**********************************************
408   *
409   * TODO:
410   * --------------------------------------------
411   * if strat->sbaOrder == 1
412   * Since we are subdividing lower index and
413   * current index reductions it is enough to
414   * look at the polynomial part of the signature
415   * for a check. This should speed-up checking
416   * a lot!
417   * if !strat->sbaOrder == 0
418   * We are not subdividing lower and current index
419   * due to the fact that we are using the induced
420   * Schreyer order
421   *
422   * nevertheless, this different behaviour is
423   * taken care of by is_sigsafe
424   * => one reduction procedure can be used for
425   * both, the incremental and the non-incremental
426   * attempt!
427   * --------------------------------------------
428   *
429   *********************************************/
430  //printf("COMPARE IDX: %ld -- %ld\n",idx,strat->currIdx);
431  if (!PW->is_sigsafe)
432  {
433    poly sigMult = pCopy(PW->sig);   // copy signature of reducer
434//#if 1
435#ifdef DEBUGF5
436    printf("IN KSREDUCEPOLYSIG: \n");
437    pWrite(pHead(f1));
438    pWrite(pHead(f2));
439    pWrite(sigMult);
440    printf("--------------\n");
441#endif
442    p_ExpVectorAddSub(sigMult,PR->GetLmCurrRing(),PW->GetLmCurrRing(),currRing);
443    //I have also to set the leading coeficient for sigMult (in the case of rings)
444    if(rField_is_Ring(currRing))
445    {
446      pSetCoeff(sigMult,nMult(nDiv(pGetCoeff(PR->p),pGetCoeff(PW->p)), pGetCoeff(sigMult)));
447      if(nIsZero(pGetCoeff(sigMult)))
448      {
449        sigMult = NULL;
450      }
451    }
452//#if 1
453#ifdef DEBUGF5
454    printf("------------------- IN KSREDUCEPOLYSIG: --------------------\n");
455    pWrite(pHead(f1));
456    pWrite(pHead(f2));
457    pWrite(sigMult);
458    pWrite(PR->sig);
459    printf("--------------\n");
460#endif
461    int sigSafe;
462    if(!rField_is_Ring(currRing))
463      sigSafe = p_LmCmp(PR->sig,sigMult,currRing);
464    // now we can delete the copied polynomial data used for checking for
465    // sig-safeness of the reduction step
466//#if 1
467#ifdef DEBUGF5
468    printf("%d -- %d sig\n",sigSafe,PW->is_sigsafe);
469
470#endif
471    if(rField_is_Ring(currRing))
472    {
473      // Set the sig
474      poly origsig = pCopy(PR->sig);
475      if(sigMult != NULL)
476        PR->sig = pHead(pSub(PR->sig, sigMult));
477      //The sigs have the same lm, have to substract
478      //It may happen that now the signature is 0 (drop)
479      if(PR->sig == NULL)
480      {
481        #ifdef ADIDEBUG
482        printf("\nPossible sigdrop in ksreducepolysig (lost signature)\n");
483        #endif
484        strat->sigdrop=TRUE;
485      }
486      else
487      {
488        if(pLtCmp(PR->sig,origsig) == 1)
489        {
490          // do not allow this reduction - it will increase it's signature
491          // and the partially standard basis is just till the old sig, not the new one
492          PR->is_redundant = TRUE;
493          pDelete(&PR->sig);
494          PR->sig = origsig;
495          strat->blockred++;
496          return 3;
497        }
498        if(pLtCmp(PR->sig,origsig) == -1)
499        {
500          #ifdef ADIDEBUG
501          printf("\nSigdrop in ksreducepolysig from * to *\n");pWrite(origsig);pWrite(PR->sig);
502          #endif
503          strat->sigdrop=TRUE;
504        }
505      }
506      pDelete(&origsig);
507    }
508    //pDelete(&f1);
509    // go on with the computations only if the signature of p2 is greater than the
510    // signature of fm*p1
511    if(sigSafe != 1 && !rField_is_Ring(currRing))
512    {
513      PR->is_redundant = TRUE;
514      return 3;
515    }
516    //PW->is_sigsafe  = TRUE;
517  }
518  PR->is_redundant = FALSE;
519  poly p1 = PR->GetLmTailRing();   // p2 | p1
520  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
521  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
522  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
523  p_CheckPolyRing(p1, tailRing);
524  p_CheckPolyRing(p2, tailRing);
525
526  pAssume1(p2 != NULL && p1 != NULL &&
527      p_DivisibleBy(p2,  p1, tailRing));
528
529  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
530      (p_GetComp(p2, tailRing) == 0 &&
531       p_MaxComp(pNext(p2),tailRing) == 0));
532
533#ifdef HAVE_PLURAL
534  if (rIsPluralRing(currRing))
535  {
536    // for the time being: we know currRing==strat->tailRing
537    // no exp-bound checking needed
538    // (only needed if exp-bound(tailring)<exp-b(currRing))
539    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
540    else
541    {
542      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
543      assume(_p != NULL);
544      nc_PolyPolyRed(_p, p2, coef, currRing);
545      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
546      PR->pLength=0; // usaully not used, GetpLength re-comoutes it if needed
547    }
548    return 0;
549  }
550#endif
551
552  if (t2==NULL)           // Divisor is just one term, therefore it will
553  {                       // just cancel the leading term
554    PR->LmDeleteAndIter();
555    if (coef != NULL) *coef = n_Init(1, tailRing);
556    return 0;
557  }
558
559  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
560
561  if (tailRing != currRing)
562  {
563    // check that reduction does not violate exp bound
564    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
565    {
566      // undo changes of lm
567      p_ExpVectorAdd(lm, p2, tailRing);
568      if (strat == NULL) return 2;
569      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
570      tailRing = strat->tailRing;
571      p1 = PR->GetLmTailRing();
572      p2 = PW->GetLmTailRing();
573      t2 = pNext(p2);
574      lm = p1;
575      p_ExpVectorSub(lm, p2, tailRing);
576      ret = 1;
577    }
578  }
579  // take care of coef buisness
580  if(rField_is_Ring(currRing))
581  {
582    p_SetCoeff(lm, nDiv(pGetCoeff(lm),pGetCoeff(p2)), tailRing);
583    if (coef != NULL) *coef = n_Init(1, tailRing);
584  }
585  else
586  {
587    if (! n_IsOne(pGetCoeff(p2), tailRing))
588    {
589      number bn = pGetCoeff(lm);
590      number an = pGetCoeff(p2);
591      int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
592      p_SetCoeff(lm, bn, tailRing);
593      if (((ct == 0) || (ct == 2)))
594        PR->Tail_Mult_nn(an);
595      if (coef != NULL) *coef = an;
596      else n_Delete(&an, tailRing);
597    }
598    else
599    {
600      if (coef != NULL) *coef = n_Init(1, tailRing);
601    }
602  }
603
604  // and finally,
605  PR->Tail_Minus_mm_Mult_qq(lm, t2, PW->GetpLength() - 1, spNoether);
606  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
607  PR->LmDeleteAndIter();
608
609  // the following is commented out: shrinking
610#ifdef HAVE_SHIFTBBA_NONEXISTENT
611  if ( (currRing->isLPring) && (!strat->homog) )
612  {
613    // assume? h->p in currRing
614    PR->GetP();
615    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
616    PR->Clear(); // does the right things
617    PR->p = qq;
618    PR->t_p = NULL;
619    PR->SetShortExpVector();
620  }
621#endif
622#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
623  if (TEST_OPT_DEBUG)
624  {
625    Print(" to: "); PR->wrp(); Print("\n");
626  }
627#endif
628  return ret;
629}
630
631/***************************************************************
632 *
633 * Creates S-Poly of p1 and p2
634 *
635 *
636 ***************************************************************/
637void ksCreateSpoly(LObject* Pair,   poly spNoether,
638                   int use_buckets, ring tailRing,
639                   poly m1, poly m2, TObject** R)
640{
641#ifdef KDEBUG
642  create_count++;
643#endif
644  kTest_L(Pair);
645  poly p1 = Pair->p1;
646  poly p2 = Pair->p2;
647  Pair->tailRing = tailRing;
648
649  assume(p1 != NULL);
650  assume(p2 != NULL);
651  assume(tailRing != NULL);
652
653  poly a1 = pNext(p1), a2 = pNext(p2);
654  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
655  int co=0/*, ct = ksCheckCoeff(&lc1, &lc2, currRing->cf)*/; // gcd and zero divisors
656  (void) ksCheckCoeff(&lc1, &lc2, currRing->cf);
657
658  int l1=0, l2=0;
659
660  if (p_GetComp(p1, currRing)!=p_GetComp(p2, currRing))
661  {
662    if (p_GetComp(p1, currRing)==0)
663    {
664      co=1;
665      p_SetCompP(p1,p_GetComp(p2, currRing), currRing, tailRing);
666    }
667    else
668    {
669      co=2;
670      p_SetCompP(p2, p_GetComp(p1, currRing), currRing, tailRing);
671    }
672  }
673
674  // get m1 = LCM(LM(p1), LM(p2))/LM(p1)
675  //     m2 = LCM(LM(p1), LM(p2))/LM(p2)
676  if (m1 == NULL)
677    k_GetLeadTerms(p1, p2, currRing, m1, m2, tailRing);
678
679  pSetCoeff0(m1, lc2);
680  pSetCoeff0(m2, lc1);  // and now, m1 * LT(p1) == m2 * LT(p2)
681
682  if (R != NULL)
683  {
684    if (Pair->i_r1 == -1)
685    {
686      l1 = pLength(p1) - 1;
687    }
688    else
689    {
690      l1 = (R[Pair->i_r1])->GetpLength() - 1;
691    }
692    if ((Pair->i_r2 == -1)||(R[Pair->i_r2]==NULL))
693    {
694      l2 = pLength(p2) - 1;
695    }
696    else
697    {
698      l2 = (R[Pair->i_r2])->GetpLength() - 1;
699    }
700  }
701
702  // get m2 * a2
703  if (spNoether != NULL)
704  {
705    l2 = -1;
706    a2 = tailRing->p_Procs->pp_Mult_mm_Noether(a2, m2, spNoether, l2, tailRing);
707    assume(l2 == pLength(a2));
708  }
709  else
710    a2 = tailRing->p_Procs->pp_Mult_mm(a2, m2, tailRing);
711#ifdef HAVE_RINGS
712  if (!(rField_is_Domain(currRing))) l2 = pLength(a2);
713#endif
714
715  Pair->SetLmTail(m2, a2, l2, use_buckets, tailRing);
716
717  // get m2*a2 - m1*a1
718  Pair->Tail_Minus_mm_Mult_qq(m1, a1, l1, spNoether);
719
720  // Clean-up time
721  Pair->LmDeleteAndIter();
722  p_LmDelete(m1, tailRing);
723
724  if (co != 0)
725  {
726    if (co==1)
727    {
728      p_SetCompP(p1,0, currRing, tailRing);
729    }
730    else
731    {
732      p_SetCompP(p2,0, currRing, tailRing);
733    }
734  }
735
736  // the following is commented out: shrinking
737#ifdef HAVE_SHIFTBBA_NONEXISTENT
738  if (currRing->isLPring)
739  {
740    // assume? h->p in currRing
741    Pair->GetP();
742    poly qq = p_Shrink(Pair->p, currRing->isLPring, currRing);
743    Pair->Clear(); // does the right things
744    Pair->p = qq;
745    Pair->t_p = NULL;
746    Pair->SetShortExpVector();
747  }
748#endif
749
750}
751
752int ksReducePolyTail(LObject* PR, TObject* PW, poly Current, poly spNoether)
753{
754  BOOLEAN ret;
755  number coef;
756  poly Lp =     PR->GetLmCurrRing();
757  poly Save =   PW->GetLmCurrRing();
758
759  kTest_L(PR);
760  kTest_T(PW);
761  pAssume(pIsMonomOf(Lp, Current));
762
763  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
764  assume(PR->bucket == NULL);
765
766  LObject Red(pNext(Current), PR->tailRing);
767  TObject With(PW, Lp == Save);
768
769  pAssume(!pHaveCommonMonoms(Red.p, With.p));
770  ret = ksReducePoly(&Red, &With, spNoether, &coef);
771
772  if (!ret)
773  {
774    if (! n_IsOne(coef, currRing))
775    {
776      pNext(Current) = NULL;
777      if (Current == PR->p && PR->t_p != NULL)
778        pNext(PR->t_p) = NULL;
779      PR->Mult_nn(coef);
780    }
781
782    n_Delete(&coef, currRing);
783    pNext(Current) = Red.GetLmTailRing();
784    if (Current == PR->p && PR->t_p != NULL)
785      pNext(PR->t_p) = pNext(Current);
786  }
787
788  if (Lp == Save)
789    With.Delete();
790
791  // the following is commented out: shrinking
792#ifdef HAVE_SHIFTBBA_NONEXISTENT
793  if (currRing->isLPring)
794  {
795    // assume? h->p in currRing
796    PR->GetP();
797    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
798    PR->Clear(); // does the right things
799    PR->p = qq;
800    PR->t_p = NULL;
801    PR->SetShortExpVector();
802  }
803#endif
804
805  return ret;
806}
807
808/***************************************************************
809 *
810 * Auxillary Routines
811 *
812 *
813 ***************************************************************/
814
815/*2
816* creates the leading term of the S-polynomial of p1 and p2
817* do not destroy p1 and p2
818* remarks:
819*   1. the coefficient is 0 (nNew)
820*   1. a) in the case of coefficient ring, the coefficient is calculated
821*   2. pNext is undefined
822*/
823//static void bbb() { int i=0; }
824poly ksCreateShortSpoly(poly p1, poly p2, ring tailRing)
825{
826  poly a1 = pNext(p1), a2 = pNext(p2);
827  long c1=p_GetComp(p1, currRing),c2=p_GetComp(p2, currRing);
828  long c;
829  poly m1,m2;
830  number t1 = NULL,t2 = NULL;
831  int cm,i;
832  BOOLEAN equal;
833
834#ifdef HAVE_RINGS
835  BOOLEAN is_Ring=rField_is_Ring(currRing);
836  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
837  if (is_Ring)
838  {
839    ksCheckCoeff(&lc1, &lc2, currRing->cf); // gcd and zero divisors
840    if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
841    if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
842    while (a1 != NULL && nIsZero(t2))
843    {
844      pIter(a1);
845      nDelete(&t2);
846      if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
847    }
848    while (a2 != NULL && nIsZero(t1))
849    {
850      pIter(a2);
851      nDelete(&t1);
852      if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
853    }
854  }
855#endif
856
857  if (a1==NULL)
858  {
859    if(a2!=NULL)
860    {
861      m2=p_Init(currRing);
862x2:
863      for (i = (currRing->N); i; i--)
864      {
865        c = p_GetExpDiff(p1, p2,i, currRing);
866        if (c>0)
867        {
868          p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)),currRing);
869        }
870        else
871        {
872          p_SetExp(m2,i,p_GetExp(a2,i,tailRing),currRing);
873        }
874      }
875      if ((c1==c2)||(c2!=0))
876      {
877        p_SetComp(m2,p_GetComp(a2,tailRing), currRing);
878      }
879      else
880      {
881        p_SetComp(m2,c1,currRing);
882      }
883      p_Setm(m2, currRing);
884#ifdef HAVE_RINGS
885      if (is_Ring)
886      {
887          nDelete(&lc1);
888          nDelete(&lc2);
889          nDelete(&t2);
890          pSetCoeff0(m2, t1);
891      }
892      else
893#endif
894        nNew(&(pGetCoeff(m2)));
895      return m2;
896    }
897    else
898    {
899#ifdef HAVE_RINGS
900      if (is_Ring)
901      {
902        nDelete(&lc1);
903        nDelete(&lc2);
904        nDelete(&t1);
905        nDelete(&t2);
906      }
907#endif
908      return NULL;
909    }
910  }
911  if (a2==NULL)
912  {
913    m1=p_Init(currRing);
914x1:
915    for (i = (currRing->N); i; i--)
916    {
917      c = p_GetExpDiff(p2, p1,i,currRing);
918      if (c>0)
919      {
920        p_SetExp(m1,i,(c+p_GetExp(a1,i, tailRing)),currRing);
921      }
922      else
923      {
924        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
925      }
926    }
927    if ((c1==c2)||(c1!=0))
928    {
929      p_SetComp(m1,p_GetComp(a1,tailRing),currRing);
930    }
931    else
932    {
933      p_SetComp(m1,c2,currRing);
934    }
935    p_Setm(m1, currRing);
936#ifdef HAVE_RINGS
937    if (is_Ring)
938    {
939      pSetCoeff0(m1, t2);
940      nDelete(&lc1);
941      nDelete(&lc2);
942      nDelete(&t1);
943    }
944    else
945#endif
946      nNew(&(pGetCoeff(m1)));
947    return m1;
948  }
949  m1 = p_Init(currRing);
950  m2 = p_Init(currRing);
951  loop
952  {
953    for (i = (currRing->N); i; i--)
954    {
955      c = p_GetExpDiff(p1, p2,i,currRing);
956      if (c > 0)
957      {
958        p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)), currRing);
959        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
960      }
961      else
962      {
963        p_SetExp(m1,i,(p_GetExp(a1,i,tailRing)-c), currRing);
964        p_SetExp(m2,i,p_GetExp(a2,i, tailRing), currRing);
965      }
966    }
967    if(c1==c2)
968    {
969      p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
970      p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
971    }
972    else
973    {
974      if(c1!=0)
975      {
976        p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
977        p_SetComp(m2,c1, currRing);
978      }
979      else
980      {
981        p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
982        p_SetComp(m1,c2, currRing);
983      }
984    }
985    p_Setm(m1,currRing);
986    p_Setm(m2,currRing);
987    cm = p_LmCmp(m1, m2,currRing);
988    if (cm!=0)
989    {
990      if(cm==1)
991      {
992        p_LmFree(m2,currRing);
993#ifdef HAVE_RINGS
994        if (is_Ring)
995        {
996          pSetCoeff0(m1, t2);
997          nDelete(&lc1);
998          nDelete(&lc2);
999          nDelete(&t1);
1000        }
1001        else
1002#endif
1003          nNew(&(pGetCoeff(m1)));
1004        return m1;
1005      }
1006      else
1007      {
1008        p_LmFree(m1,currRing);
1009#ifdef HAVE_RINGS
1010        if (is_Ring)
1011        {
1012          pSetCoeff0(m2, t1);
1013          nDelete(&lc1);
1014          nDelete(&lc2);
1015          nDelete(&t2);
1016        }
1017        else
1018#endif
1019          nNew(&(pGetCoeff(m2)));
1020        return m2;
1021      }
1022    }
1023#ifdef HAVE_RINGS
1024    if (is_Ring)
1025    {
1026      equal = nEqual(t1,t2);
1027    }
1028    else
1029#endif
1030    {
1031      t1 = nMult(pGetCoeff(a2),pGetCoeff(p1));
1032      t2 = nMult(pGetCoeff(a1),pGetCoeff(p2));
1033      equal = nEqual(t1,t2);
1034      nDelete(&t2);
1035      nDelete(&t1);
1036    }
1037    if (!equal)
1038    {
1039      p_LmFree(m2,currRing);
1040#ifdef HAVE_RINGS
1041      if (is_Ring)
1042      {
1043          pSetCoeff0(m1, nSub(t1, t2));
1044          nDelete(&lc1);
1045          nDelete(&lc2);
1046          nDelete(&t1);
1047          nDelete(&t2);
1048      }
1049      else
1050#endif
1051        nNew(&(pGetCoeff(m1)));
1052      return m1;
1053    }
1054    pIter(a1);
1055    pIter(a2);
1056#ifdef HAVE_RINGS
1057    if (is_Ring)
1058    {
1059      if (a2 != NULL)
1060      {
1061        nDelete(&t1);
1062        t1 = nMult(pGetCoeff(a2),lc1);
1063      }
1064      if (a1 != NULL)
1065      {
1066        nDelete(&t2);
1067        t2 = nMult(pGetCoeff(a1),lc2);
1068      }
1069      while ((a1 != NULL) && nIsZero(t2))
1070      {
1071        pIter(a1);
1072        if (a1 != NULL)
1073        {
1074          nDelete(&t2);
1075          t2 = nMult(pGetCoeff(a1),lc2);
1076        }
1077      }
1078      while ((a2 != NULL) && nIsZero(t1))
1079      {
1080        pIter(a2);
1081        if (a2 != NULL)
1082        {
1083          nDelete(&t1);
1084          t1 = nMult(pGetCoeff(a2),lc1);
1085        }
1086      }
1087    }
1088#endif
1089    if (a2==NULL)
1090    {
1091      p_LmFree(m2,currRing);
1092      if (a1==NULL)
1093      {
1094#ifdef HAVE_RINGS
1095        if (is_Ring)
1096        {
1097          nDelete(&lc1);
1098          nDelete(&lc2);
1099          nDelete(&t1);
1100          nDelete(&t2);
1101        }
1102#endif
1103        p_LmFree(m1,currRing);
1104        return NULL;
1105      }
1106      goto x1;
1107    }
1108    if (a1==NULL)
1109    {
1110      p_LmFree(m1,currRing);
1111      goto x2;
1112    }
1113  }
1114}
Note: See TracBrowser for help on using the repository browser.