source: git/kernel/kspoly.cc @ 0758b5

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