source: git/Singular/pcv.cc

spielwiese
Last change on this file was a3f0fea, checked in by Reimer Behrends <behrends@…>, 6 years ago
Modify variable declarions for pSingular.
  • Property mode set to 100644
File size: 8.9 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/*
5* ABSTRACT: conversion between polys and coef vectors
6*/
7
8
9
10
11#include "kernel/mod2.h"
12
13#ifdef HAVE_PCV
14
15#include "Singular/tok.h"
16#include "Singular/ipid.h"
17#include "coeffs/numbers.h"
18#include "kernel/polys.h"
19#include "kernel/ideals.h"
20#include "Singular/lists.h"
21#include "polys/matpol.h"
22#include "Singular/ipshell.h"
23#include "Singular/pcv.h"
24
25STATIC_VAR int pcvMaxDegree;
26STATIC_VAR int pcvTableSize;
27STATIC_VAR int pcvIndexSize;
28STATIC_VAR unsigned* pcvTable=NULL;
29STATIC_VAR unsigned** pcvIndex=NULL;
30
31lists pcvLAddL(lists l1,lists l2)
32{
33  lists l0=(lists)omAllocBin(slists_bin);
34  int i=l1->nr;
35  if(l1->nr<l2->nr) i=l2->nr;
36  l0->Init(i+1);
37  for(;i>=0;i--)
38  {
39    if(i<=l1->nr
40    &&(l1->m[i].rtyp==POLY_CMD||l1->m[i].rtyp==BUCKET_CMD||l1->m[i].rtyp==VECTOR_CMD))
41    {
42      if (l1->m[i].rtyp==BUCKET_CMD)
43      {
44        poly p=sBucketPeek((sBucket_pt)l1->m[i].data);
45        l0->m[i].rtyp=POLY_CMD;
46        l0->m[i].data=pCopy(p);
47      }
48      else
49      {
50        l0->m[i].rtyp=l1->m[i].rtyp;
51        l0->m[i].data=pCopy((poly)l1->m[i].data);
52      }
53      if(i<=l2->nr&&l2->m[i].rtyp==l0->m[i].rtyp)
54        l0->m[i].data=pAdd((poly)l0->m[i].data,pCopy((poly)l2->m[i].data));
55    }
56    else
57    if(i<=l2->nr
58    &&(l2->m[i].rtyp==POLY_CMD||l2->m[i].rtyp==BUCKET_CMD||l2->m[i].rtyp==VECTOR_CMD))
59    {
60      if (l2->m[i].rtyp==BUCKET_CMD)
61      {
62        poly p=sBucketPeek((sBucket_pt)l2->m[i].data);
63        l0->m[i].rtyp=POLY_CMD;
64        l0->m[i].data=pCopy(p);
65      }
66      else
67      {
68        l0->m[i].rtyp=l2->m[i].rtyp;
69        l0->m[i].data=pCopy((poly)l2->m[i].data);
70      }
71    }
72  }
73  return(l0);
74}
75
76lists pcvPMulL(poly p,lists l1)
77{
78  lists l0=(lists)omAllocBin(slists_bin);
79  l0->Init(l1->nr+1);
80  for(int i=l1->nr;i>=0;i--)
81  {
82    if(l1->m[i].rtyp==POLY_CMD)
83    {
84      l0->m[i].rtyp=POLY_CMD;
85      l0->m[i].data=ppMult_qq(p,(poly)l1->m[i].data);
86    }
87    else if(l1->m[i].rtyp==BUCKET_CMD)
88    {
89      l0->m[i].rtyp=POLY_CMD;
90      l0->m[i].data=ppMult_qq(p,sBucketPeek((sBucket_pt)l1->m[i].data));
91    }
92  }
93  return(l0);
94}
95
96BOOLEAN pcvLAddL(leftv res,leftv h)
97{
98  const short t[]={2,LIST_CMD,LIST_CMD};
99  if (iiCheckTypes(h,t,1))
100  {
101    lists l1=(lists)h->Data();
102    h=h->next;
103    lists l2=(lists)h->Data();
104    res->rtyp=LIST_CMD;
105    res->data=(void*)pcvLAddL(l1,l2);
106    return FALSE;
107  }
108  return TRUE;
109}
110
111BOOLEAN pcvPMulL(leftv res,leftv h)
112{
113  const short t1[]={2,POLY_CMD,LIST_CMD};
114  const short t2[]={2,BUCKET_CMD,LIST_CMD};
115  if (iiCheckTypes(h,t2,0)
116  ||iiCheckTypes(h,t1,1))
117  {
118    poly p=(poly)h->Data();
119    h=h->next;
120    lists l=(lists)h->Data();
121    res->rtyp=LIST_CMD;
122    res->data=(void*)pcvPMulL(p,l);
123    return FALSE;
124  }
125  return TRUE;
126}
127
128int pcvDeg(poly p)
129{
130  int d=0;
131  for(int i=currRing->N;i>=1;i--) d+=pGetExp(p,i);
132  return d;
133}
134
135int pcvMinDeg(poly p)
136{
137  if(!p) return -1;
138  int md=pcvDeg(p);
139  pIter(p);
140  while(p)
141  {
142    int d=pcvDeg(p);
143    if(d<md) md=d;
144    pIter(p);
145  }
146  return md;
147}
148
149int pcvMinDeg(matrix m)
150{
151  int i,j,d;
152  int md=-1;
153  for(i=1;i<=MATROWS(m);i++)
154  {
155    for(j=1;j<=MATCOLS(m);j++)
156    {
157      d=pcvMinDeg(MATELEM(m,i,j));
158      if((d>=0&&md>d)||md==-1) md=d;
159    }
160  }
161  return(md);
162}
163
164BOOLEAN pcvMinDeg(leftv res,leftv h)
165{
166  if(h)
167  {
168    if(h->Typ()==BUCKET_CMD)
169    {
170      res->rtyp=INT_CMD;
171      res->data=(void*)(long)pcvMinDeg(sBucketPeek((sBucket_pt)h->Data()));
172      return FALSE;
173    }
174    else
175    if(h->Typ()==POLY_CMD)
176    {
177      res->rtyp=INT_CMD;
178      res->data=(void*)(long)pcvMinDeg((poly)h->Data());
179      return FALSE;
180    }
181    else
182    if(h->Typ()==MATRIX_CMD)
183    {
184      res->rtyp=INT_CMD;
185      res->data=(void*)(long)pcvMinDeg((matrix)h->Data());
186      return FALSE;
187    }
188  }
189  WerrorS("<poly> expected");
190  return TRUE;
191}
192
193void pcvInit(int d)
194{
195  if(d<0) d=1;
196  pcvMaxDegree=d+1;
197  pcvTableSize=currRing->N*pcvMaxDegree*sizeof(unsigned);
198  pcvTable=(unsigned*)omAlloc0(pcvTableSize);
199  pcvIndexSize=currRing->N*sizeof(unsigned*);
200  pcvIndex=(unsigned**)omAlloc(pcvIndexSize);
201  for(int i=0;i<currRing->N;i++)
202    pcvIndex[i]=pcvTable+i*pcvMaxDegree;
203  for(int i=0;i<pcvMaxDegree;i++)
204    pcvIndex[0][i]=i;
205  unsigned k,l;
206  for(int i=1;i<currRing->N;i++)
207  {
208    k=0;
209    for(int j=0;j<pcvMaxDegree;j++)
210    {
211      l=pcvIndex[i-1][j];
212      if(l>unsigned(~0)-k)
213      {
214        j=pcvMaxDegree;
215        i=currRing->N;
216        WerrorS("unsigned overflow");
217      }
218      else pcvIndex[i][j]=k+=l;
219    }
220  }
221}
222
223void pcvClean()
224{
225  if(pcvTable)
226  {
227    omFreeSize(pcvTable,pcvTableSize);
228    pcvTable=NULL;
229  }
230  if(pcvIndex)
231  {
232    omFreeSize(pcvIndex,pcvIndexSize);
233    pcvIndex=NULL;
234  }
235}
236
237int pcvM2N(poly m)
238{
239  unsigned n=0,dn,d=0;
240  for(int i=0;i<currRing->N;i++)
241  {
242    d+=pGetExp(m,i+1);
243    dn=pcvIndex[i][d];
244    if(dn>MAX_INT_VAL-n)
245    {
246      i=currRing->N;
247      WerrorS("component overflow");
248    }
249    else n+=dn;
250  }
251  return n+1;
252}
253
254poly pcvN2M(int n)
255{
256  n--;
257  poly m=pOne();
258  int i,j=0,k;
259  for(i=currRing->N-1;i>=0;i--)
260  {
261    k=j;
262    for(j=0; (j<pcvMaxDegree) && (pcvIndex[i][j]<=(unsigned)n); j++);
263    j--;
264    n-=pcvIndex[i][j];
265    if(i<currRing->N-1) pSetExp(m,i+2,k-j);
266  }
267  if(n==0)
268  {
269    pSetExp(m,1,j);
270    pSetm(m);
271    return m;
272  }
273  else
274  {
275    pLmDelete(&m);
276    return NULL;
277  }
278}
279
280poly pcvP2CV(poly p,int d0,int d1)
281{
282  poly cv=NULL;
283  while(p)
284  {
285    int d=pcvDeg(p);
286    if(d0<=d&&d<d1)
287    {
288      poly c=pNSet(nCopy(pGetCoeff(p)));
289      pSetComp(c,pcvM2N(p));
290      cv=pAdd(cv,c);
291    }
292    pIter(p);
293  }
294  return cv;
295}
296
297poly pcvCV2P(poly cv,int d0,int d1)
298{
299  poly p=NULL;
300  while(cv)
301  {
302    poly m=pcvN2M(pGetComp(cv));
303    if(m)
304    {
305      int d=pcvDeg(m);
306      if(d0<=d&&d<d1)
307      {
308        pSetCoeff(m,nCopy(pGetCoeff(cv)));
309        p=pAdd(p,m);
310      }
311    }
312    pIter(cv);
313  }
314  return p;
315}
316
317lists pcvP2CV(lists pl,int d0,int d1)
318{
319  lists cvl=(lists)omAllocBin(slists_bin);
320  cvl->Init(pl->nr+1);
321  pcvInit(d1);
322  for(int i=pl->nr;i>=0;i--)
323  {
324    if(pl->m[i].rtyp==POLY_CMD)
325    {
326      cvl->m[i].rtyp=VECTOR_CMD;
327      cvl->m[i].data=pcvP2CV((poly)pl->m[i].data,d0,d1);
328    }
329    else if(pl->m[i].rtyp==BUCKET_CMD)
330    {
331      cvl->m[i].rtyp=VECTOR_CMD;
332      cvl->m[i].data=pcvP2CV(sBucketPeek((sBucket_pt)pl->m[i].data),d0,d1);
333    }
334  }
335  pcvClean();
336  return cvl;
337}
338
339lists pcvCV2P(lists cvl,int d0,int d1)
340{
341  lists pl=(lists)omAllocBin(slists_bin);
342  pl->Init(cvl->nr+1);
343  pcvInit(d1);
344  for(int i=cvl->nr;i>=0;i--)
345  {
346    if(cvl->m[i].rtyp==VECTOR_CMD)
347    {
348      pl->m[i].rtyp=POLY_CMD;
349      pl->m[i].data=pcvCV2P((poly)cvl->m[i].data,d0,d1);
350    }
351  }
352  pcvClean();
353  return pl;
354}
355
356BOOLEAN pcvP2CV(leftv res,leftv h)
357{
358  if(currRing)
359  {
360    const short t[]={3,LIST_CMD,INT_CMD,INT_CMD};
361    if (iiCheckTypes(h,t,1))
362    {
363      lists p=(lists)h->Data();
364      h=h->next;
365      int d0=(int)(long)h->Data();
366      h=h->next;
367      int d1=(int)(long)h->Data();
368      res->rtyp=LIST_CMD;
369      res->data=(void*)pcvP2CV(p,d0,d1);
370      return FALSE;
371    }
372    return TRUE;
373  }
374  WerrorS("no ring active");
375  return TRUE;
376}
377
378BOOLEAN pcvCV2P(leftv res,leftv h)
379{
380  if(currRing)
381  {
382    const short t[]={3,LIST_CMD,INT_CMD,INT_CMD};
383    if (iiCheckTypes(h,t,1))
384    {
385      lists pl=(lists)h->Data();
386      h=h->next;
387      int d0=(int)(long)h->Data();
388      h=h->next;
389      int d1=(int)(long)h->Data();
390      res->rtyp=LIST_CMD;
391      res->data=(void*)pcvCV2P(pl,d0,d1);
392      return FALSE;
393    }
394    return TRUE;
395  }
396  WerrorS("no ring active");
397  return TRUE;
398}
399
400int pcvDim(int d0,int d1)
401{
402  if(d0<0) d0=0;
403  if(d1<0) d1=0;
404  pcvInit(d1);
405  int d=pcvIndex[currRing->N-1][d1]-pcvIndex[currRing->N-1][d0];
406  pcvClean();
407  return d;
408}
409
410BOOLEAN pcvDim(leftv res,leftv h)
411{
412  if(currRing)
413  {
414    const short t[]={2,INT_CMD,INT_CMD};
415    if (iiCheckTypes(h,t,1))
416    {
417      int d0=(int)(long)h->Data();
418      h=h->next;
419      int d1=(int)(long)h->Data();
420      res->rtyp=INT_CMD;
421      res->data=(void*)(long)pcvDim(d0,d1);
422      return FALSE;
423    }
424    return TRUE;
425  }
426  WerrorS("no ring active");
427  return TRUE;
428}
429
430int pcvBasis(lists b,int i,poly m,int d,int n)
431{
432  if(n<currRing->N)
433  {
434    for(int k=0,l=d;k<=l;k++,d--)
435    {
436      pSetExp(m,n,k);
437      i=pcvBasis(b,i,m,d,n+1);
438    }
439  }
440  else
441  {
442    pSetExp(m,n,d);
443    pSetm(m);
444    b->m[i].rtyp=POLY_CMD;
445    b->m[i++].data=pCopy(m);
446  }
447  return i;
448}
449
450lists pcvBasis(int d0,int d1)
451{
452  if(d0<0) d0=0;
453  if(d1<0) d1=0;
454  lists b=(lists)omAllocBin(slists_bin);
455  b->Init(pcvDim(d0,d1));
456  poly m=pOne();
457  for(int d=d0,i=0;d<d1;d++)
458    i=pcvBasis(b,i,m,d,1);
459  pLmDelete(&m);
460  return b;
461}
462
463BOOLEAN pcvBasis(leftv res,leftv h)
464{
465  if(currRing)
466  {
467    const short t[]={2,INT_CMD,INT_CMD};
468    if (iiCheckTypes(h,t,1))
469    {
470      int d0=(int)(long)h->Data();
471      h=h->next;
472      int d1=(int)(long)h->Data();
473      res->rtyp=LIST_CMD;
474      res->data=(void*)pcvBasis(d0,d1);
475      return FALSE;
476    }
477    return TRUE;
478  }
479  WerrorS("no ring active");
480  return TRUE;
481}
482
483#endif /* HAVE_PCV */
Note: See TracBrowser for help on using the repository browser.