source: git/kernel/kspoly.cc @ 3901ebf

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