source: git/kernel/kspoly.cc @ cea6f3

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