source: git/kernel/kspoly.cc @ 073e96

spielwiese
Last change on this file since 073e96 was 073e96, checked in by Burcin Erocal <burcin@…>, 13 years ago
Fix kernel/kspoly.cc.
  • Property mode set to 100644
File size: 14.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6*  ABSTRACT -  Routines for Spoly creation and reductions
7*/
8
9// #define PDEBUG 2
10#include <kernel/mod2.h>
11#include <misc/options.h>
12#include <kernel/kutil.h>
13#include <coeffs/numbers.h>
14#include <polys/monomials/p_polys.h>
15#include <polys/templates/p_Procs.h>
16#include <polys/nc/nc.h>
17#ifdef KDEBUG
18#include <kernel/febase.h>
19#endif
20#ifdef HAVE_RINGS
21#include <polys/polys.h>
22#endif
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 ***************************************************************/
37int ksReducePoly(LObject* PR,
38                 TObject* PW,
39                 poly spNoether,
40                 number *coef,
41                 kStrategy strat)
42{
43#ifdef KDEBUG
44  red_count++;
45#ifdef TEST_OPT_DEBUG_RED
46  if (TEST_OPT_DEBUG)
47  {
48    Print("Red %d:", red_count); PR->wrp(); Print(" with:");
49    PW->wrp();
50  }
51#endif
52#endif
53  int ret = 0;
54  ring tailRing = PR->tailRing;
55  kTest_L(PR);
56  kTest_T(PW);
57
58  poly p1 = PR->GetLmTailRing();   // p2 | p1
59  poly p2 = PW->GetLmTailRing();   // i.e. will reduce p1 with p2; lm = LT(p1) / LM(p2)
60  poly t2 = pNext(p2), lm = p1;    // t2 = p2 - LT(p2); really compute P = LC(p2)*p1 - LT(p1)/LM(p2)*p2
61  assume(p1 != NULL && p2 != NULL);// Attention, we have rings and there LC(p2) and LC(p1) are special
62  p_CheckPolyRing(p1, tailRing);
63  p_CheckPolyRing(p2, tailRing);
64
65  pAssume1(p2 != NULL && p1 != NULL &&
66           p_DivisibleBy(p2,  p1, tailRing));
67
68  pAssume1(p_GetComp(p1, tailRing) == p_GetComp(p2, tailRing) ||
69           (p_GetComp(p2, tailRing) == 0 &&
70            p_MaxComp(pNext(p2),tailRing) == 0));
71
72#ifdef HAVE_PLURAL
73  if (rIsPluralRing(currRing))
74  {
75    // for the time being: we know currRing==strat->tailRing
76    // no exp-bound checking needed
77    // (only needed if exp-bound(tailring)<exp-b(currRing))
78    if (PR->bucket!=NULL)  nc_kBucketPolyRed(PR->bucket, p2,coef);
79    else 
80    {
81      poly _p = (PR->t_p != NULL ? PR->t_p : PR->p);
82      assume(_p != NULL);
83      nc_PolyPolyRed(_p, p2,coef, currRing);
84      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
85      PR->pLength=0; // usually not used, GetpLength re-computes it if needed
86    }
87    return 0;
88  }
89#endif
90
91  if (t2==NULL)           // Divisor is just one term, therefore it will
92  {                       // just cancel the leading term
93    PR->LmDeleteAndIter();
94    if (coef != NULL) *coef = n_Init(1, tailRing);
95    return 0;
96  }
97
98  p_ExpVectorSub(lm, p2, tailRing); // Calculate the Monomial we must multiply to p2
99
100  if (tailRing != currRing)
101  {
102    // check that reduction does not violate exp bound
103    while (PW->max != NULL && !p_LmExpVectorAddIsOk(lm, PW->max, tailRing))
104    {
105      // undo changes of lm
106      p_ExpVectorAdd(lm, p2, tailRing);
107      if (strat == NULL) return 2;
108      if (! kStratChangeTailRing(strat, PR, PW)) return -1;
109      tailRing = strat->tailRing;
110      p1 = PR->GetLmTailRing();
111      p2 = PW->GetLmTailRing();
112      t2 = pNext(p2);
113      lm = p1;
114      p_ExpVectorSub(lm, p2, tailRing);
115      ret = 1;
116    }
117  }
118
119  // take care of coef buisness
120  if (! n_IsOne(pGetCoeff(p2), tailRing))
121  {
122    number bn = pGetCoeff(lm);
123    number an = pGetCoeff(p2);
124    int ct = ksCheckCoeff(&an, &bn, tailRing->cf);    // Calculate special LC
125    p_SetCoeff(lm, bn, tailRing);
126    if ((ct == 0) || (ct == 2))
127      PR->Tail_Mult_nn(an);
128    if (coef != NULL) *coef = an;
129    else n_Delete(&an, tailRing);
130  }
131  else
132  {
133    if (coef != NULL) *coef = n_Init(1, tailRing);
134  }
135
136
137  // and finally,
138  PR->Tail_Minus_mm_Mult_qq(lm, t2, PW->GetpLength() - 1, spNoether);
139  assume(PW->GetpLength() == pLength(PW->p != NULL ? PW->p : PW->t_p));
140  PR->LmDeleteAndIter();
141
142  // the following is commented out: shrinking
143#ifdef HAVE_SHIFTBBA_NONEXISTENT
144  if ( (currRing->isLPring) && (!strat->homog) )
145  {
146    // assume? h->p in currRing
147    PR->GetP();
148    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
149    PR->Clear(); // does the right things
150    PR->p = qq; 
151    PR->t_p = NULL;
152    PR->SetShortExpVector();
153  }
154#endif
155 
156#if defined(KDEBUG) && defined(TEST_OPT_DEBUG_RED)
157  if (TEST_OPT_DEBUG)
158  {
159    Print(" to: "); PR->wrp(); Print("\n");
160  }
161#endif
162  return ret;
163}
164
165/***************************************************************
166 *
167 * Creates S-Poly of p1 and p2
168 *
169 *
170 ***************************************************************/
171void ksCreateSpoly(LObject* Pair,   poly spNoether,
172                   int use_buckets, ring tailRing,
173                   poly m1, poly m2, TObject** R)
174{
175#ifdef KDEBUG
176  create_count++;
177#endif
178  kTest_L(Pair);
179  poly p1 = Pair->p1;
180  poly p2 = Pair->p2;
181  poly last;
182  Pair->tailRing = tailRing;
183
184  assume(p1 != NULL);
185  assume(p2 != NULL);
186  assume(tailRing != NULL);
187
188  poly a1 = pNext(p1), a2 = pNext(p2);
189  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
190  int co=0, ct = ksCheckCoeff(&lc1, &lc2, currRing->cf); // gcd and zero divisors
191
192  int l1=0, l2=0;
193
194  if (p_GetComp(p1, currRing)!=p_GetComp(p2, currRing))
195  {
196    if (p_GetComp(p1, currRing)==0)
197    {
198      co=1;
199      p_SetCompP(p1,p_GetComp(p2, currRing), currRing, tailRing);
200    }
201    else
202    {
203      co=2;
204      p_SetCompP(p2, p_GetComp(p1, currRing), currRing, tailRing);
205    }
206  }
207
208  // get m1 = LCM(LM(p1), LM(p2))/LM(p1)
209  //     m2 = LCM(LM(p1), LM(p2))/LM(p2)
210  if (m1 == NULL)
211    k_GetLeadTerms(p1, p2, currRing, m1, m2, tailRing);
212
213  pSetCoeff0(m1, lc2);
214  pSetCoeff0(m2, lc1);  // and now, m1 * LT(p1) == m2 * LT(p2)
215
216  if (R != NULL)
217  {
218    if (Pair->i_r1 == -1)
219    {
220      l1 = pLength(p1) - 1;
221    }
222    else
223    {
224      l1 = (R[Pair->i_r1])->GetpLength() - 1;
225    }
226    if (Pair->i_r2 == -1)
227    {
228      l2 = pLength(p2) - 1;
229    }
230    else
231    {
232      l2 = (R[Pair->i_r2])->GetpLength() - 1;
233    }
234  }
235
236  // get m2 * a2
237  if (spNoether != NULL)
238  {
239    l2 = -1;
240    a2 = tailRing->p_Procs->pp_Mult_mm_Noether(a2, m2, spNoether, l2, tailRing,last);
241    assume(l2 == pLength(a2));
242  }
243  else
244    a2 = tailRing->p_Procs->pp_Mult_mm(a2, m2, tailRing,last);
245#ifdef HAVE_RINGS
246  if (!(rField_is_Domain(currRing))) l2 = pLength(a2);
247#endif
248
249  Pair->SetLmTail(m2, a2, l2, use_buckets, tailRing, last);
250
251  // get m2*a2 - m1*a1
252  Pair->Tail_Minus_mm_Mult_qq(m1, a1, l1, spNoether);
253
254  // Clean-up time
255  Pair->LmDeleteAndIter();
256  p_LmDelete(m1, tailRing);
257
258  if (co != 0)
259  {
260    if (co==1)
261    {
262      p_SetCompP(p1,0, currRing, tailRing);
263    }
264    else
265    {
266      p_SetCompP(p2,0, currRing, tailRing);
267    }
268  }
269
270  // the following is commented out: shrinking
271#ifdef HAVE_SHIFTBBA_NONEXISTENT
272  if (currRing->isLPring)
273  {
274    // assume? h->p in currRing
275    Pair->GetP();
276    poly qq = p_Shrink(Pair->p, currRing->isLPring, currRing);
277    Pair->Clear(); // does the right things
278    Pair->p = qq; 
279    Pair->t_p = NULL;
280    Pair->SetShortExpVector();
281  }
282#endif
283
284}
285
286int ksReducePolyTail(LObject* PR, TObject* PW, poly Current, poly spNoether)
287{
288  BOOLEAN ret;
289  number coef;
290  poly Lp =     PR->GetLmCurrRing();
291  poly Save =   PW->GetLmCurrRing();
292
293  kTest_L(PR);
294  kTest_T(PW);
295  pAssume(pIsMonomOf(Lp, Current));
296
297  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
298  assume(PR->bucket == NULL);
299
300  LObject Red(pNext(Current), PR->tailRing);
301  TObject With(PW, Lp == Save);
302
303  pAssume(!pHaveCommonMonoms(Red.p, With.p));
304  ret = ksReducePoly(&Red, &With, spNoether, &coef);
305
306  if (!ret)
307  {
308    if (! n_IsOne(coef, currRing))
309    {
310      pNext(Current) = NULL;
311      if (Current == PR->p && PR->t_p != NULL)
312        pNext(PR->t_p) = NULL;
313      PR->Mult_nn(coef);
314    }
315
316    n_Delete(&coef, currRing);
317    pNext(Current) = Red.GetLmTailRing();
318    if (Current == PR->p && PR->t_p != NULL)
319      pNext(PR->t_p) = pNext(Current);
320  }
321
322  if (Lp == Save)
323    With.Delete();
324
325  // the following is commented out: shrinking
326#ifdef HAVE_SHIFTBBA_NONEXISTENT
327  if (currRing->isLPring)
328  {
329    // assume? h->p in currRing
330    PR->GetP();
331    poly qq = p_Shrink(PR->p, currRing->isLPring, currRing);
332    PR->Clear(); // does the right things
333    PR->p = qq; 
334    PR->t_p = NULL;
335    PR->SetShortExpVector();
336  }
337#endif
338
339  return ret;
340}
341
342/***************************************************************
343 *
344 * Auxillary Routines
345 *
346 *
347 ***************************************************************/
348
349/*2
350* creates the leading term of the S-polynomial of p1 and p2
351* do not destroy p1 and p2
352* remarks:
353*   1. the coefficient is 0 (nNew)
354*   1. a) in the case of coefficient ring, the coefficient is calculated
355*   2. pNext is undefined
356*/
357//static void bbb() { int i=0; }
358poly ksCreateShortSpoly(poly p1, poly p2, ring tailRing)
359{
360  poly a1 = pNext(p1), a2 = pNext(p2);
361  long c1=p_GetComp(p1, currRing),c2=p_GetComp(p2, currRing);
362  long c;
363  poly m1,m2;
364  number t1 = NULL,t2 = NULL;
365  int cm,i;
366  BOOLEAN equal;
367
368#ifdef HAVE_RINGS
369  BOOLEAN is_Ring=rField_is_Ring(currRing);
370  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
371  if (is_Ring)
372  {
373    ksCheckCoeff(&lc1, &lc2, currRing->cf); // gcd and zero divisors
374    if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
375    if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
376    while (a1 != NULL && nIsZero(t2))
377    {
378      pIter(a1);
379      nDelete(&t2);
380      if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
381    }
382    while (a2 != NULL && nIsZero(t1))
383    {
384      pIter(a2);
385      nDelete(&t1);
386      if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
387    }
388  }
389#endif
390
391  if (a1==NULL)
392  {
393    if(a2!=NULL)
394    {
395      m2=p_Init(currRing);
396x2:
397      for (i = (currRing->N); i; i--)
398      {
399        c = p_GetExpDiff(p1, p2,i, currRing);
400        if (c>0)
401        {
402          p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)),currRing);
403        }
404        else
405        {
406          p_SetExp(m2,i,p_GetExp(a2,i,tailRing),currRing);
407        }
408      }
409      if ((c1==c2)||(c2!=0))
410      {
411        p_SetComp(m2,p_GetComp(a2,tailRing), currRing);
412      }
413      else
414      {
415        p_SetComp(m2,c1,currRing);
416      }
417      p_Setm(m2, currRing);
418#ifdef HAVE_RINGS
419      if (is_Ring)
420      {
421          nDelete(&lc1);
422          nDelete(&lc2);
423          nDelete(&t2);
424          pSetCoeff0(m2, t1);
425      }
426      else
427#endif
428        nNew(&(pGetCoeff(m2)));
429      return m2;
430    }
431    else
432    {
433#ifdef HAVE_RINGS
434      if (is_Ring)
435      {
436        nDelete(&lc1);
437        nDelete(&lc2);
438        nDelete(&t1);
439        nDelete(&t2);
440      }
441#endif
442      return NULL;
443    }
444  }
445  if (a2==NULL)
446  {
447    m1=p_Init(currRing);
448x1:
449    for (i = (currRing->N); i; i--)
450    {
451      c = p_GetExpDiff(p2, p1,i,currRing);
452      if (c>0)
453      {
454        p_SetExp(m1,i,(c+p_GetExp(a1,i, tailRing)),currRing);
455      }
456      else
457      {
458        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
459      }
460    }
461    if ((c1==c2)||(c1!=0))
462    {
463      p_SetComp(m1,p_GetComp(a1,tailRing),currRing);
464    }
465    else
466    {
467      p_SetComp(m1,c2,currRing);
468    }
469    p_Setm(m1, currRing);
470#ifdef HAVE_RINGS
471    if (is_Ring)
472    {
473      pSetCoeff0(m1, t2);
474      nDelete(&lc1);
475      nDelete(&lc2);
476      nDelete(&t1);
477    }
478    else
479#endif
480      nNew(&(pGetCoeff(m1)));
481    return m1;
482  }
483  m1 = p_Init(currRing);
484  m2 = p_Init(currRing);
485  loop
486  {
487    for (i = (currRing->N); i; i--)
488    {
489      c = p_GetExpDiff(p1, p2,i,currRing);
490      if (c > 0)
491      {
492        p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)), currRing);
493        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
494      }
495      else
496      {
497        p_SetExp(m1,i,(p_GetExp(a1,i,tailRing)-c), currRing);
498        p_SetExp(m2,i,p_GetExp(a2,i, tailRing), currRing);
499      }
500    }
501    if(c1==c2)
502    {
503      p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
504      p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
505    }
506    else
507    {
508      if(c1!=0)
509      {
510        p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
511        p_SetComp(m2,c1, currRing);
512      }
513      else
514      {
515        p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
516        p_SetComp(m1,c2, currRing);
517      }
518    }
519    p_Setm(m1,currRing);
520    p_Setm(m2,currRing);
521    cm = p_LmCmp(m1, m2,currRing);
522    if (cm!=0)
523    {
524      if(cm==1)
525      {
526        p_LmFree(m2,currRing);
527#ifdef HAVE_RINGS
528        if (is_Ring)
529        {
530          pSetCoeff0(m1, t2);
531          nDelete(&lc1);
532          nDelete(&lc2);
533          nDelete(&t1);
534        }
535        else
536#endif
537          nNew(&(pGetCoeff(m1)));
538        return m1;
539      }
540      else
541      {
542        p_LmFree(m1,currRing);
543#ifdef HAVE_RINGS
544        if (is_Ring)
545        {
546          pSetCoeff0(m2, t1);
547          nDelete(&lc1);
548          nDelete(&lc2);
549          nDelete(&t2);
550        }
551        else
552#endif
553          nNew(&(pGetCoeff(m2)));
554        return m2;
555      }
556    }
557#ifdef HAVE_RINGS
558    if (is_Ring)
559    {
560      equal = nEqual(t1,t2);
561    }
562    else
563#endif
564    {
565      t1 = nMult(pGetCoeff(a2),pGetCoeff(p1));
566      t2 = nMult(pGetCoeff(a1),pGetCoeff(p2));
567      equal = nEqual(t1,t2);
568      nDelete(&t2);
569      nDelete(&t1);
570    }
571    if (!equal)
572    {
573      p_LmFree(m2,currRing);
574#ifdef HAVE_RINGS
575      if (is_Ring)
576      {
577          pSetCoeff0(m1, nSub(t1, t2));
578          nDelete(&lc1);
579          nDelete(&lc2);
580          nDelete(&t1);
581          nDelete(&t2);
582      }
583      else
584#endif
585        nNew(&(pGetCoeff(m1)));
586      return m1;
587    }
588    pIter(a1);
589    pIter(a2);
590#ifdef HAVE_RINGS
591    if (is_Ring)
592    {
593      if (a2 != NULL)
594      {
595        nDelete(&t1);
596        t1 = nMult(pGetCoeff(a2),lc1);
597      }
598      if (a1 != NULL)
599      {
600        nDelete(&t2);
601        t2 = nMult(pGetCoeff(a1),lc2);
602      }
603      while ((a1 != NULL) && nIsZero(t2))
604      {
605        pIter(a1);
606        if (a1 != NULL)
607        {
608          nDelete(&t2);
609          t2 = nMult(pGetCoeff(a1),lc2);
610        }
611      }
612      while ((a2 != NULL) && nIsZero(t1))
613      {
614        pIter(a2);
615        if (a2 != NULL)
616        {
617          nDelete(&t1);
618          t1 = nMult(pGetCoeff(a2),lc1);
619        }
620      }
621    }
622#endif
623    if (a2==NULL)
624    {
625      p_LmFree(m2,currRing);
626      if (a1==NULL)
627      {
628#ifdef HAVE_RINGS
629        if (is_Ring)
630        {
631          nDelete(&lc1);
632          nDelete(&lc2);
633          nDelete(&t1);
634          nDelete(&t2);
635        }
636#endif
637        p_LmFree(m1,currRing);
638        return NULL;
639      }
640      goto x1;
641    }
642    if (a1==NULL)
643    {
644      p_LmFree(m1,currRing);
645      goto x2;
646    }
647  }
648}
Note: See TracBrowser for help on using the repository browser.