source: git/kernel/pInline2.h @ 171950

spielwiese
Last change on this file since 171950 was 9fb610, checked in by Viktor Levandovskyy <levandov@…>, 17 years ago
*levandov: updates rational ncGB git-svn-id: file:///usr/local/Singular/svn/trunk@9936 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.2 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.12 2007-03-10 15:41:49 levandov 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, const 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}
149// partial compare exponent
150// r->VarOffset encodes the position in p->exp (lower 24 bits)
151// and number of bits to shift to the right in the upper 8 bits
152PINLINE2 int p_Comp_k_n(poly a, poly b, int k, ring r)
153{
154  if ((a==NULL) || (b==NULL) ) return FALSE;
155  p_LmCheckPolyRing2(a, r);
156  p_LmCheckPolyRing2(b, r);
157  pAssume2(k > 0 && k <= r->N);
158  int i=k;
159  for(;i<=r->N;i++)
160  {
161    if (p_GetExp(a,i,r) != p_GetExp(b,i,r)) return FALSE;
162    //    if (a->exp[(r->VarOffset[i] & 0xffffff)] != b->exp[(r->VarOffset[i] & 0xffffff)]) return FALSE;
163  }
164  return TRUE;
165}
166PINLINE2 int p_SetExp(poly p, int v, int e, ring r)
167{
168  p_LmCheckPolyRing2(p, r);
169  pAssume2(v>0 && v <= r->N);
170  pAssume2(e>=0);
171  pAssume2((unsigned int) e<=r->bitmask);
172
173  // shift e to the left:
174  register int shift = r->VarOffset[v] >> 24;
175  unsigned long ee = ((unsigned long)e) << shift /*(r->VarOffset[v] >> 24)*/;
176  // find the bits in the exponent vector
177  register int offset = (r->VarOffset[v] & 0xffffff);
178  // clear the bits in the exponent vector:
179  p->exp[offset]  &= ~( r->bitmask << shift );
180  // insert e with |
181  p->exp[ offset ] |= ee;
182  return e;
183}
184
185// the following should be implemented more efficiently
186PINLINE2  int p_IncrExp(poly p, int v, ring r)
187{
188  p_LmCheckPolyRing2(p, r);
189  int e = p_GetExp(p,v,r);
190  e++;
191  return p_SetExp(p,v,e,r);
192}
193PINLINE2  int p_DecrExp(poly p, int v, ring r)
194{
195  p_LmCheckPolyRing2(p, r);
196  int e = p_GetExp(p,v,r);
197  pAssume2(e > 0);
198  e--;
199  return p_SetExp(p,v,e,r);
200}
201PINLINE2  int p_AddExp(poly p, int v, int ee, ring r)
202{
203  p_LmCheckPolyRing2(p, r);
204  int e = p_GetExp(p,v,r);
205  e += ee;
206  return p_SetExp(p,v,e,r);
207}
208PINLINE2  int p_SubExp(poly p, int v, int ee, ring r)
209{
210  p_LmCheckPolyRing2(p, r);
211  int e = p_GetExp(p,v,r);
212  pAssume2(e >= ee);
213  e -= ee;
214  return p_SetExp(p,v,e,r);
215}
216PINLINE2  int p_MultExp(poly p, int v, int ee, ring r)
217{
218  p_LmCheckPolyRing2(p, r);
219  int e = p_GetExp(p,v,r);
220  e *= ee;
221  return p_SetExp(p,v,e,r);
222}
223
224PINLINE2 int p_GetExpSum(poly p1, poly p2, int i, ring r)
225{
226  p_LmCheckPolyRing2(p1, r);
227  p_LmCheckPolyRing2(p2, r);
228  return p_GetExp(p1,i,r) + p_GetExp(p2,i,r);
229}
230PINLINE2 int p_GetExpDiff(poly p1, poly p2, int i, ring r)
231{
232  return p_GetExp(p1,i,r) - p_GetExp(p2,i,r);
233}
234
235
236/***************************************************************
237 *
238 * Allocation/Initalization/Deletion
239 *
240 ***************************************************************/
241PINLINE2 poly p_New(ring r, omBin bin)
242{
243  p_CheckRing2(r);
244  pAssume2(bin != NULL && r->PolyBin->sizeW == bin->sizeW);
245  poly p;
246  omTypeAllocBin(poly, p, bin);
247  p_SetRingOfLm(p, r);
248  return p;
249}
250
251PINLINE2 poly p_New(ring r)
252{
253  return p_New(r, r->PolyBin);
254}
255
256PINLINE2 void p_DeleteLm(poly *p, ring r)
257{
258  pIfThen2(*p != NULL, p_LmCheckPolyRing2(*p, r));
259  poly h = *p;
260  if (h != NULL)
261  {
262    n_Delete(&_pGetCoeff(h), r);
263    *p = _pNext(h);
264    omFreeBinAddr(h);
265  }
266}
267PINLINE2 void p_DeleteLm(poly p, ring r)
268{
269  pIfThen2(p != NULL, p_LmCheckPolyRing2(p, r));
270  if (p != NULL)
271  {
272    n_Delete(&_pGetCoeff(p), r);
273    omFreeBinAddr(p);
274  }
275}
276PINLINE2 void p_LmFree(poly p, ring r)
277{
278  p_LmCheckPolyRing2(p, r);
279  omFreeBinAddr(p);
280}
281PINLINE2 void p_LmFree(poly *p, ring r)
282{
283  p_LmCheckPolyRing2(*p, r);
284  poly h = *p;
285  *p = pNext(h);
286  omFreeBinAddr(h);
287}
288PINLINE2 poly p_LmFreeAndNext(poly p, ring r)
289{
290  p_LmCheckPolyRing2(p, r);
291  poly pnext = pNext(p);
292  omFreeBinAddr(p);
293  return pnext;
294}
295PINLINE2 void p_LmDelete(poly p, ring r)
296{
297  p_LmCheckPolyRing2(p, r);
298  n_Delete(&_pGetCoeff(p), r);
299  omFreeBinAddr(p);
300}
301PINLINE2 void p_LmDelete(poly *p, ring r)
302{
303  p_LmCheckPolyRing2(*p, r);
304  poly h = *p;
305  *p = pNext(h);
306  n_Delete(&pGetCoeff(h), r);
307  omFreeBinAddr(h);
308}
309PINLINE2 poly p_LmDeleteAndNext(poly p, ring r)
310{
311  p_LmCheckPolyRing2(p, r);
312  poly pnext = _pNext(p);
313  n_Delete(&_pGetCoeff(p), r);
314  omFreeBinAddr(p);
315  return pnext;
316}
317
318/***************************************************************
319 *
320 * Misc routines
321 *
322 ***************************************************************/
323PINLINE2 int p_Cmp(poly p1, poly p2, ring r)
324{
325  if (p2==NULL)
326    return 1;
327  if (p1==NULL)
328    return -1;
329  return p_LmCmp(p1,p2,r);
330}
331
332PINLINE2 unsigned long p_GetMaxExp(poly p, ring r)
333{
334  return p_GetMaxExp(p_GetMaxExpL(p, r), r);
335}
336
337PINLINE2 unsigned long
338p_GetMaxExp(const unsigned long l, const ring r, const int number_of_exps)
339{
340  unsigned long bitmask = r->bitmask;
341  unsigned long max = (l & bitmask);
342  unsigned long j = number_of_exps - 1;
343
344  if (j > 0)
345  {
346    unsigned long i = r->BitsPerExp;
347    long e;
348    loop
349    {
350      e = ((l >> i) & bitmask);
351      if ((unsigned long) e > max)
352        max = e;
353      j--;
354      if (j==0) break;
355      i += r->BitsPerExp;
356    }
357  }
358  return max;
359}
360
361PINLINE2 unsigned long p_GetMaxExp(const unsigned long l, const ring r)
362{
363  return p_GetMaxExp(l, r, r->ExpPerLong);
364}
365
366PINLINE2 unsigned long
367p_GetTotalDegree(const unsigned long l, const ring r, const int number_of_exps)
368{
369  const unsigned long bitmask = r->bitmask;
370  unsigned long sum = (l & bitmask);
371  unsigned long j = number_of_exps - 1;
372
373  if (j > 0)
374  {
375    unsigned long i = r->BitsPerExp;
376    loop
377    {
378      sum += ((l >> i) & bitmask);
379      j--;
380      if (j==0) break;
381      i += r->BitsPerExp;
382    }
383  }
384  return sum;
385}
386
387PINLINE2 unsigned long
388p_GetTotalDegree(const unsigned long l, const ring r)
389{
390  return p_GetTotalDegree(l, r, r->ExpPerLong);
391}
392
393/***************************************************************
394 *
395 * Dispatcher to r->p_Procs, they do the tests/checks
396 *
397 ***************************************************************/
398// returns a copy of p
399PINLINE2 poly p_Copy(poly p, const ring r)
400{
401#ifdef PDEBUG
402  poly pp= r->p_Procs->p_Copy(p, r);
403  p_Test(pp,r);
404  return pp;
405#else
406  return r->p_Procs->p_Copy(p, r);
407#endif
408}
409
410PINLINE2 poly p_Copy(poly p, const ring lmRing, const ring tailRing)
411{
412#ifndef PDEBUG
413  if (tailRing == lmRing)
414    return tailRing->p_Procs->p_Copy(p, tailRing);
415#endif
416  if (p != NULL)
417  {
418    poly pres = p_Head(p, lmRing);
419    pNext(pres) = tailRing->p_Procs->p_Copy(pNext(p), tailRing);
420    return pres;
421  }
422  else
423    return NULL;
424}
425
426// deletes *p, and sets *p to NULL
427PINLINE2 void p_Delete(poly *p, const ring r)
428{
429  r->p_Procs->p_Delete(p, r);
430}
431
432PINLINE2 void p_Delete(poly *p,  const ring lmRing, const ring tailRing)
433{
434#ifndef PDEBUG
435  if (tailRing == lmRing)
436  {
437    tailRing->p_Procs->p_Delete(p, tailRing);
438    return;
439  }
440#endif
441  if (*p != NULL)
442  {
443    if (pNext(*p) != NULL)
444      tailRing->p_Procs->p_Delete(&pNext(*p), tailRing);
445    p_LmDelete(p, lmRing);
446  }
447}
448
449PINLINE2 poly p_ShallowCopyDelete(poly p, const ring r, omBin bin)
450{
451  p_LmCheckPolyRing2(p, r);
452  pAssume2(r->PolyBin->sizeW == bin->sizeW);
453  return r->p_Procs->p_ShallowCopyDelete(p, r, bin);
454}
455
456// returns p+q, destroys p and q
457PINLINE2 poly p_Add_q(poly p, poly q, const ring r)
458{
459  int shorter;
460  return r->p_Procs->p_Add_q(p, q, shorter, r);
461}
462
463PINLINE2 poly p_Add_q(poly p, poly q, int &lp, int lq, const ring r)
464{
465  int shorter;
466  poly res = r->p_Procs->p_Add_q(p, q, shorter, r);
467  lp = (lp + lq) - shorter;
468  return res;
469}
470
471// returns p*n, destroys p
472PINLINE2 poly p_Mult_nn(poly p, number n, const ring r)
473{
474  if (n_IsOne(n, r))
475    return p;
476  else
477    return r->p_Procs->p_Mult_nn(p, n, r);
478}
479
480PINLINE2 poly p_Mult_nn(poly p, number n, const ring lmRing,
481                        const ring tailRing)
482{
483#ifndef PDEBUG
484  if (lmRing == tailRing)
485  {
486    return p_Mult_nn(p, n, tailRing);
487  }
488#endif
489  poly pnext = pNext(p);
490  pNext(p) = NULL;
491  p = lmRing->p_Procs->p_Mult_nn(p, n, lmRing);
492  pNext(p) = tailRing->p_Procs->p_Mult_nn(pnext, n, tailRing);
493  return p;
494}
495
496// returns p*n, does not destroy p
497PINLINE2 poly pp_Mult_nn(poly p, number n, const ring r)
498{
499  if (n_IsOne(n, r))
500    return p_Copy(p, r);
501  else
502    return r->p_Procs->pp_Mult_nn(p, n, r);
503}
504
505// returns Copy(p)*m, does neither destroy p nor m
506PINLINE2 poly pp_Mult_mm(poly p, poly m, const ring r)
507{
508  if (p_LmIsConstant(m, r))
509    return pp_Mult_nn(p, pGetCoeff(m), r);
510  else
511  {
512    poly last;
513    return r->p_Procs->pp_Mult_mm(p, m, r, last);
514  }
515}
516
517// returns p*m, destroys p, const: m
518PINLINE2 poly p_Mult_mm(poly p, poly m, const ring r)
519{
520  if (p_LmIsConstant(m, r))
521    return p_Mult_nn(p, pGetCoeff(m), r);
522  else
523    return r->p_Procs->p_Mult_mm(p, m, r);
524}
525
526// return p - m*Copy(q), destroys p; const: p,m
527PINLINE2 poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
528{
529#ifdef HAVE_PLURAL
530  if (rIsPluralRing(r))
531  {
532    int lp, lq;
533    poly spNoether;
534    return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
535  }
536#endif
537
538  int shorter;
539  poly last;
540
541  return r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last); // !!!
542}
543
544PINLINE2 poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
545                                 poly spNoether, const ring r)
546{
547#ifdef HAVE_PLURAL
548  if (rIsPluralRing(r))
549     return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
550#endif
551
552  int shorter;
553  poly last,res;
554  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, spNoether, r, last);
555  lp = (lp + lq) - shorter;
556  return res;
557}
558
559PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
560{
561  int shorter;
562  return r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
563}
564
565PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, int &lp, const poly m, const ring r)
566{
567  int shorter;
568  poly pp = r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
569  lp -= shorter;
570  return pp;
571}
572
573// returns -p, destroys p
574PINLINE2 poly p_Neg(poly p, const ring r)
575{
576  return r->p_Procs->p_Neg(p, r);
577}
578
579extern poly  _p_Mult_q(poly p, poly q, const int copy, const ring r);
580// returns p*q, destroys p and q
581PINLINE2 poly p_Mult_q(poly p, poly q, const ring r)
582{
583  if (p == NULL)
584  {
585    r->p_Procs->p_Delete(&q, r);
586    return NULL;
587  }
588  if (q == NULL)
589  {
590    r->p_Procs->p_Delete(&p, r);
591    return NULL;
592  }
593
594  if (pNext(p) == NULL)
595  {
596#ifdef HAVE_PLURAL
597    if (rIsPluralRing(r))
598      q = nc_mm_Mult_p(p, q, r);
599    else
600#endif /* HAVE_PLURAL */
601      q = r->p_Procs->p_Mult_mm(q, p, r);
602
603    r->p_Procs->p_Delete(&p, r);
604    return q;
605  }
606
607  if (pNext(q) == NULL)
608  {
609  // NEEDED
610#ifdef HAVE_PLURAL
611/*    if (rIsPluralRing(r))
612      p = gnc_p_Mult_mm(p, q, r); // ???
613    else*/
614#endif /* HAVE_PLURAL */
615      p = r->p_Procs->p_Mult_mm(p, q, r);
616
617    r->p_Procs->p_Delete(&q, r);
618    return p;
619  }
620#ifdef HAVE_PLURAL
621  if (rIsPluralRing(r))
622    return _nc_p_Mult_q(p, q, r);
623  else
624#endif
625  return _p_Mult_q(p, q, 0, r);
626}
627
628// returns p*q, does neither destroy p nor q
629PINLINE2 poly pp_Mult_qq(poly p, poly q, const ring r)
630{
631  poly last;
632  if (p == NULL || q == NULL) return NULL;
633
634  if (pNext(p) == NULL)
635  {
636#ifdef HAVE_PLURAL
637    if (rIsPluralRing(r))
638      return nc_mm_Mult_pp(p, q, r);
639#endif
640    return r->p_Procs->pp_Mult_mm(q, p, r, last);
641  }
642
643  if (pNext(q) == NULL)
644  {
645    return r->p_Procs->pp_Mult_mm(p, q, r, last);
646  }
647
648  poly qq = q;
649  if (p == q)
650    qq = p_Copy(q, r);
651
652  poly res;
653#ifdef HAVE_PLURAL
654  if (rIsPluralRing(r))
655    res = _nc_pp_Mult_qq(p, qq, r);
656  else
657#endif
658    res = _p_Mult_q(p, qq, 1, r);
659
660  if (qq != q)
661    p_Delete(&qq, r);
662  return res;
663}
664
665// returns p + m*q destroys p, const: q, m
666// this should be implemented more efficiently
667PINLINE2 poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
668                                const ring r)
669{
670#ifdef HAVE_PLURAL
671  if (rIsPluralRing(r))
672    return nc_p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
673#endif
674
675  poly res, last;
676  int shorter;
677  number n_old = pGetCoeff(m);
678  number n_neg = n_Copy(n_old, r);
679  n_neg = n_Neg(n_neg, r);
680  pSetCoeff0(m, n_neg);
681  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last);
682  lp = (lp + lq) - shorter;
683  pSetCoeff0(m, n_old);
684  n_Delete(&n_neg, r);
685  return res;
686}
687
688PINLINE2 poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
689{
690  int lp = 0, lq = 0;
691  return p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
692}
693
694PINLINE2 poly p_Merge_q(poly p, poly q, const ring r)
695{
696  return r->p_Procs->p_Merge_q(p, q, r);
697}
698
699PINLINE2 poly p_SortAdd(poly p, const ring r, BOOLEAN revert)
700{
701  if (revert) p = pReverse(p);
702  return sBucketSortAdd(p, r);
703}
704
705PINLINE2 poly p_SortMerge(poly p, const ring r, BOOLEAN revert)
706{
707  if (revert) p = pReverse(p);
708  return sBucketSortMerge(p, r);
709}
710
711/***************************************************************
712 *
713 * I/O
714 *
715 ***************************************************************/
716PINLINE2 char*     p_String(poly p, ring p_ring)
717{
718  return p_String(p, p_ring, p_ring);
719}
720PINLINE2 char*     p_String0(poly p, ring p_ring)
721{
722  return p_String0(p, p_ring, p_ring);
723}
724PINLINE2 void      p_Write(poly p, ring p_ring)
725{
726  p_Write(p, p_ring, p_ring);
727}
728PINLINE2 void      p_Write0(poly p, ring p_ring)
729{
730  p_Write0(p, p_ring, p_ring);
731}
732PINLINE2 void      p_wrp(poly p, ring p_ring)
733{
734  p_wrp(p, p_ring, p_ring);
735}
736#endif // !defined(NO_PINLINE2) || defined(POLYS_IMPL_CC)
737#endif // PINLINE2_H
Note: See TracBrowser for help on using the repository browser.