source: git/Singular/kbPolyProcs.cc @ c232af

spielwiese
Last change on this file since c232af was c232af, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* omalloc stuff git-svn-id: file:///usr/local/Singular/svn/trunk@4524 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.0 KB
Line 
1#include "mod2.h"
2#include <omalloc.h>
3#include "polys.h"
4#include "febase.h"
5#include "kutil.h"
6#include "kbuckets.h"
7#include "ring.h"
8#include "ipid.h"
9#include "numbers.h"
10#include "kbPolyProcs.h"
11#include "modulop.h"
12#include "spolys0.h"
13#include "polys-comp.h"
14
15#define HAVE_COMP_MACROS
16
17/***************************************************************
18 *
19 * Reduces PR with PW
20 * Assumes PR != NULL, PW != NULL, Lm(PR) divides Lm(PW)
21 *
22 ***************************************************************/
23void kbReducePoly(LObject* PR,
24                  TObject* PW,
25                  poly spNoether)
26{
27  assume(kTest_L(PR));
28  assume(kTest_T(PW));
29
30  poly p1 = PR->p;
31  poly p2 = PW->p;
32  assume(p2 != NULL && p1 != NULL && pDivisibleBy(p2,  p1));
33
34  poly a2 = pNext(p2), lm = p1;
35  p1 = pNext(p1);
36  PR->p = p1;
37  BOOLEAN reset_vec=FALSE;
38 
39  if (a2==NULL)
40  {
41    pDelete1(&lm);
42    return;
43  }
44
45  if (! nIsOne(pGetCoeff(p2)))
46  {
47    number an = pGetCoeff(p2), bn = pGetCoeff(lm);
48    int ct = spCheckCoeff(&an, &bn);
49    pSetCoeff(lm, bn);
50    if ((ct == 0) || (ct == 2)) kb_n_Mult_p(an, PR->p);
51    nDelete(&an);
52  }
53 
54  if (pGetComp(p2) != pGetComp(lm))
55  {
56    pSetCompP(p2, pGetComp(lm));
57    reset_vec = TRUE;
58  }
59
60  pMonSubFast(lm, p2);
61  assume(pGetComp(lm) == 0);
62  lm->next = NULL;
63  pTest(lm);
64  int l1, l2;
65  kb_p_Minus_m_Mult_q(&(PR->p), &l1, 
66                      lm, 
67                      a2, l2, spNoether);
68 
69  pDelete1(&lm);
70  if (reset_vec) spModuleToPoly(p2);
71}
72
73/***************************************************************
74 *
75 * Creates S-Poly of p1 and p2
76 *
77 *
78 ***************************************************************/
79void kbCreateSpoly(LObject* Pair, 
80                   poly spNoether)
81{
82  assume(kTest_L(Pair));
83  poly p1 = Pair->p1;
84  poly p2 = Pair->p2;
85
86  assume(p1 != NULL);
87  assume(p2 != NULL);
88
89  poly a1 = pNext(p1), a2 = pNext(p2);
90  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
91  poly m1, m2;
92  int co=0, ct = spCheckCoeff(&lc1, &lc2);
93  Exponent_t x;
94  int l1;
95
96  if (pGetComp(p1)!=pGetComp(p2))
97  {
98    if (pGetComp(p1)==0)
99    {
100      co=1;
101      pSetCompP(p1,pGetComp(p2));
102    }
103    else
104    {
105      co=2;
106      pSetCompP(p2,pGetComp(p1));
107    }
108  }
109
110  // get m1 = LCM(LM(p1), LM(p2))/LM(p1)
111  //     m2 = LCM(LM(p1), LM(p2))/LM(p2)
112  m1 = pInit();
113  m2 = pInit();
114  for (int i = pVariables; i; i--)
115  {
116    x = pGetExpDiff(p1, p2, i);
117    if (x > 0)
118    {
119      pSetExp(m2,i,x);
120      pSetExp(m1,i,0);
121    }
122    else
123    {
124      pSetExp(m1,i,-x);
125      pSetExp(m2,i,0);
126    }
127  }
128  pSetm(m1);
129  pSetm(m2);           // now we have m1 * LM(p1) == m2 * LM(p2)
130  pSetCoeff0(m1, lc2);
131  pSetCoeff0(m2, lc1); // and now, m1 * LT(p1) == m2 * LT(p2)
132
133  // get m2 * a2
134  Pair->p = kb_p_Mult_m(a2, m2);
135
136  int l2;
137  // and, finally, the spoly
138  kb_p_Minus_m_Mult_q(&(Pair->p), &l2, m1, a1, l1);
139 
140  // Clean-up time
141  pDelete1(&m1);
142  pDelete1(&m2);
143 
144  if (co != 0)
145  {
146    if (co==1)
147    {
148      spModuleToPoly(p1);
149    }
150    else
151    {
152      spModuleToPoly(p2);
153    }
154  }
155}
156
157/***************************************************************
158 *
159 * General kb_n_Mult_p which always works
160 * Realizes p = n*p
161 *
162 ***************************************************************/
163
164void kb_n_Mult_p_General(number n, poly p)
165{
166  while (p != NULL)
167  {
168    number nc = pGetCoeff(p);
169    pSetCoeff0(p, nMult(n, nc));
170    nDelete(&nc);
171    pIter(p);
172  }
173}
174
175/***************************************************************
176 *
177 * General kb_p_Add_q which always works
178 *
179 * assume pLength(*p) == *lp && pLength(*q) == *q
180 * *p and *q from heap
181 * Destroys *p and *q
182 * On return, *p == *p + *q, *q == NULL,
183 *            *lp == pLength(*p + *q), *lq == 0
184 *
185 ***************************************************************/
186
187void kb_p_Add_q_General(poly *p, int *lp,
188                        poly *q, int *lq, 
189                        omBin heap)
190{
191#ifdef KB_USE_HEAPS
192  assume(heap != NULL);
193#else
194  assume(heap == NULL);
195#endif
196#ifdef HAVE_LENGTH
197  assume(pLength(*p) == *lp && pLength(*q) == *lq);
198#endif
199 
200  number t, n1, n2;
201  unsigned int l = *lp + *lq;
202  spolyrec rp;
203  poly a = &rp, a1 = *p, a2 = *q;
204
205  if (a2 == NULL) return;
206
207  *q = NULL;
208  *lq = 0;
209 
210  if (a1 == NULL) 
211  {
212    *p = a2;
213    *lp = l;
214    return;
215  }
216
217  Top:     // compare a1 and a2 w.r.t. monomial ordering
218  register long d;
219#ifdef HAVE_COMP_MACROS
220  _pMonComp(a1, a2, d, goto NotEqual , goto Equal);
221#else
222  if ((d = pComp0(a1, a2))) goto NotEqual; else goto Equal;
223#endif
224
225  Equal:
226  assume(pComp0(a1, a2) == 0);
227
228  n1 = pGetCoeff(a1);
229  n2 = pGetCoeff(a2);
230  t = nAdd(n1,n2);
231  nDelete(&n1);
232  nDelete(&n2);
233  kb_pFree1AndAdvance(a2, heap);
234 
235  if (nIsZero(t))
236  {
237    l -= 2;
238    nDelete(&t);
239    kb_pFree1AndAdvance(a1, heap);
240  }
241  else
242  {
243    l--;
244    pSetCoeff0(a1,t);
245    a = pNext(a) = a1;
246    pIter(a1);
247  }
248  if (a1==NULL)
249  {
250    pNext(a) = a2;
251    goto Finish;
252  }
253  else if (a2==NULL)
254  {
255    pNext(a) = a1;
256    goto Finish;
257  }
258  goto Top;
259     
260  NotEqual:
261  if (d < 0)
262  {
263    assume(pComp0(a1, a2) == -1);
264    a = pNext(a) = a2;
265    pIter(a2);
266    if (a2==NULL)
267    {
268      pNext(a) = a1;
269      goto Finish;
270    }
271  }
272  else 
273  {
274    assume(pComp0(a1, a2) == 1);
275    a = pNext(a) = a1;
276    pIter(a1);
277    if (a1==NULL)
278    {
279      pNext(a) = a2;
280      goto Finish;
281    }
282  }
283  goto Top;
284 
285
286  Finish:
287#ifdef HAVE_LENGTH
288  assume(pLength(pNext(&rp)) == (int) l);
289#endif
290  *lp = l;
291  *p  = pNext(&rp);
292}
293
294/***************************************************************
295 *
296 * General kb_p_Mult_m which always works
297 *
298 * assume pLength(m) == 1
299 * Returns m*a1, monoms are allocated from heap
300 *
301 ***************************************************************/
302 
303poly  kb_p_Mult_m_General(poly p,
304                          poly m, 
305                          poly spNoether,
306                          omBin heap)
307{
308#ifdef KB_USE_HEAPS
309  assume(heap != NULL);
310#else
311  assume(heap == NULL);
312#endif
313  assume(pGetComp(m) == 0);
314
315  if (p == NULL) return NULL;
316  spolyrec rp;
317  poly q = &rp;
318  number ln = pGetCoeff(m);
319
320  while (p != NULL)
321  {
322    kb_pNew(pNext(q), heap);
323    q = pNext(q);
324
325    pSetCoeff0(q, nMult(ln, pGetCoeff(p)));
326
327    memaddW((unsigned long*) &(q->exp.l[0]),
328            (unsigned long*) &(p->exp.l[0]),
329            (unsigned long*) &(m->exp.l[0]),
330            currRing->ExpLSize);
331   
332    p = pNext(p);
333  }
334  pNext(q) = NULL;
335  return rp.next;
336}
337
338/***************************************************************
339 *
340 * General spoly loop which always works
341 *
342 * assume(pLength(*pp) == *lpp)
343 * assume(pLength(q) == lq)
344 * assume(pLength(m) == NULL)
345 * assume: Monoms of *pp are from heap
346 *
347 * Destroys *pp
348 * Does not touch q
349 * On return *pp == *pp - m*q
350 *           *lpp == pLength(*pp - m*q)
351 * New monoms are allocated from heap
352 *
353 ***************************************************************/
354void kb_p_Minus_m_Mult_q_General (poly *pp, int *lpp, 
355                                  poly m,
356                                  poly q, int lq,
357                                  poly spNoether,
358                                  omBin heap)
359{ 
360#ifdef KB_USE_HEAPS
361  assume(heap != NULL);
362#else
363  assume(heap == NULL);
364#endif
365#ifdef HAVE_LENGTH
366  assume(pLength(q) == lq);
367  assume(pLength(*pp) == *lpp);
368  pTest(*pp);
369  pTest(q);
370#endif
371  assume(pGetComp(m) == 0);
372  // we are done if q == NULL
373  if (q == NULL || m == NULL) return;
374 
375  poly a = m,                         // collects the result
376       qm = NULL,                     // stores q*m
377       c,                             // used for temporary storage
378       p = *pp;
379
380  number tm   = pGetCoeff(m),       // coefficient of m
381         tneg = nNeg(nCopy(tm)),    // - (coefficient of m)
382         tb,                        // used for tm*coeff(a1)
383         tc;                        // used as intermediate number
384
385  unsigned int lp = *lpp + lq;
386
387  if (p == NULL) goto Finish; // we are done if p is 0
388
389  kb_pNew(qm, heap);
390  memaddW((unsigned long*) &(qm->exp.l[0]),
391          (unsigned long*) &(q->exp.l[0]),
392          (unsigned long*) &(m->exp.l[0]),
393          currRing->ExpLSize);
394 
395  // MAIN LOOP:
396  Top:     // compare qm = m*q and p w.r.t. monomial ordering
397//    register long d;
398  long d;
399#ifdef HAVE_COMP_MACROS
400    _pMonComp(qm, p, d, goto NotEqual, goto Equal);
401#else
402    if ((d = pComp0(qm, p))) goto NotEqual; else goto Equal;
403#endif
404
405  Equal:   // qm equals p
406    tb = nMult(pGetCoeff(q), tm);
407    tc = pGetCoeff(p);
408    if (!nEqual(tc, tb))
409    {
410      lp--;
411      tc = nSub(tc, tb);
412      nDelete(&(pGetCoeff(p)));
413      pSetCoeff0(p,tc); // adjust coeff of p
414      a = pNext(a) = p; // append p to result and advance p
415      pIter(p);
416    }
417    else
418    { // coeffs are equal, so their difference is 0:
419      lp -= 2;
420      nDelete(&tc);
421      kb_pFree1AndAdvance(p, heap);
422    }
423    nDelete(&tb);
424    pIter(q);
425    if (q == NULL || p == NULL) goto Finish; // are we done ?
426    // no, so update qm
427    memaddW((unsigned long*) &(qm->exp.l[0]),
428            (unsigned long*) &(q->exp.l[0]),
429            (unsigned long*) &(m->exp.l[0]),
430            currRing->ExpLSize);
431    goto Top;
432
433  NotEqual:     // qm != p
434    if (d < 0)  // qm < p:
435    {
436      a = pNext(a) = p;// append p to result and advance p
437      pIter(p);
438      if (p == NULL) goto Finish;;
439      goto Top;
440    }
441    else // now d >= 0, i.e., qm > p
442    {
443      pSetCoeff0(qm,nMult(pGetCoeff(q), tneg));
444      a = pNext(a) = qm;       // append qm to result and advance q
445      pIter(q);
446      if (q == NULL) // are we done?
447      {
448        qm = NULL;
449        goto Finish; 
450      }
451      // construct new qm
452      kb_pNew(qm, heap);
453      memaddW((unsigned long*) &(qm->exp.l[0]),
454              (unsigned long*) &(q->exp.l[0]),
455              (unsigned long*) &(m->exp.l[0]),
456              currRing->ExpLSize);
457      goto Top;
458    }
459 
460 Finish: // q or p is NULL: Clean-up time
461   if (q == NULL) // append rest of p to result
462     pNext(a) = p;
463   else  // append (- q*m) to result
464   {
465     pSetCoeff0(m, tneg);
466     pNext(a) = kb_p_Mult_m(q, m, spNoether, heap);
467     pSetCoeff0(m, tm);
468   }
469   
470   nDelete(&tneg);
471   if (qm != NULL) kb_pFree1(qm, heap);
472
473   *pp = pNext(m);
474   *lpp = lp;
475   pNext(m) = NULL;
476#ifdef HAVE_LENGTH
477   assume(pLength(*pp) == *lpp);
478   pTest(*pp);
479#endif
480} 
481
482
483/***************************************************************
484 *
485 * fast poly proc business
486 *
487 ***************************************************************/
488
489#ifdef FAST_POLY_PROCS
490#define NonZeroA(d, multiplier, actionE)        \
491{                                               \
492  d ^= multiplier;                              \
493  actionE;                                      \
494}                                               \
495
496#define NonZeroTestA(d, multiplier, actionE)    \
497do                                              \
498{                                               \
499  if (d)                                        \
500  {                                             \
501    d ^= multiplier;                            \
502    actionE;                                    \
503  }                                             \
504}                                               \
505while(0)
506
507#include "kbPolyProcs.pin"
508#endif
509
510void kbSetPolyProcs(kbPolyProcs_pt pprocs,
511                    ring r, rOrderType_t rot, BOOLEAN homog)
512{
513  assume(pprocs != NULL);
514#ifdef FAST_POLY_PROCS
515  Characteristics ch = chGEN;
516  OrderingTypes ot = otGEN;
517  Homogs hom = homGEN;
518  NumWords nw = nwGEN;
519  int Variables1W;
520
521  // set characterisic
522  if (r->ch > 1) ch = chMODP;
523 
524  // set Ordering Type
525  switch (rot)
526  {   
527      case rOrderType_Exp:
528        ot = otEXP;
529        break;
530   
531      case rOrderType_CompExp:
532        ot = otCOMPEXP;
533        break;
534   
535      case rOrderType_ExpComp:
536        ot = otEXPCOMP;
537        break;
538
539      case rOrderType_Syz2dpc:
540        ot = otSYZDPC;
541        break;
542       
543      case rOrderType_ExpNoComp:
544        ot = otExpNoComp;
545       
546      default:
547        ot = otGEN;
548        break;
549  }
550 
551  // set homogenous
552  if (homog) hom = homYES;
553 
554  // set NumWords
555  if ((((r->N+1)*sizeof(Exponent_t)) % sizeof(void*)) == 0)
556    Variables1W = (r->N+1)*sizeof(Exponent_t) / sizeof(void*);
557  else
558    Variables1W = ((r->N+1)*sizeof(Exponent_t) / sizeof(void*)) + 1;
559  if (Variables1W > 2)
560  {
561    if (Variables1W & 1) nw = nwODD;
562    else nw = nwEVEN;
563  }
564  else 
565  {
566    if (Variables1W == 2) nw = nwTWO;
567    else nw = nwONE;
568  }
569 
570  // Get the nPoly Procs
571  pprocs->p_Add_q = Getkb_p_Add_q(ch, ot, hom, nw);
572  if (pprocs->p_Add_q == NULL)
573    pprocs->p_Add_q = kb_p_Add_q_General;
574
575  pprocs->p_Mult_m = Getkb_p_Mult_m(ch, ot, hom, nw);
576  if (pprocs->p_Mult_m == NULL)
577    pprocs->p_Mult_m = kb_p_Mult_m_General;
578
579  pprocs->p_Minus_m_Mult_q = Getkb_p_Minus_m_Mult_q(ch, ot, hom, nw);
580  if (pprocs->p_Minus_m_Mult_q == NULL)
581    pprocs->p_Minus_m_Mult_q = kb_p_Minus_m_Mult_q_General;
582
583  pprocs->n_Mult_p = Getkb_n_Mult_p(ch, ot, hom, nw);
584  if (pprocs->n_Mult_p == NULL)
585    pprocs->n_Mult_p = kb_n_Mult_p_General;
586#else // FAST_POLYPROCS
587  pprocs->p_Add_q = kb_p_Add_q_General;
588  pprocs->p_Mult_m = kb_p_Mult_m_General;
589  pprocs->p_Minus_m_Mult_q = kb_p_Minus_m_Mult_q_General;
590  pprocs->n_Mult_p = kb_n_Mult_p_General;
591#endif 
592}
Note: See TracBrowser for help on using the repository browser.