source: git/Singular/kbPolyProcs.cc @ c06a32

spielwiese
Last change on this file since c06a32 was c06a32, checked in by Olaf Bachmann <obachman@…>, 25 years ago
* New Handling of Command-line options git-svn-id: file:///usr/local/Singular/svn/trunk@3623 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.2 KB
Line 
1#include "mod2.h"
2#include "mmemory.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 "modulop.h"
10#include "numbers.h"
11#include "polys-comp.h"
12#include "kbPolyProcs.h"
13
14// define to enable fast poly procs
15// #define FAST_POLY_PROCS
16
17/***************************************************************
18 *
19 * General kb_n_Mult_p which always works
20 * Realizes p = n*p
21 *
22 ***************************************************************/
23
24void kb_n_Mult_p_General(number n, poly p)
25{
26  while (p != NULL)
27  {
28    number nc = pGetCoeff(p);
29    pSetCoeff0(p, nMult(n, nc));
30    nDelete(&nc);
31    pIter(p);
32  }
33}
34
35/***************************************************************
36 *
37 * General kb_p_Add_q which always works
38 *
39 * assume pLength(*p) == *lp && pLength(*q) == *q
40 * *p and *q from heap
41 * Destroys *p and *q
42 * On return, *p == *p + *q, *q == NULL,
43 *            *lp == pLength(*p + *q), *lq == 0
44 *
45 ***************************************************************/
46
47void kb_p_Add_q_General(poly *p, int *lp,
48                        poly *q, int *lq, 
49                        memHeap heap)
50{
51#ifdef KB_USE_HEAPS
52  assume(heap != NULL);
53#else
54  assume(heap == NULL);
55#endif
56  assume(pLength(*p) == *lp && pLength(*q) == *lq);
57 
58  number t, n1, n2;
59  unsigned int l = *lp + *lq;
60  spolyrec rp;
61  poly a = &rp, a1 = *p, a2 = *q;
62
63  if (a2 == NULL) return;
64
65  *q = NULL;
66  *lq = 0;
67 
68  if (a1 == NULL) 
69  {
70    *p = a2;
71    *lp = l;
72    return;
73  }
74
75  Top:     // compare a1 and a2 w.r.t. monomial ordering
76  register long d;
77  if ((d = pComp0(a1, a2))) goto NotEqual; else goto Equal;
78
79  Equal:
80  assume(pComp0(a1, a2) == 0);
81
82  n1 = pGetCoeff(a1);
83  n2 = pGetCoeff(a2);
84  t = nAdd(n1,n2);
85  nDelete(&n1);
86  nDelete(&n2);
87  kb_pFree1AndAdvance(a2, heap);
88 
89  if (nIsZero(t))
90  {
91    l -= 2;
92    nDelete(&t);
93    kb_pFree1AndAdvance(a1, heap);
94  }
95  else
96  {
97    l--;
98    pSetCoeff0(a1,t);
99    a = pNext(a) = a1;
100    pIter(a1);
101  }
102  if (a1==NULL)
103  {
104    pNext(a) = a2;
105    goto Finish;
106  }
107  else if (a2==NULL)
108  {
109    pNext(a) = a1;
110    goto Finish;
111  }
112  goto Top;
113     
114  NotEqual:
115  if (d < 0)
116  {
117    assume(pComp0(a1, a2) == -1);
118    a = pNext(a) = a2;
119    pIter(a2);
120    if (a2==NULL)
121    {
122      pNext(a) = a1;
123      goto Finish;
124    }
125  }
126  else 
127  {
128    assume(pComp0(a1, a2) == 1);
129    a = pNext(a) = a1;
130    pIter(a1);
131    if (a1==NULL)
132    {
133      pNext(a) = a2;
134      goto Finish;
135    }
136  }
137  goto Top;
138 
139
140  Finish: 
141  assume(pLength(pNext(&rp)) == (int) l);
142  *lp = l;
143  *p  = pNext(&rp);
144}
145
146/***************************************************************
147 *
148 * General kb_p_Mult_m which always works
149 *
150 * assume pLength(m) == 1
151 * Returns m*a1, monoms are allocated from heap
152 *
153 ***************************************************************/
154 
155poly  kb_p_Mult_m_General(poly p,
156                          poly m, 
157                          poly spNoether,
158                          memHeap heap)
159{
160#ifdef KB_USE_HEAPS
161  assume(heap != NULL);
162#else
163  assume(heap == NULL);
164#endif
165  if (p == NULL) return NULL;
166  spolyrec rp;
167  poly q = &rp;
168  number ln = pGetCoeff(m);
169  int comp = pGetComp(m);
170
171  pSetComp(m, 0);
172
173  while (p != NULL)
174  {
175    kb_pNew(pNext(q), heap);
176    q = pNext(q);
177
178    pSetCoeff0(q, nMult(ln, pGetCoeff(p)));
179
180    q->Order = p->Order + m->Order;
181    memaddW((unsigned long*) &(q->exp[0]),
182            (unsigned long*) &(p->exp[0]),
183            (unsigned long*) &(m->exp[0]),
184            pVariables1W);
185   
186    p = pNext(p);
187  }
188  pNext(q) = NULL;
189  pSetComp(m, comp);
190  return rp.next;
191}
192
193/***************************************************************
194 *
195 * General spoly loop which always works
196 *
197 * assume(pLength(*pp) == *lpp)
198 * assume(pLength(q) == lq)
199 * assume(pLength(m) == NULL)
200 * assume: Monoms of *pp are from heap
201 *
202 * Destroys *pp
203 * Does not touch q
204 * On return *pp == *pp - m*q
205 *           *lpp == pLength(*pp - m*q)
206 * New monoms are allocated from heap
207 *
208 ***************************************************************/
209void kb_p_Minus_m_Mult_q_General (poly *pp, int *lpp, 
210                                  poly m,
211                                  poly q, int lq,
212                                  poly spNoether,
213                                  kb_p_Mult_m_Proc kb_p_Mult_m,
214                                  memHeap heap)
215{ 
216#ifdef KB_USE_HEAPS
217  assume(heap != NULL);
218#else
219  assume(heap == NULL);
220#endif
221  assume(pLength(q) == lq);
222  assume(pLength(*pp) == *lpp);
223 
224  // we are done if q == NULL
225  if (q == NULL || m == NULL) return;
226 
227  poly a = m,                         // collects the result
228       qm = NULL,                     // stores q*m
229       c,                             // used for temporary storage
230       p = *pp;
231
232  number tm   = pGetCoeff(m),       // coefficient of m
233         tneg = nNeg(nCopy(tm)),    // - (coefficient of m)
234         tb,                        // used for tm*coeff(a1)
235         tc;                        // used as intermediate number
236
237  unsigned int lp = *lpp + lq;
238
239  int comp = pGetComp(m);
240  pSetComp(m, 0);
241 
242  if (p == NULL) goto Finish; // we are done if p is 0
243
244  kb_pNew(qm, heap);
245  qm->Order = q->Order + m->Order;
246  memaddW((unsigned long*) &(qm->exp[0]),
247          (unsigned long*) &(q->exp[0]),
248          (unsigned long*) &(m->exp[0]),
249          pVariables1W);
250 
251  // MAIN LOOP:
252  Top:     // compare qm = m*q and p w.r.t. monomial ordering
253    register long d;
254    if ((d = pComp0(qm, p))) goto NotEqual; else goto Equal;
255
256  Equal:   // qm equals p
257    tb = nMult(pGetCoeff(q), tm);
258    tc = pGetCoeff(p);
259    if (!nEqual(tc, tb))
260    {
261      lp--;
262      tc = nSub(tc, tb);
263      nDelete(&(pGetCoeff(p)));
264      pSetCoeff0(p,tc); // adjust coeff of p
265      a = pNext(a) = p; // append p to result and advance p
266      pIter(p);
267    }
268    else
269    { // coeffs are equal, so their difference is 0:
270      lp -= 2;
271      nDelete(&tc);
272      kb_pFree1AndAdvance(p, heap);
273    }
274    nDelete(&tb);
275    pIter(q);
276    if (q == NULL || p == NULL) goto Finish; // are we done ?
277    // no, so update qm
278    qm->Order = q->Order + m->Order;
279    memaddW((unsigned long*) &(qm->exp[0]),
280            (unsigned long*) &(q->exp[0]),
281            (unsigned long*) &(m->exp[0]),
282            pVariables1W);
283    goto Top;
284
285  NotEqual:     // qm != p
286    if (d < 0)  // qm < p:
287    {
288      a = pNext(a) = p;// append p to result and advance p
289      pIter(p);
290      if (p == NULL) goto Finish;;
291      goto Top;
292    }
293    else // now d >= 0, i.e., qm > p
294    {
295      pSetCoeff0(qm,nMult(pGetCoeff(q), tneg));
296      a = pNext(a) = qm;       // append qm to result and advance q
297      pIter(q);
298      if (q == NULL) // are we done?
299      {
300        qm = NULL;
301        goto Finish; 
302      }
303      // construct new qm
304      kb_pNew(qm, heap);
305      qm->Order = q->Order + m->Order;
306      memaddW((unsigned long*) &(qm->exp[0]),
307              (unsigned long*) &(q->exp[0]),
308              (unsigned long*) &(m->exp[0]),
309              pVariables1W);
310      goto Top;
311    }
312 
313 Finish: // q or p is NULL: Clean-up time
314   pSetComp(m, comp);
315   if (q == NULL) // append rest of p to result
316     pNext(a) = p;
317   else  // append (- q*m) to result
318   {
319     pSetCoeff0(m, tneg);
320     pNext(a) = kb_p_Mult_m(q, m, spNoether, heap);
321     pSetCoeff0(m, tm);
322   }
323   
324   nDelete(&tneg);
325   if (qm != NULL) kb_pFree1(qm, heap);
326
327   *pp = pNext(m);
328   *lpp = lp;
329   pNext(m) = NULL;
330   
331   assume(pLength(*pp) == *lpp);
332   pHeapTest(*pp, heap);
333} 
334
335
336/***************************************************************
337 *
338 * fast poly proc business
339 *
340 ***************************************************************/
341
342#ifdef FAST_POLY_PROCS
343#define NonZeroA(d, multiplier, actionE)        \
344{                                               \
345  d ^= multiplier;                              \
346  actionE;                                      \
347}                                               \
348
349#define NonZeroTestA(d, multiplier, actionE)    \
350do                                              \
351{                                               \
352  if (d)                                        \
353  {                                             \
354    d ^= multiplier;                            \
355    actionE;                                    \
356  }                                             \
357}                                               \
358while(0)
359
360#include "kbPolyProcs.pin"
361#endif
362
363void kbSetPolyProcs(kbPolyProcs_pt pprocs,
364                    ring r, rOrderType_t rot, BOOLEAN homog)
365{
366  assume(pprocs != NULL);
367#ifdef FAST_POLY_PROCS
368  Characteristics ch = chGEN;
369  OrderingTypes ot = otGEN;
370  Homogs hom = homGEN;
371  NumWords nw = nwGEN;
372  int Variables1W;
373
374  // set characterisic
375  if (rField_is_Zp(r)) ch = chMODP;
376 
377  // set Ordering Type
378  switch (rot)
379  {   
380      case rOrderType_Exp:
381        ot = otEXP;
382        break;
383   
384      case rOrderType_CompExp:
385        ot = otCOMPEXP;
386        break;
387   
388      case rOrderType_ExpComp:
389        ot = otEXPCOMP;
390        break;
391
392      case rOrderType_Syz2dpc:
393        ot = otSYZDPC;
394        break;
395       
396      default:
397        ot = otGEN;
398        break;
399  }
400 
401  // set homogenous
402  if (homog) hom = homYES;
403 
404  // set NumWords
405  if ((((r->N+1)*sizeof(Exponent_t)) % sizeof(void*)) == 0)
406    Variables1W = (r->N+1)*sizeof(Exponent_t) / sizeof(void*);
407  else
408    Variables1W = ((r->N+1)*sizeof(Exponent_t) / sizeof(void*)) + 1;
409  if (Variables1W > 2)
410  {
411    if (Variables1W & 1) nw = nwODD;
412    else nw = nwEVEN;
413  }
414  else 
415  {
416    if (Variables1W == 2) nw = nwTWO;
417    else nw = nwONE;
418  }
419 
420  // Get the nPoly Procs
421  pprocs->p_Add_q = Getkb_p_Add_q(ch, ot, hom, nw);
422  if (pprocs->p_Add_q == NULL)
423    pprocs->p_Add_q = kb_p_Add_q_General;
424
425  pprocs->p_Mult_m = Getkb_p_Mult_m(ch, ot, hom, nw);
426  if (pprocs->p_Mult_m == NULL)
427    pprocs->p_Mult_m = kb_p_Mult_m_General;
428
429  pprocs->p_Minus_m_Mult_q = Getkb_p_Minus_m_Mult_q(ch, ot, hom, nw);
430  if (pprocs->p_Minus_m_Mult_q == NULL)
431    pprocs->p_Minus_m_Mult_q = kb_p_Minus_m_Mult_q_General;
432
433  pprocs->n_Mult_p = Getkb_n_Mult_p(ch, ot, hom, nw);
434  if (pprocs->n_Mult_p == NULL)
435    pprocs->n_Mult_p = kb_n_Mult_p_General;
436#else // FAST_POLYPROCS
437  pprocs->p_Add_q = kb_p_Add_q_General;
438  pprocs->p_Mult_m = kb_p_Mult_m_General;
439  pprocs->p_Minus_m_Mult_q = kb_p_Minus_m_Mult_q_General;
440  pprocs->n_Mult_p = kb_n_Mult_p_General;
441#endif 
442}
Note: See TracBrowser for help on using the repository browser.