source: git/kernel/pInline2.h @ 321283

spielwiese
Last change on this file since 321283 was 5be91d, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: fix nc_* git-svn-id: file:///usr/local/Singular/svn/trunk@9255 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    pInline2.h
6 *  Purpose: implementation of poly procs which are of constant time
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 8/00
9 *  Version: $Id: pInline2.h,v 1.5 2006-06-22 16:59:13 Singular Exp $
10 *******************************************************************/
11#ifndef PINLINE2_H
12#define PINLINE2_H
13
14/***************************************************************
15 *
16 * Primitives for accessing and setting fields of a poly
17 *
18 ***************************************************************/
19#if !defined(NO_PINLINE2) || defined(PINLINE2_CC)
20
21#include "structs.h"
22#include "omalloc.h"
23#include "numbers.h"
24#include "p_Procs.h"
25#include "sbuckets.h"
26#ifdef HAVE_PLURAL
27#include "gring.h"
28#include "ring.h"
29#endif
30
31PINLINE2 number p_SetCoeff(poly p, number n, ring r)
32{
33  p_LmCheckPolyRing2(p, r);
34  n_Delete(&(p->coef), r);
35  (p)->coef=n;
36  return n;
37}
38
39// order
40PINLINE2 Order_t p_GetOrder(poly p, ring r)
41{
42  p_LmCheckPolyRing2(p, r);
43  if (r->typ==NULL) return ((p)->exp[r->pOrdIndex]);
44  int i=0;
45  loop
46  {
47    switch(r->typ[i].ord_typ)
48    {
49      case ro_wp_neg:
50        return (((long)((p)->exp[r->pOrdIndex]))-POLY_NEGWEIGHT_OFFSET);
51      case ro_syzcomp:
52      case ro_syz:
53      case ro_cp:
54        i++;
55        break;
56      //case ro_dp:
57      //case ro_wp:
58      default:
59        return ((p)->exp[r->pOrdIndex]);
60    }
61  }
62}
63
64PINLINE2 Order_t p_SetOrder(poly p, long o, ring r)
65{
66  p_LmCheckPolyRing2(p, r);
67  pAssume2(o >= 0);
68  if (r->typ==NULL) return ((p)->exp[r->pOrdIndex]=o);
69  int i=0;
70  loop
71  {
72    switch(r->typ[i].ord_typ)
73    {
74      case ro_wp_neg:
75        return (p)->exp[r->pOrdIndex]=o+POLY_NEGWEIGHT_OFFSET;
76      case ro_syzcomp:
77      case ro_syz:
78      case ro_cp:
79        i++;
80        break;
81      //case ro_dp:
82      //case ro_wp:
83      default:
84        return (p)->exp[r->pOrdIndex] = o;
85    }
86  }
87}
88
89// Setm
90PINLINE2 void p_Setm(poly p, ring r)
91{
92  p_CheckRing2(r);
93  r->p_Setm(p, r);
94}
95
96// component
97PINLINE2  unsigned long p_SetComp(poly p, unsigned long c, ring r)
98{
99  p_LmCheckPolyRing2(p, r);
100  pAssume2(rRing_has_Comp(r));
101  __p_GetComp(p,r) = c;
102  return c;
103}
104PINLINE2 unsigned long p_IncrComp(poly p, ring r)
105{
106  p_LmCheckPolyRing2(p, r);
107  pAssume2(rRing_has_Comp(r));
108  return ++(__p_GetComp(p,r));
109}
110PINLINE2 unsigned long p_DecrComp(poly p, ring r)
111{
112  p_LmCheckPolyRing2(p, r);
113  pAssume2(rRing_has_Comp(r));
114  pPolyAssume2(__p_GetComp(p,r) > 0);
115  return --(__p_GetComp(p,r));
116}
117PINLINE2 unsigned long p_AddComp(poly p, unsigned long v, ring r)
118{
119  p_LmCheckPolyRing2(p, r);
120  pAssume2(rRing_has_Comp(r));
121  return __p_GetComp(p,r) += v;
122}
123PINLINE2 unsigned long p_SubComp(poly p, unsigned long v, ring r)
124{
125  p_LmCheckPolyRing2(p, r);
126  pAssume2(rRing_has_Comp(r));
127  pPolyAssume2(__p_GetComp(p,r) >= v);
128  return __p_GetComp(p,r) -= v;
129}
130
131// exponent
132// r->VarOffset encodes the position in p->exp (lower 24 bits)
133// and number of bits to shift to the right in the upper 8 bits
134PINLINE2 int p_GetExp(poly p, int v, ring r)
135{
136  p_LmCheckPolyRing2(p, r);
137  pAssume2(v > 0 && v <= r->N);
138#if 0
139  int pos=(r->VarOffset[v] & 0xffffff);
140  int bitpos=(r->VarOffset[v] >> 24);
141  int exp=(p->exp[pos] >> bitmask) & r->bitmask;
142  return exp;
143#else
144  return (int)
145         ((p->exp[(r->VarOffset[v] & 0xffffff)] >> (r->VarOffset[v] >> 24))
146          & r->bitmask);
147#endif
148}
149PINLINE2 int p_SetExp(poly p, int v, int e, ring r)
150{
151  p_LmCheckPolyRing2(p, r);
152  pAssume2(v>0 && v <= r->N);
153  pAssume2(e>=0);
154  pAssume2((unsigned int) e<=r->bitmask);
155
156  // shift e to the left:
157  register int shift = r->VarOffset[v] >> 24;
158  unsigned long ee = ((unsigned long)e) << shift /*(r->VarOffset[v] >> 24)*/;
159  // find the bits in the exponent vector
160  register int offset = (r->VarOffset[v] & 0xffffff);
161  // clear the bits in the exponent vector:
162  p->exp[offset]  &= ~( r->bitmask << shift );
163  // insert e with |
164  p->exp[ offset ] |= ee;
165  return e;
166}
167
168// the following should be implemented more efficiently
169PINLINE2  int p_IncrExp(poly p, int v, ring r)
170{
171  p_LmCheckPolyRing2(p, r);
172  int e = p_GetExp(p,v,r);
173  e++;
174  return p_SetExp(p,v,e,r);
175}
176PINLINE2  int p_DecrExp(poly p, int v, ring r)
177{
178  p_LmCheckPolyRing2(p, r);
179  int e = p_GetExp(p,v,r);
180  pAssume2(e > 0);
181  e--;
182  return p_SetExp(p,v,e,r);
183}
184PINLINE2  int p_AddExp(poly p, int v, int ee, ring r)
185{
186  p_LmCheckPolyRing2(p, r);
187  int e = p_GetExp(p,v,r);
188  e += ee;
189  return p_SetExp(p,v,e,r);
190}
191PINLINE2  int p_SubExp(poly p, int v, int ee, ring r)
192{
193  p_LmCheckPolyRing2(p, r);
194  int e = p_GetExp(p,v,r);
195  pAssume2(e >= ee);
196  e -= ee;
197  return p_SetExp(p,v,e,r);
198}
199PINLINE2  int p_MultExp(poly p, int v, int ee, ring r)
200{
201  p_LmCheckPolyRing2(p, r);
202  int e = p_GetExp(p,v,r);
203  e *= ee;
204  return p_SetExp(p,v,e,r);
205}
206
207PINLINE2 int p_GetExpSum(poly p1, poly p2, int i, ring r)
208{
209  p_LmCheckPolyRing2(p1, r);
210  p_LmCheckPolyRing2(p2, r);
211  return p_GetExp(p1,i,r) + p_GetExp(p2,i,r);
212}
213PINLINE2 int p_GetExpDiff(poly p1, poly p2, int i, ring r)
214{
215  return p_GetExp(p1,i,r) - p_GetExp(p2,i,r);
216}
217
218
219/***************************************************************
220 *
221 * Allocation/Initalization/Deletion
222 *
223 ***************************************************************/
224PINLINE2 poly p_New(ring r, omBin bin)
225{
226  p_CheckRing2(r);
227  pAssume2(bin != NULL && r->PolyBin->sizeW == bin->sizeW);
228  poly p;
229  omTypeAllocBin(poly, p, bin);
230  p_SetRingOfLm(p, r);
231  return p;
232}
233
234PINLINE2 poly p_New(ring r)
235{
236  return p_New(r, r->PolyBin);
237}
238
239PINLINE2 void p_DeleteLm(poly *p, ring r)
240{
241  pIfThen2(*p != NULL, p_LmCheckPolyRing2(*p, r));
242  poly h = *p;
243  if (h != NULL)
244  {
245    n_Delete(&_pGetCoeff(h), r);
246    *p = _pNext(h);
247    omFreeBinAddr(h);
248  }
249}
250PINLINE2 void p_DeleteLm(poly p, ring r)
251{
252  pIfThen2(p != NULL, p_LmCheckPolyRing2(p, r));
253  if (p != NULL)
254  {
255    n_Delete(&_pGetCoeff(p), r);
256    omFreeBinAddr(p);
257  }
258}
259PINLINE2 void p_LmFree(poly p, ring r)
260{
261  p_LmCheckPolyRing2(p, r);
262  omFreeBinAddr(p);
263}
264PINLINE2 void p_LmFree(poly *p, ring r)
265{
266  p_LmCheckPolyRing2(*p, r);
267  poly h = *p;
268  *p = pNext(h);
269  omFreeBinAddr(h);
270}
271PINLINE2 poly p_LmFreeAndNext(poly p, ring r)
272{
273  p_LmCheckPolyRing2(p, r);
274  poly pnext = pNext(p);
275  omFreeBinAddr(p);
276  return pnext;
277}
278PINLINE2 void p_LmDelete(poly p, ring r)
279{
280  p_LmCheckPolyRing2(p, r);
281  n_Delete(&_pGetCoeff(p), r);
282  omFreeBinAddr(p);
283}
284PINLINE2 void p_LmDelete(poly *p, ring r)
285{
286  p_LmCheckPolyRing2(*p, r);
287  poly h = *p;
288  *p = pNext(h);
289  n_Delete(&pGetCoeff(h), r);
290  omFreeBinAddr(h);
291}
292PINLINE2 poly p_LmDeleteAndNext(poly p, ring r)
293{
294  p_LmCheckPolyRing2(p, r);
295  poly pnext = _pNext(p);
296  n_Delete(&_pGetCoeff(p), r);
297  omFreeBinAddr(p);
298  return pnext;
299}
300
301/***************************************************************
302 *
303 * Misc routines
304 *
305 ***************************************************************/
306PINLINE2 int p_Cmp(poly p1, poly p2, ring r)
307{
308  if (p2==NULL)
309    return 1;
310  if (p1==NULL)
311    return -1;
312  return p_LmCmp(p1,p2,r);
313}
314
315PINLINE2 unsigned long p_GetMaxExp(poly p, ring r)
316{
317  return p_GetMaxExp(p_GetMaxExpL(p, r), r);
318}
319
320PINLINE2 unsigned long
321p_GetMaxExp(const unsigned long l, const ring r, const int number_of_exps)
322{
323  unsigned long bitmask = r->bitmask;
324  unsigned long max = (l & bitmask);
325  unsigned long j = number_of_exps - 1;
326
327  if (j > 0)
328  {
329    unsigned long i = r->BitsPerExp;
330    long e;
331    loop
332    {
333      e = ((l >> i) & bitmask);
334      if ((unsigned long) e > max)
335        max = e;
336      j--;
337      if (j==0) break;
338      i += r->BitsPerExp;
339    }
340  }
341  return max;
342}
343
344PINLINE2 unsigned long p_GetMaxExp(const unsigned long l, const ring r)
345{
346  return p_GetMaxExp(l, r, r->ExpPerLong);
347}
348
349PINLINE2 unsigned long
350p_GetTotalDegree(const unsigned long l, const ring r, const int number_of_exps)
351{
352  const unsigned long bitmask = r->bitmask;
353  unsigned long sum = (l & bitmask);
354  unsigned long j = number_of_exps - 1;
355
356  if (j > 0)
357  {
358    unsigned long i = r->BitsPerExp;
359    loop
360    {
361      sum += ((l >> i) & bitmask);
362      j--;
363      if (j==0) break;
364      i += r->BitsPerExp;
365    }
366  }
367  return sum;
368}
369
370PINLINE2 unsigned long
371p_GetTotalDegree(const unsigned long l, const ring r)
372{
373  return p_GetTotalDegree(l, r, r->ExpPerLong);
374}
375
376/***************************************************************
377 *
378 * Dispatcher to r->p_Procs, they do the tests/checks
379 *
380 ***************************************************************/
381// returns a copy of p
382PINLINE2 poly p_Copy(poly p, const ring r)
383{
384#ifdef PDEBUG
385  poly pp= r->p_Procs->p_Copy(p, r);
386  p_Test(pp,r);
387  return pp;
388#else
389  return r->p_Procs->p_Copy(p, r);
390#endif
391}
392
393PINLINE2 poly p_Copy(poly p, const ring lmRing, const ring tailRing)
394{
395#ifndef PDEBUG
396  if (tailRing == lmRing)
397    return tailRing->p_Procs->p_Copy(p, tailRing);
398#endif
399  if (p != NULL)
400  {
401    poly pres = p_Head(p, lmRing);
402    pNext(pres) = tailRing->p_Procs->p_Copy(pNext(p), tailRing);
403    return pres;
404  }
405  else
406    return NULL;
407}
408
409// deletes *p, and sets *p to NULL
410PINLINE2 void p_Delete(poly *p, const ring r)
411{
412  r->p_Procs->p_Delete(p, r);
413}
414
415PINLINE2 void p_Delete(poly *p,  const ring lmRing, const ring tailRing)
416{
417#ifndef PDEBUG
418  if (tailRing == lmRing)
419  {
420    tailRing->p_Procs->p_Delete(p, tailRing);
421    return;
422  }
423#endif
424  if (*p != NULL)
425  {
426    if (pNext(*p) != NULL)
427      tailRing->p_Procs->p_Delete(&pNext(*p), tailRing);
428    p_LmDelete(p, lmRing);
429  }
430}
431
432PINLINE2 poly p_ShallowCopyDelete(poly p, const ring r, omBin bin)
433{
434  p_LmCheckPolyRing2(p, r);
435  pAssume2(r->PolyBin->sizeW == bin->sizeW);
436  return r->p_Procs->p_ShallowCopyDelete(p, r, bin);
437}
438
439// returns p+q, destroys p and q
440PINLINE2 poly p_Add_q(poly p, poly q, const ring r)
441{
442  int shorter;
443  return r->p_Procs->p_Add_q(p, q, shorter, r);
444}
445
446PINLINE2 poly p_Add_q(poly p, poly q, int &lp, int lq, const ring r)
447{
448  int shorter;
449  poly res = r->p_Procs->p_Add_q(p, q, shorter, r);
450  lp = (lp + lq) - shorter;
451  return res;
452}
453
454// returns p*n, destroys p
455PINLINE2 poly p_Mult_nn(poly p, number n, const ring r)
456{
457  if (n_IsOne(n, r))
458    return p;
459  else
460    return r->p_Procs->p_Mult_nn(p, n, r);
461}
462
463PINLINE2 poly p_Mult_nn(poly p, number n, const ring lmRing,
464                        const ring tailRing)
465{
466#ifndef PDEBUG
467  if (lmRing == tailRing)
468  {
469    return p_Mult_nn(p, n, tailRing);
470  }
471#endif
472  poly pnext = pNext(p);
473  pNext(p) = NULL;
474  p = lmRing->p_Procs->p_Mult_nn(p, n, lmRing);
475  pNext(p) = tailRing->p_Procs->p_Mult_nn(pnext, n, tailRing);
476  return p;
477}
478
479// returns p*n, does not destroy p
480PINLINE2 poly pp_Mult_nn(poly p, number n, const ring r)
481{
482  if (n_IsOne(n, r))
483    return p_Copy(p, r);
484  else
485    return r->p_Procs->pp_Mult_nn(p, n, r);
486}
487
488// returns Copy(p)*m, does neither destroy p nor m
489PINLINE2 poly pp_Mult_mm(poly p, poly m, const ring r)
490{
491  if (p_LmIsConstant(m, r))
492    return pp_Mult_nn(p, pGetCoeff(m), r);
493  else
494  {
495    poly last;
496    return r->p_Procs->pp_Mult_mm(p, m, r, last);
497  }
498}
499
500// returns p*m, destroys p, const: m
501PINLINE2 poly p_Mult_mm(poly p, poly m, const ring r)
502{
503  if (p_LmIsConstant(m, r))
504    return p_Mult_nn(p, pGetCoeff(m), r);
505  else
506    return r->p_Procs->p_Mult_mm(p, m, r);
507}
508
509// return p - m*Copy(q), destroys p; const: p,m
510PINLINE2 poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
511{
512  int shorter;
513  poly last;
514  return r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last);
515}
516PINLINE2 poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
517                                 poly spNoether, const ring r)
518{
519  int shorter;
520  poly last,res;
521#ifdef HAVE_PLURAL
522  if (rIsPluralRing(r))
523  {
524     res = nc_p_Minus_mm_Mult_qq(p, m, q, r);
525     lp = pLength(res);
526  }
527  else
528#endif
529  {
530    res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, spNoether, r, last);
531    lp = (lp + lq) - shorter;
532  }
533  return res;
534}
535
536PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
537{
538  int shorter;
539  return r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
540}
541
542PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, int &lp, const poly m, const ring r)
543{
544  int shorter;
545  poly pp = r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
546  lp -= shorter;
547  return pp;
548}
549
550// returns -p, destroys p
551PINLINE2 poly p_Neg(poly p, const ring r)
552{
553  return r->p_Procs->p_Neg(p, r);
554}
555
556extern poly  _p_Mult_q(poly p, poly q, const int copy, const ring r);
557// returns p*q, destroys p and q
558PINLINE2 poly p_Mult_q(poly p, poly q, const ring r)
559{
560  if (p == NULL)
561  {
562    r->p_Procs->p_Delete(&q, r);
563    return NULL;
564  }
565  if (q == NULL)
566  {
567    r->p_Procs->p_Delete(&p, r);
568    return NULL;
569  }
570
571  if (pNext(p) == NULL)
572  {
573#ifdef HAVE_PLURAL
574    if (rIsPluralRing(r))
575      q = nc_mm_Mult_p(p, q, r);
576    else
577#endif /* HAVE_PLURAL */
578      q = r->p_Procs->p_Mult_mm(q, p, r);
579
580    r->p_Procs->p_Delete(&p, r);
581    return q;
582  }
583
584  if (pNext(q) == NULL)
585  {
586  // NEEDED
587#ifdef HAVE_PLURAL
588    if (rIsPluralRing(r))
589      p = nc_p_Mult_mm(p, q, r);
590    else
591#endif /* HAVE_PLURAL */
592      p = r->p_Procs->p_Mult_mm(p, q, r);
593
594    r->p_Procs->p_Delete(&q, r);
595    return p;
596  }
597#ifdef HAVE_PLURAL
598  if (rIsPluralRing(r))
599    return _nc_p_Mult_q(p, q, 0, r);
600  else
601#endif
602  return _p_Mult_q(p, q, 0, r);
603}
604
605// returns p*q, does neither destroy p nor q
606PINLINE2 poly pp_Mult_qq(poly p, poly q, const ring r)
607{
608  poly last;
609  if (p == NULL || q == NULL) return NULL;
610
611  if (pNext(p) == NULL)
612  {
613#ifdef HAVE_PLURAL
614    if (rIsPluralRing(r))
615      return nc_mm_Mult_p(p, p_Copy(q,r), r);
616#endif
617    return r->p_Procs->pp_Mult_mm(q, p, r, last);
618  }
619
620  if (pNext(q) == NULL)
621  {
622    return r->p_Procs->pp_Mult_mm(p, q, r, last);
623  }
624
625  poly qq = q;
626  if (p == q)
627    qq = p_Copy(q, r);
628
629  poly res;
630#ifdef HAVE_PLURAL
631  if (rIsPluralRing(r))
632    res = _nc_p_Mult_q(p, qq, 1, r);
633  else
634#endif
635    res = _p_Mult_q(p, qq, 1, r);
636
637  if (qq != q)
638    p_Delete(&qq, r);
639  return res;
640}
641
642// returns p + m*q destroys p, const: q, m
643// this should be implemented more efficiently
644PINLINE2 poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
645                                const ring r)
646{
647  poly res, last;
648  int shorter;
649  number n_old = pGetCoeff(m);
650  number n_neg = n_Copy(n_old, r);
651  n_neg = n_Neg(n_neg, r);
652  pSetCoeff0(m, n_neg);
653#ifdef HAVE_PLURAL
654  if (rIsPluralRing(r))
655    res = nc_p_Minus_mm_Mult_qq(p, m, q, r);
656  else
657#endif
658    res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last);
659
660  lp = (lp + lq) - shorter;
661  pSetCoeff0(m, n_old);
662  n_Delete(&n_neg, r);
663  return res;
664}
665
666PINLINE2 poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
667{
668  int lp = 0, lq = 0;
669  return p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
670}
671
672PINLINE2 poly p_Merge_q(poly p, poly q, const ring r)
673{
674  return r->p_Procs->p_Merge_q(p, q, r);
675}
676
677PINLINE2 poly p_SortAdd(poly p, const ring r, BOOLEAN revert)
678{
679  if (revert) p = pReverse(p);
680  return sBucketSortAdd(p, r);
681}
682
683PINLINE2 poly p_SortMerge(poly p, const ring r, BOOLEAN revert)
684{
685  if (revert) p = pReverse(p);
686  return sBucketSortMerge(p, r);
687}
688
689/***************************************************************
690 *
691 * I/O
692 *
693 ***************************************************************/
694PINLINE2 char*     p_String(poly p, ring p_ring)
695{
696  return p_String(p, p_ring, p_ring);
697}
698PINLINE2 char*     p_String0(poly p, ring p_ring)
699{
700  return p_String0(p, p_ring, p_ring);
701}
702PINLINE2 void      p_Write(poly p, ring p_ring)
703{
704  p_Write(p, p_ring, p_ring);
705}
706PINLINE2 void      p_Write0(poly p, ring p_ring)
707{
708  p_Write0(p, p_ring, p_ring);
709}
710PINLINE2 void      p_wrp(poly p, ring p_ring)
711{
712  p_wrp(p, p_ring, p_ring);
713}
714#endif // !defined(NO_PINLINE2) || defined(POLYS_IMPL_CC)
715#endif // PINLINE2_H
Note: See TracBrowser for help on using the repository browser.