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

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