source: git/kernel/kspoly.cc @ 3c473c

spielwiese
Last change on this file since 3c473c was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 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 <kernel/mod2.h>
11#include <kernel/options.h>
12#include <kernel/kutil.h>
13#include <kernel/numbers.h>
14#include <kernel/p_polys.h>
15#include <kernel/p_Procs.h>
16#include <kernel/gring.h>
17#ifdef KDEBUG
18#include <kernel/febase.h>
19#endif
20#ifdef HAVE_RINGS
21#include <kernel/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);
84      if (PR->t_p!=NULL) PR->t_p=_p; else PR->p=_p;
85      PR->pLength=0; // usaully not used, GetpLength re-comoutes 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);    // 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); // 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/*
350* input - output: a, b
351* returns:
352*   a := a/gcd(a,b), b := b/gcd(a,b)
353*   and return value
354*       0  ->  a != 1,  b != 1
355*       1  ->  a == 1,  b != 1
356*       2  ->  a != 1,  b == 1
357*       3  ->  a == 1,  b == 1
358*   this value is used to control the spolys
359*/
360int ksCheckCoeff(number *a, number *b)
361{
362  int c = 0;
363  number an = *a, bn = *b;
364  nTest(an);
365  nTest(bn);
366
367  number cn = nGcd(an, bn, currRing);
368
369  if(nIsOne(cn))
370  {
371    an = nCopy(an);
372    bn = nCopy(bn);
373  }
374  else
375  {
376    an = nIntDiv(an, cn);
377    bn = nIntDiv(bn, cn);
378  }
379  nDelete(&cn);
380  if (nIsOne(an))
381  {
382    c = 1;
383  }
384  if (nIsOne(bn))
385  {
386    c += 2;
387  }
388  *a = an;
389  *b = bn;
390  return c;
391}
392
393/*2
394* creates the leading term of the S-polynomial of p1 and p2
395* do not destroy p1 and p2
396* remarks:
397*   1. the coefficient is 0 (nNew)
398*   1. a) in the case of coefficient ring, the coefficient is calculated
399*   2. pNext is undefined
400*/
401//static void bbb() { int i=0; }
402poly ksCreateShortSpoly(poly p1, poly p2, ring tailRing)
403{
404  poly a1 = pNext(p1), a2 = pNext(p2);
405  long c1=p_GetComp(p1, currRing),c2=p_GetComp(p2, currRing);
406  long c;
407  poly m1,m2;
408  number t1 = NULL,t2 = NULL;
409  int cm,i;
410  BOOLEAN equal;
411
412#ifdef HAVE_RINGS
413  BOOLEAN is_Ring=rField_is_Ring(currRing);
414  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
415  if (is_Ring)
416  {
417    ksCheckCoeff(&lc1, &lc2); // gcd and zero divisors
418    if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
419    if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
420    while (a1 != NULL && nIsZero(t2))
421    {
422      pIter(a1);
423      nDelete(&t2);
424      if (a1 != NULL) t2 = nMult(pGetCoeff(a1),lc2);
425    }
426    while (a2 != NULL && nIsZero(t1))
427    {
428      pIter(a2);
429      nDelete(&t1);
430      if (a2 != NULL) t1 = nMult(pGetCoeff(a2),lc1);
431    }
432  }
433#endif
434
435  if (a1==NULL)
436  {
437    if(a2!=NULL)
438    {
439      m2=p_Init(currRing);
440x2:
441      for (i = pVariables; i; i--)
442      {
443        c = p_GetExpDiff(p1, p2,i, currRing);
444        if (c>0)
445        {
446          p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)),currRing);
447        }
448        else
449        {
450          p_SetExp(m2,i,p_GetExp(a2,i,tailRing),currRing);
451        }
452      }
453      if ((c1==c2)||(c2!=0))
454      {
455        p_SetComp(m2,p_GetComp(a2,tailRing), currRing);
456      }
457      else
458      {
459        p_SetComp(m2,c1,currRing);
460      }
461      p_Setm(m2, currRing);
462#ifdef HAVE_RINGS
463      if (is_Ring)
464      {
465          nDelete(&lc1);
466          nDelete(&lc2);
467          nDelete(&t2);
468          pSetCoeff0(m2, t1);
469      }
470      else
471#endif
472        nNew(&(pGetCoeff(m2)));
473      return m2;
474    }
475    else
476    {
477#ifdef HAVE_RINGS
478      if (is_Ring)
479      {
480        nDelete(&lc1);
481        nDelete(&lc2);
482        nDelete(&t1);
483        nDelete(&t2);
484      }
485#endif
486      return NULL;
487    }
488  }
489  if (a2==NULL)
490  {
491    m1=p_Init(currRing);
492x1:
493    for (i = pVariables; i; i--)
494    {
495      c = p_GetExpDiff(p2, p1,i,currRing);
496      if (c>0)
497      {
498        p_SetExp(m1,i,(c+p_GetExp(a1,i, tailRing)),currRing);
499      }
500      else
501      {
502        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
503      }
504    }
505    if ((c1==c2)||(c1!=0))
506    {
507      p_SetComp(m1,p_GetComp(a1,tailRing),currRing);
508    }
509    else
510    {
511      p_SetComp(m1,c2,currRing);
512    }
513    p_Setm(m1, currRing);
514#ifdef HAVE_RINGS
515    if (is_Ring)
516    {
517      pSetCoeff0(m1, t2);
518      nDelete(&lc1);
519      nDelete(&lc2);
520      nDelete(&t1);
521    }
522    else
523#endif
524      nNew(&(pGetCoeff(m1)));
525    return m1;
526  }
527  m1 = p_Init(currRing);
528  m2 = p_Init(currRing);
529  loop
530  {
531    for (i = pVariables; i; i--)
532    {
533      c = p_GetExpDiff(p1, p2,i,currRing);
534      if (c > 0)
535      {
536        p_SetExp(m2,i,(c+p_GetExp(a2,i,tailRing)), currRing);
537        p_SetExp(m1,i,p_GetExp(a1,i, tailRing), currRing);
538      }
539      else
540      {
541        p_SetExp(m1,i,(p_GetExp(a1,i,tailRing)-c), currRing);
542        p_SetExp(m2,i,p_GetExp(a2,i, tailRing), currRing);
543      }
544    }
545    if(c1==c2)
546    {
547      p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
548      p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
549    }
550    else
551    {
552      if(c1!=0)
553      {
554        p_SetComp(m1,p_GetComp(a1, tailRing), currRing);
555        p_SetComp(m2,c1, currRing);
556      }
557      else
558      {
559        p_SetComp(m2,p_GetComp(a2, tailRing), currRing);
560        p_SetComp(m1,c2, currRing);
561      }
562    }
563    p_Setm(m1,currRing);
564    p_Setm(m2,currRing);
565    cm = p_LmCmp(m1, m2,currRing);
566    if (cm!=0)
567    {
568      if(cm==1)
569      {
570        p_LmFree(m2,currRing);
571#ifdef HAVE_RINGS
572        if (is_Ring)
573        {
574          pSetCoeff0(m1, t2);
575          nDelete(&lc1);
576          nDelete(&lc2);
577          nDelete(&t1);
578        }
579        else
580#endif
581          nNew(&(pGetCoeff(m1)));
582        return m1;
583      }
584      else
585      {
586        p_LmFree(m1,currRing);
587#ifdef HAVE_RINGS
588        if (is_Ring)
589        {
590          pSetCoeff0(m2, t1);
591          nDelete(&lc1);
592          nDelete(&lc2);
593          nDelete(&t2);
594        }
595        else
596#endif
597          nNew(&(pGetCoeff(m2)));
598        return m2;
599      }
600    }
601#ifdef HAVE_RINGS
602    if (is_Ring)
603    {
604      equal = nEqual(t1,t2);
605    }
606    else
607#endif
608    {
609      t1 = nMult(pGetCoeff(a2),pGetCoeff(p1));
610      t2 = nMult(pGetCoeff(a1),pGetCoeff(p2));
611      equal = nEqual(t1,t2);
612      nDelete(&t2);
613      nDelete(&t1);
614    }
615    if (!equal)
616    {
617      p_LmFree(m2,currRing);
618#ifdef HAVE_RINGS
619      if (is_Ring)
620      {
621          pSetCoeff0(m1, nSub(t1, t2));
622          nDelete(&lc1);
623          nDelete(&lc2);
624          nDelete(&t1);
625          nDelete(&t2);
626      }
627      else
628#endif
629        nNew(&(pGetCoeff(m1)));
630      return m1;
631    }
632    pIter(a1);
633    pIter(a2);
634#ifdef HAVE_RINGS
635    if (is_Ring)
636    {
637      if (a2 != NULL)
638      {
639        nDelete(&t1);
640        t1 = nMult(pGetCoeff(a2),lc1);
641      }
642      if (a1 != NULL)
643      {
644        nDelete(&t2);
645        t2 = nMult(pGetCoeff(a1),lc2);
646      }
647      while ((a1 != NULL) && nIsZero(t2))
648      {
649        pIter(a1);
650        if (a1 != NULL)
651        {
652          nDelete(&t2);
653          t2 = nMult(pGetCoeff(a1),lc2);
654        }
655      }
656      while ((a2 != NULL) && nIsZero(t1))
657      {
658        pIter(a2);
659        if (a2 != NULL)
660        {
661          nDelete(&t1);
662          t1 = nMult(pGetCoeff(a2),lc1);
663        }
664      }
665    }
666#endif
667    if (a2==NULL)
668    {
669      p_LmFree(m2,currRing);
670      if (a1==NULL)
671      {
672#ifdef HAVE_RINGS
673        if (is_Ring)
674        {
675          nDelete(&lc1);
676          nDelete(&lc2);
677          nDelete(&t1);
678          nDelete(&t2);
679        }
680#endif
681        p_LmFree(m1,currRing);
682        return NULL;
683      }
684      goto x1;
685    }
686    if (a1==NULL)
687    {
688      p_LmFree(m1,currRing);
689      goto x2;
690    }
691  }
692}
Note: See TracBrowser for help on using the repository browser.