source: git/Singular/pcv.cc @ e4e36c

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