source: git/kernel/kspoly.cc @ 690e21e

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