source: git/kernel/kspoly.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was c91ffe3, checked in by Christian Eder, 9 years ago
fix for deep sba tail reduction
  • Property mode set to 100644
File size: 20.5 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#ifdef HAVE_CONFIG_H
10#include "singularconfig.h"
11#endif /* HAVE_CONFIG_H */
12#include <kernel/mod2.h>
13#include <misc/options.h>
14#include <kernel/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#include <kernel/febase.h>
21#endif
22#ifdef HAVE_RINGS
23#include <kernel/polys.h>
24#endif
25
26#ifdef KDEBUG
27int red_count = 0;
28int create_count = 0;
29// define this if reductions are reported on TEST_OPT_DEBUG
30#define TEST_OPT_DEBUG_RED
31#endif
32
33/***************************************************************
34 *
35 * Reduces PR with PW
36 * Assumes PR != NULL, PW != NULL, Lm(PW) divides Lm(PR)
37 *
38 ***************************************************************/
39int ksReducePoly(LObject* PR,
40                 TObject* PW,
41                 poly spNoether,
42                 number *coef,
43                 kStrategy strat)
44{
45#ifdef KDEBUG
46  red_count++;
47#ifdef TEST_OPT_DEBUG_RED
48  if (TEST_OPT_DEBUG)
49  {
50    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
51    PW->wrp();
52  }
53#endif
54#endif
55  int ret = 0;
56  ring tailRing = PR->tailRing;
57  assume(kTest_L(PR));
58  assume(kTest_T(PW));
59
60  poly p1 = PR->GetLmTailRing();   // p2 | p1
61  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
62  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
63  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
64  p_CheckPolyRing(p1, tailRing);
65  p_CheckPolyRing(p2, tailRing);
66
67  pAssume1(p2 != NULL && p1 != NULL &&
68           p_DivisibleBy(p2,  p1, tailRing));
69
70  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
71           (p_GetComp(p2, tailRing) == 0 &&
72            p_MaxComp(pNext(p2),tailRing) == 0));
73
74#ifdef HAVE_PLURAL
75  if (rIsPluralRing(currRing))
76  {
77    // for the time being: we know currRing==strat->tailRing
78    // no exp-bound checking needed
79    // (only needed if exp-bound(tailring)<exp-b(currRing))
80    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
81    else
82    {
83      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
84      assume(_p != NULL);
85      nc_PolyPolyRed(_p, p2,coef, currRing);
86      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
87      PR->pLength=0; // usually not used, GetpLength re-computes it if needed
88    }
89    return 0;
90  }
91#endif
92
93  if (t2==NULL)           // Divisor is just one term, therefore it will
94  {                       // just cancel the leading term
95    PR->LmDeleteAndIter();
96    if (coef != NULL) *coef = n_Init(1, tailRing);
97    return 0;
98  }
99
100  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
101
102  if (tailRing != currRing)
103  {
104    // check that reduction does not violate exp bound
105    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
106    {
107      // undo changes of lm
108      p_ExpVectorAdd(lm, p2, tailRing);
109      if (strat == NULL) return 2;
110      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
111      tailRing = strat->tailRing;
112      p1 = PR->GetLmTailRing();
113      p2 = PW->GetLmTailRing();
114      t2 = pNext(p2);
115      lm = p1;
116      p_ExpVectorSub(lm, p2, tailRing);
117      ret = 1;
118    }
119  }
120
121  // take care of coef buisness
122  if (! n_IsOne(pGetCoeff(p2), tailRing))
123  {
124    number bn = pGetCoeff(lm);
125    number an = pGetCoeff(p2);
126    int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
127    p_SetCoeff(lm, bn, tailRing);
128    if ((ct == 0) || (ct == 2))
129      PR->Tail_Mult_nn(an);
130    if (coef != NULL) *coef = an;
131    else n_Delete(&an, tailRing);
132  }
133  else
134  {
135    if (coef != NULL) *coef = n_Init(1, tailRing);
136  }
137
138
139  // and finally,
140  PR->Tail_Minus_mm_Mult_qq(lm, t2, PW->GetpLength() - 1, spNoether);
141  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
142  PR->LmDeleteAndIter();
143
144  // the following is commented out: shrinking
145#ifdef HAVE_SHIFTBBA_NONEXISTENT
146  if ( (currRing->isLPring) && (!strat->homog) )
147  {
148    // assume? h->p in currRing
149    PR->GetP();
150    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
151    PR->Clear(); // does the right things
152    PR->p = qq;
153    PR->t_p = NULL;
154    PR->SetShortExpVector();
155  }
156#endif
157
158#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
159  if (TEST_OPT_DEBUG)
160  {
161    Print(" to: "); PR->wrp(); Print("\n");
162  }
163#endif
164  return ret;
165}
166
167/***************************************************************
168 *
169 * Reduces PR with PW
170 * Assumes PR != NULL, PW != NULL, Lm(PW) divides Lm(PR)
171 *
172 ***************************************************************/
173int ksReducePolySig(LObject* PR,
174                 TObject* PW,
175                 long /*idx*/,
176                 poly spNoether,
177                 number *coef,
178                 kStrategy strat)
179{
180#ifdef KDEBUG
181  red_count++;
182#ifdef TEST_OPT_DEBUG_RED
183  if (TEST_OPT_DEBUG)
184  {
185    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
186    PW->wrp();
187  }
188#endif
189#endif
190  int ret = 0;
191  ring tailRing = PR->tailRing;
192  assume(kTest_L(PR));
193  assume(kTest_T(PW));
194
195  // signature-based stuff:
196  // checking for sig-safeness first
197  // NOTE: This has to be done in the current ring
198  //
199  /**********************************************
200   *
201   * TODO:
202   * --------------------------------------------
203   * if strat->incremental
204   * Since we are subdividing lower index and
205   * current index reductions it is enough to
206   * look at the polynomial part of the signature
207   * for a check. This should speed-up checking
208   * a lot!
209   * if !strat->incremental
210   * We are not subdividing lower and current index
211   * due to the fact that we are using the induced
212   * Schreyer order
213   *
214   * nevertheless, this different behaviour is
215   * taken care of by is_sigsafe
216   * => one reduction procedure can be used for
217   * both, the incremental and the non-incremental
218   * attempt!
219   * --------------------------------------------
220   *
221   *********************************************/
222  //printf("COMPARE IDX: %ld -- %ld\n",idx,strat->currIdx);
223  if (!PW->is_sigsafe)
224  {
225    PR->SetLmCurrRing();
226    //poly f1 = pCopy(PR->GetLmTailRing());   // p2 | p1
227    //poly f2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
228    poly f1 = p_Copy(PR->GetLmCurrRing(),currRing);
229    poly f2 = PW->GetLmCurrRing();
230    poly sigMult = pCopy(PW->sig);   // copy signature of reducer
231    p_ExpVectorSub(f1, f2, currRing); // Calculate the Monomial we must multiply to p2
232//#if 1
233#ifdef DEBUGF5
234    printf("IN KSREDUCEPOLYSIG: \n");
235    pWrite(pHead(f1));
236    pWrite(pHead(f2));
237    pWrite(sigMult);
238    printf("--------------\n");
239#endif
240    sigMult = pp_Mult_qq(f1,sigMult,currRing);
241//#if 1
242#ifdef DEBUGF5
243    printf("------------------- IN KSREDUCEPOLYSIG: --------------------\n");
244    pWrite(pHead(f1));
245    pWrite(pHead(f2));
246    pWrite(sigMult);
247    pWrite(PR->sig);
248    printf("--------------\n");
249#endif
250    int sigSafe = p_LmCmp(PR->sig,sigMult,currRing);
251    // now we can delete the copied polynomial data used for checking for
252    // sig-safeness of the reduction step
253//#if 1
254#ifdef DEBUGF5
255    printf("%d -- %d sig\n",sigSafe,PW->is_sigsafe);
256
257#endif
258    pDelete(&f1);
259    pDelete(&sigMult);
260    // go on with the computations only if the signature of p2 is greater than the
261    // signature of fm*p1
262    if(sigSafe != 1)
263    {
264      PR->is_redundant = TRUE;
265      return 3;
266    }
267    //PW->is_sigsafe  = TRUE;
268  }
269  PR->is_redundant = FALSE;
270  poly p1 = PR->GetLmTailRing();   // p2 | p1
271  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
272  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
273  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
274  p_CheckPolyRing(p1, tailRing);
275  p_CheckPolyRing(p2, tailRing);
276
277  pAssume1(p2 != NULL && p1 != NULL &&
278      p_DivisibleBy(p2,  p1, tailRing));
279
280  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
281      (p_GetComp(p2, tailRing) == 0 &&
282       p_MaxComp(pNext(p2),tailRing) == 0));
283
284#ifdef HAVE_PLURAL
285  if (rIsPluralRing(currRing))
286  {
287    // for the time being: we know currRing==strat->tailRing
288    // no exp-bound checking needed
289    // (only needed if exp-bound(tailring)<exp-b(currRing))
290    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
291    else
292    {
293      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
294      assume(_p != NULL);
295      nc_PolyPolyRed(_p, p2, coef, currRing);
296      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
297      PR->pLength=0; // usaully not used, GetpLength re-comoutes it if needed
298    }
299    return 0;
300  }
301#endif
302
303  if (t2==NULL)           // Divisor is just one term, therefore it will
304  {                       // just cancel the leading term
305    PR->LmDeleteAndIter();
306    if (coef != NULL) *coef = n_Init(1, tailRing);
307    return 0;
308  }
309
310  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
311
312  if (tailRing != currRing)
313  {
314    // check that reduction does not violate exp bound
315    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
316    {
317      // undo changes of lm
318      p_ExpVectorAdd(lm, p2, tailRing);
319      if (strat == NULL) return 2;
320      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
321      tailRing = strat->tailRing;
322      p1 = PR->GetLmTailRing();
323      p2 = PW->GetLmTailRing();
324      t2 = pNext(p2);
325      lm = p1;
326      p_ExpVectorSub(lm, p2, tailRing);
327      ret = 1;
328    }
329  }
330
331  // take care of coef buisness
332  if (! n_IsOne(pGetCoeff(p2), tailRing))
333  {
334    number bn = pGetCoeff(lm);
335    number an = pGetCoeff(p2);
336    int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
337    p_SetCoeff(lm, bn, tailRing);
338    if ((ct == 0) || (ct == 2))
339      PR->Tail_Mult_nn(an);
340    if (coef != NULL) *coef = an;
341    else n_Delete(&an, tailRing);
342  }
343  else
344  {
345    if (coef != NULL) *coef = n_Init(1, tailRing);
346  }
347
348
349  // and finally,
350  PR->Tail_Minus_mm_Mult_qq(lm, t2, PW->GetpLength() - 1, spNoether);
351  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
352  PR->LmDeleteAndIter();
353
354  // the following is commented out: shrinking
355#ifdef HAVE_SHIFTBBA_NONEXISTENT
356  if ( (currRing->isLPring) && (!strat->homog) )
357  {
358    // assume? h->p in currRing
359    PR->GetP();
360    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
361    PR->Clear(); // does the right things
362    PR->p = qq;
363    PR->t_p = NULL;
364    PR->SetShortExpVector();
365  }
366#endif
367
368#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
369  if (TEST_OPT_DEBUG)
370  {
371    Print(" to: "); PR->wrp(); Print("\n");
372  }
373#endif
374  return ret;
375}
376
377/***************************************************************
378 *
379 * Creates S-Poly of p1 and p2
380 *
381 *
382 ***************************************************************/
383void ksCreateSpoly(LObject* Pair,   poly spNoether,
384                   int use_buckets, ring tailRing,
385                   poly m1, poly m2, TObject** R)
386{
387#ifdef KDEBUG
388  create_count++;
389#endif
390  assume(kTest_L(Pair));
391  poly p1 = Pair->p1;
392  poly p2 = Pair->p2;
393  Pair->tailRing = tailRing;
394
395  assume(p1 != NULL);
396  assume(p2 != NULL);
397  assume(tailRing != NULL);
398
399  poly a1 = pNext(p1), a2 = pNext(p2);
400  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
401  int co=0/*, ct = ksCheckCoeff(&lc1, &lc2, currRing->cf)*/; // gcd and zero divisors
402  (void) ksCheckCoeff(&lc1, &lc2, currRing->cf);
403
404  int l1=0, l2=0;
405
406  if (p_GetComp(p1, currRing)!=p_GetComp(p2, currRing))
407  {
408    if (p_GetComp(p1, currRing)==0)
409    {
410      co=1;
411      p_SetCompP(p1,p_GetComp(p2, currRing), currRing, tailRing);
412    }
413    else
414    {
415      co=2;
416      p_SetCompP(p2, p_GetComp(p1, currRing), currRing, tailRing);
417    }
418  }
419
420  // get m1 = LCM(LM(p1), LM(p2))/LM(p1)
421  //     m2 = LCM(LM(p1), LM(p2))/LM(p2)
422  if (m1 == NULL)
423    k_GetLeadTerms(p1, p2, currRing, m1, m2, tailRing);
424
425  pSetCoeff0(m1, lc2);
426  pSetCoeff0(m2, lc1);  // and now, m1 * LT(p1) == m2 * LT(p2)
427
428  if (R != NULL)
429  {
430    if (Pair->i_r1 == -1)
431    {
432      l1 = pLength(p1) - 1;
433    }
434    else
435    {
436      l1 = (R[Pair->i_r1])->GetpLength() - 1;
437    }
438    if (Pair->i_r2 == -1)
439    {
440      l2 = pLength(p2) - 1;
441    }
442    else
443    {
444      l2 = (R[Pair->i_r2])->GetpLength() - 1;
445    }
446  }
447
448  // get m2 * a2
449  if (spNoether != NULL)
450  {
451    l2 = -1;
452    a2 = tailRing->p_Procs->pp_Mult_mm_Noether(a2, m2, spNoether, l2, tailRing);
453    assume(l2 == pLength(a2));
454  }
455  else
456    a2 = tailRing->p_Procs->pp_Mult_mm(a2, m2, tailRing);
457#ifdef HAVE_RINGS
458  if (!(rField_is_Domain(currRing))) l2 = pLength(a2);
459#endif
460
461  Pair->SetLmTail(m2, a2, l2, use_buckets, tailRing);
462
463  // get m2*a2 - m1*a1
464  Pair->Tail_Minus_mm_Mult_qq(m1, a1, l1, spNoether);
465
466  // Clean-up time
467  Pair->LmDeleteAndIter();
468  p_LmDelete(m1, tailRing);
469
470  if (co != 0)
471  {
472    if (co==1)
473    {
474      p_SetCompP(p1,0, currRing, tailRing);
475    }
476    else
477    {
478      p_SetCompP(p2,0, currRing, tailRing);
479    }
480  }
481
482  // the following is commented out: shrinking
483#ifdef HAVE_SHIFTBBA_NONEXISTENT
484  if (currRing->isLPring)
485  {
486    // assume? h->p in currRing
487    Pair->GetP();
488    poly qq = p_Shrink(Pair->p, currRing->isLPring, currRing);
489    Pair->Clear(); // does the right things
490    Pair->p = qq;
491    Pair->t_p = NULL;
492    Pair->SetShortExpVector();
493  }
494#endif
495
496}
497
498int ksReducePolyTail(LObject* PR, TObject* PW, poly Current, poly spNoether)
499{
500  BOOLEAN ret;
501  number coef;
502  poly Lp =     PR->GetLmCurrRing();
503  poly Save =   PW->GetLmCurrRing();
504
505  assume(kTest_L(PR));
506  assume(kTest_T(PW));
507  pAssume(pIsMonomOf(Lp, Current));
508
509  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
510  assume(PR->bucket == NULL);
511
512  LObject Red(pNext(Current), PR->tailRing);
513  TObject With(PW, Lp == Save);
514
515  pAssume(!pHaveCommonMonoms(Red.p, With.p));
516  ret = ksReducePoly(&Red, &With, spNoether, &coef);
517
518  if (!ret)
519  {
520    if (! n_IsOne(coef, currRing))
521    {
522      pNext(Current) = NULL;
523      if (Current == PR->p && PR->t_p != NULL)
524        pNext(PR->t_p) = NULL;
525      PR->Mult_nn(coef);
526    }
527
528    n_Delete(&coef, currRing);
529    pNext(Current) = Red.GetLmTailRing();
530    if (Current == PR->p && PR->t_p != NULL)
531      pNext(PR->t_p) = pNext(Current);
532  }
533
534  if (Lp == Save)
535    With.Delete();
536
537  // the following is commented out: shrinking
538#ifdef HAVE_SHIFTBBA_NONEXISTENT
539  if (currRing->isLPring)
540  {
541    // assume? h->p in currRing
542    PR->GetP();
543    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
544    PR->Clear(); // does the right things
545    PR->p = qq;
546    PR->t_p = NULL;
547    PR->SetShortExpVector();
548  }
549#endif
550
551  return ret;
552}
553
554/***************************************************************
555 *
556 * Auxillary Routines
557 *
558 *
559 ***************************************************************/
560
561/*2
562* creates the leading term of the S-polynomial of p1 and p2
563* do not destroy p1 and p2
564* remarks:
565*   1. the coefficient is 0 (nNew)
566*   1. a) in the case of coefficient ring, the coefficient is calculated
567*   2. pNext is undefined
568*/
569//static void bbb() { int i=0; }
570poly ksCreateShortSpoly(poly p1, poly p2, ring tailRing)
571{
572  poly a1 = pNext(p1), a2 = pNext(p2);
573  long c1=p_GetComp(p1, currRing),c2=p_GetComp(p2, currRing);
574  long c;
575  poly m1,m2;
576  number t1 = NULL,t2 = NULL;
577  int cm,i;
578  BOOLEAN equal;
579
580#ifdef HAVE_RINGS
581  BOOLEAN is_Ring=rField_is_Ring(currRing);
582  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
583  if (is_Ring)
584  {
585    ksCheckCoeff(&lc1, &lc2, currRing->cf); // gcd and zero divisors
586    if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
587    if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
588    while (a1 != NULL && nIsZero(t2))
589    {
590      pIter(a1);
591      nDelete(&t2);
592      if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
593    }
594    while (a2 != NULL && nIsZero(t1))
595    {
596      pIter(a2);
597      nDelete(&t1);
598      if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
599    }
600  }
601#endif
602
603  if (a1==NULL)
604  {
605    if(a2!=NULL)
606    {
607      m2=p_Init(currRing);
608x2:
609      for (i = (currRing->N); i; i--)
610      {
611        c = p_GetExpDiff(p1, p2,i, currRing);
612        if (c>0)
613        {
614          p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)),currRing);
615        }
616        else
617        {
618          p_SetExp(m2,i,p_GetExp(a2,i,tailRing),currRing);
619        }
620      }
621      if ((c1==c2)||(c2!=0))
622      {
623        p_SetComp(m2,p_GetComp(a2,tailRing), currRing);
624      }
625      else
626      {
627        p_SetComp(m2,c1,currRing);
628      }
629      p_Setm(m2, currRing);
630#ifdef HAVE_RINGS
631      if (is_Ring)
632      {
633          nDelete(&lc1);
634          nDelete(&lc2);
635          nDelete(&t2);
636          pSetCoeff0(m2, t1);
637      }
638      else
639#endif
640        nNew(&(pGetCoeff(m2)));
641      return m2;
642    }
643    else
644    {
645#ifdef HAVE_RINGS
646      if (is_Ring)
647      {
648        nDelete(&lc1);
649        nDelete(&lc2);
650        nDelete(&t1);
651        nDelete(&t2);
652      }
653#endif
654      return NULL;
655    }
656  }
657  if (a2==NULL)
658  {
659    m1=p_Init(currRing);
660x1:
661    for (i = (currRing->N); i; i--)
662    {
663      c = p_GetExpDiff(p2, p1,i,currRing);
664      if (c>0)
665      {
666        p_SetExp(m1,i,(c+p_GetExp(a1,i, tailRing)),currRing);
667      }
668      else
669      {
670        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
671      }
672    }
673    if ((c1==c2)||(c1!=0))
674    {
675      p_SetComp(m1,p_GetComp(a1,tailRing),currRing);
676    }
677    else
678    {
679      p_SetComp(m1,c2,currRing);
680    }
681    p_Setm(m1, currRing);
682#ifdef HAVE_RINGS
683    if (is_Ring)
684    {
685      pSetCoeff0(m1, t2);
686      nDelete(&lc1);
687      nDelete(&lc2);
688      nDelete(&t1);
689    }
690    else
691#endif
692      nNew(&(pGetCoeff(m1)));
693    return m1;
694  }
695  m1 = p_Init(currRing);
696  m2 = p_Init(currRing);
697  loop
698  {
699    for (i = (currRing->N); i; i--)
700    {
701      c = p_GetExpDiff(p1, p2,i,currRing);
702      if (c > 0)
703      {
704        p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)), currRing);
705        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
706      }
707      else
708      {
709        p_SetExp(m1,i,(p_GetExp(a1,i,tailRing)-c), currRing);
710        p_SetExp(m2,i,p_GetExp(a2,i, tailRing), currRing);
711      }
712    }
713    if(c1==c2)
714    {
715      p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
716      p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
717    }
718    else
719    {
720      if(c1!=0)
721      {
722        p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
723        p_SetComp(m2,c1, currRing);
724      }
725      else
726      {
727        p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
728        p_SetComp(m1,c2, currRing);
729      }
730    }
731    p_Setm(m1,currRing);
732    p_Setm(m2,currRing);
733    cm = p_LmCmp(m1, m2,currRing);
734    if (cm!=0)
735    {
736      if(cm==1)
737      {
738        p_LmFree(m2,currRing);
739#ifdef HAVE_RINGS
740        if (is_Ring)
741        {
742          pSetCoeff0(m1, t2);
743          nDelete(&lc1);
744          nDelete(&lc2);
745          nDelete(&t1);
746        }
747        else
748#endif
749          nNew(&(pGetCoeff(m1)));
750        return m1;
751      }
752      else
753      {
754        p_LmFree(m1,currRing);
755#ifdef HAVE_RINGS
756        if (is_Ring)
757        {
758          pSetCoeff0(m2, t1);
759          nDelete(&lc1);
760          nDelete(&lc2);
761          nDelete(&t2);
762        }
763        else
764#endif
765          nNew(&(pGetCoeff(m2)));
766        return m2;
767      }
768    }
769#ifdef HAVE_RINGS
770    if (is_Ring)
771    {
772      equal = nEqual(t1,t2);
773    }
774    else
775#endif
776    {
777      t1 = nMult(pGetCoeff(a2),pGetCoeff(p1));
778      t2 = nMult(pGetCoeff(a1),pGetCoeff(p2));
779      equal = nEqual(t1,t2);
780      nDelete(&t2);
781      nDelete(&t1);
782    }
783    if (!equal)
784    {
785      p_LmFree(m2,currRing);
786#ifdef HAVE_RINGS
787      if (is_Ring)
788      {
789          pSetCoeff0(m1, nSub(t1, t2));
790          nDelete(&lc1);
791          nDelete(&lc2);
792          nDelete(&t1);
793          nDelete(&t2);
794      }
795      else
796#endif
797        nNew(&(pGetCoeff(m1)));
798      return m1;
799    }
800    pIter(a1);
801    pIter(a2);
802#ifdef HAVE_RINGS
803    if (is_Ring)
804    {
805      if (a2 != NULL)
806      {
807        nDelete(&t1);
808        t1 = nMult(pGetCoeff(a2),lc1);
809      }
810      if (a1 != NULL)
811      {
812        nDelete(&t2);
813        t2 = nMult(pGetCoeff(a1),lc2);
814      }
815      while ((a1 != NULL) && nIsZero(t2))
816      {
817        pIter(a1);
818        if (a1 != NULL)
819        {
820          nDelete(&t2);
821          t2 = nMult(pGetCoeff(a1),lc2);
822        }
823      }
824      while ((a2 != NULL) && nIsZero(t1))
825      {
826        pIter(a2);
827        if (a2 != NULL)
828        {
829          nDelete(&t1);
830          t1 = nMult(pGetCoeff(a2),lc1);
831        }
832      }
833    }
834#endif
835    if (a2==NULL)
836    {
837      p_LmFree(m2,currRing);
838      if (a1==NULL)
839      {
840#ifdef HAVE_RINGS
841        if (is_Ring)
842        {
843          nDelete(&lc1);
844          nDelete(&lc2);
845          nDelete(&t1);
846          nDelete(&t2);
847        }
848#endif
849        p_LmFree(m1,currRing);
850        return NULL;
851      }
852      goto x1;
853    }
854    if (a1==NULL)
855    {
856      p_LmFree(m1,currRing);
857      goto x2;
858    }
859  }
860}
Note: See TracBrowser for help on using the repository browser.