source: git/kernel/GBEngine/kspoly.cc @ 70a107

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