source: git/Singular/maps_ip.cc @ 93fdae

spielwiese
Last change on this file since 93fdae was ca0d043, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix tr. 317 git-svn-id: file:///usr/local/Singular/svn/trunk@13949 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT - the mapping of polynomials to other rings
7*/
8
9#include <kernel/mod2.h>
10#include <Singular/tok.h>
11#include <kernel/febase.h>
12#include <kernel/polys.h>
13#include <kernel/numbers.h>
14#include <kernel/ring.h>
15#include <kernel/ideals.h>
16#include <kernel/matpol.h>
17#include <omalloc/omalloc.h>
18#include <kernel/kstd1.h>
19#include <kernel/longalg.h>
20#include <Singular/lists.h>
21#include <kernel/maps.h>
22#include <Singular/maps_ip.h>
23#include <kernel/prCopy.h>
24
25// define this if you want to use the fast_map routine for mapping ideals
26//#define FAST_MAP
27
28#ifdef FAST_MAP
29#include <kernel/maps.h>
30#endif
31
32
33/*2
34* maps the expression w to res,
35* switch what: MAP_CMD: use theMap for mapping, N for preimage ring
36*              //FETCH_CMD: use pOrdPoly for mapping
37*              IMAP_CMD: use perm for mapping, N for preimage ring
38*              default: map only poly-structures,
39*                       use perm and par_perm, N and P,
40*/
41BOOLEAN maApplyFetch(int what,map theMap,leftv res, leftv w, ring preimage_r,
42                     int *perm, int *par_perm, int P, nMapFunc nMap)
43{
44  int i;
45  int N = preimage_r->N;
46  //Print("N=%d what=%s ",N,Tok2Cmdname(what));
47  //if (perm!=NULL) for(i=1;i<=N;i++) Print("%d -> %d ",i,perm[i]);
48  //PrintS("\n");
49  //Print("P=%d ",P);
50  //if (par_perm!=NULL) for(i=0;i<P;i++) Print("%d -> %d ",i,par_perm[i]);
51  //PrintS("\n");
52  void *data=w->Data();
53  res->rtyp = w->rtyp;
54  switch (w->rtyp)
55  {
56    case NUMBER_CMD:
57      if (P!=0)
58      {
59        res->data=(void *)napPermNumber((number)data,par_perm,P, preimage_r);
60        res->rtyp=POLY_CMD;
61        if (currRing->minpoly!=NULL)
62          res->data=(void *)pMinPolyNormalize((poly)res->data);
63        pTest((poly) res->data);
64      }
65      else
66      {
67        res->data=(void *)nMap((number)data);
68        if (currRing->minpoly!=NULL)
69        {
70          number a=(number)res->data;
71          number one=nInit(1);
72          number product = nMult(a,one );
73          nDelete(&one);
74          nDelete(&a);
75          res->data=(void *)product;
76        }
77        #ifdef LDEBUG
78        nTest((number) res->data);
79        #endif
80      }
81      break;
82    case POLY_CMD:
83    case VECTOR_CMD:
84      if ((what==FETCH_CMD)&& (nMap==nCopy))
85        res->data=(void *)prCopyR( (poly)data, preimage_r);
86      else
87      if ((what==IMAP_CMD) || ((what==FETCH_CMD) /* && (nMap!=nCopy)*/))
88        res->data=(void *)pPermPoly((poly)data,perm,preimage_r,nMap,par_perm,P);
89      else /*if (what==MAP_CMD)*/
90      {
91        matrix s=mpNew(N,maMaxDeg_P((poly)data, preimage_r));
92        res->data=(void *)maEval(theMap,(poly)data,preimage_r,nMap,s);
93        idDelete((ideal *)&s);
94      }
95      if (currRing->minpoly!=NULL)
96        res->data=(void *)pMinPolyNormalize((poly)res->data);
97      pTest((poly)res->data);
98      break;
99    case MODUL_CMD:
100    case MATRIX_CMD:
101    case IDEAL_CMD:
102    case MAP_CMD:
103    {
104      int C=((matrix)data)->cols();
105      int R;
106      if (w->rtyp==MAP_CMD) R=1;
107      else R=((matrix)data)->rows();
108      matrix m=mpNew(R,C);
109      char *tmpR=NULL;
110      if(w->rtyp==MAP_CMD)
111      {
112        tmpR=((map)data)->preimage;
113        ((matrix)data)->rank=((matrix)data)->rows();
114      }
115      if ((what==FETCH_CMD)&& (nMap==nCopy))
116      {
117        for (i=R*C-1;i>=0;i--)
118        {
119          m->m[i]=prCopyR(((ideal)data)->m[i], preimage_r);
120          pTest(m->m[i]);
121        }
122      }
123      else
124      if ((what==IMAP_CMD) || ((what==FETCH_CMD) /* && (nMap!=nCopy)*/))
125      {
126        for (i=R*C-1;i>=0;i--)
127        {
128          m->m[i]=pPermPoly(((ideal)data)->m[i],perm,preimage_r,nMap,par_perm,P);
129          pTest(m->m[i]);
130        }
131      }
132      else /* if(what==MAP_CMD) */
133      {
134        matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r));
135        for (i=R*C-1;i>=0;i--)
136        {
137          m->m[i]=maEval(theMap,((ideal)data)->m[i],preimage_r,nMap,s);
138          pTest(m->m[i]);
139        }
140        idDelete((ideal *)&s);
141      }
142      if (currRing->minpoly!=NULL)
143      {
144        for (i=R*C-1;i>=0;i--)
145        {
146          m->m[i]=pMinPolyNormalize(m->m[i]);
147          pTest(m->m[i]);
148        }
149      }
150      if(w->rtyp==MAP_CMD)
151      {
152        ((map)data)->preimage=tmpR;
153        ((map)m)->preimage=omStrDup(tmpR);
154      }
155      else
156      {
157        m->rank=((matrix)data)->rank;
158      }
159      res->data=(char *)m;
160      idTest((ideal) m);
161      break;
162    }
163
164    case LIST_CMD:
165    {
166      lists l=(lists)data;
167      lists ml=(lists)omAllocBin(slists_bin);
168      ml->Init(l->nr+1);
169      for(i=0;i<=l->nr;i++)
170      {
171        if (((l->m[i].rtyp>BEGIN_RING)&&(l->m[i].rtyp<END_RING))
172        ||(l->m[i].rtyp==LIST_CMD))
173        {
174          if (maApplyFetch(what,theMap,&ml->m[i],&l->m[i],
175                           preimage_r,perm,par_perm,P,nMap))
176          {
177            ml->Clean();
178            omFreeBin((ADDRESS)ml, slists_bin);
179            res->rtyp=0;
180            return TRUE;
181          }
182        }
183        else
184        {
185          ml->m[i].Copy(&l->m[i]);
186        }
187      }
188      res->data=(char *)ml;
189      break;
190    }
191    default:
192    {
193      return TRUE;
194    }
195  }
196  return FALSE;
197}
198
199/*2
200* substitutes the parameter par (from 1..N) by image,
201* does not destroy p and  image
202*/
203poly pSubstPar(poly p, int par, poly image)
204{
205  ideal theMapI=idInit(rPar(currRing),1);
206  nMapFunc nMap=nSetMap(currRing->algring);
207
208  int i;
209  poly pp;
210  for(i=rPar(currRing);i>0;i--)
211  {
212    if (i!=par)
213    {
214      pp=theMapI->m[i-1]=pOne();
215      lnumber n=(lnumber)pGetCoeff(pp);
216      p_SetExp(n->z,i,1,currRing->algring);
217      p_Setm(n->z,currRing->algring);
218    }
219    else
220      theMapI->m[i-1]=pCopy(image);
221  }
222
223  map theMap=(map)theMapI;
224  theMap->preimage=NULL;
225
226  leftv v=(leftv)omAllocBin(sleftv_bin);
227  sleftv tmpW;
228  poly res=NULL;
229  while (p!=NULL)
230  {
231    memset(&tmpW,0,sizeof(sleftv));
232    memset(v,0,sizeof(sleftv));
233    tmpW.rtyp=POLY_CMD;
234    lnumber n=(lnumber)pGetCoeff(p);
235    tmpW.data=n->z;
236    if (n->n!=NULL) WarnS("ignoring denominators of coefficients...");
237    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing->algring,NULL,NULL,0,nMap))
238    {
239      WerrorS("map failed");
240      v->data=NULL;
241    }
242    pp=pHead(p);
243    //PrintS("map:");pWrite(pp);
244    pSetCoeff(pp,nInit(1));
245    //PrintS("->");pWrite((poly)(v->data));
246    poly ppp=pMult((poly)(v->data),pp);
247    //PrintS("->");pWrite(ppp);
248    res=pAdd(res,ppp);
249    pIter(p);
250  }
251  idDelete((ideal *)(&theMap));
252  omFreeBin((ADDRESS)v, sleftv_bin);
253  return res;
254}
255
256/*2
257* substitute the n-th parameter by the poly e in id
258* does not destroy id and e
259*/
260ideal  idSubstPar(ideal id, int n, poly e)
261{
262  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
263  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
264
265  res->rank = id->rank;
266  for(k--;k>=0;k--)
267  {
268    res->m[k]=pSubstPar(id->m[k],n,e);
269  }
270  return res;
271}
272
273/*2
274* substitutes the variable var (from 1..N) by image,
275* does not destroy p and  image
276*/
277poly pSubstPoly(poly p, int var, poly image)
278{
279  if (p==NULL) return NULL;
280#ifdef HAVE_PLURAL
281  if (rIsPluralRing(currRing))
282  {
283    return pSubst(pCopy(p),var,image);
284  }
285#endif
286  map theMap=(map)idMaxIdeal(1);
287  theMap->preimage=NULL;
288  pDelete(&(theMap->m[var-1]));
289  theMap->m[var-1]=pCopy(image);
290
291  poly res=NULL;
292#ifdef FAST_MAP
293  if (pGetComp(p)==0)
294  {
295    ideal src_id=idInit(1,1);
296    src_id->m[0]=p;
297    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
298    res=res_id->m[0];
299    res_id->m[0]=NULL; idDelete(&res_id);
300    src_id->m[0]=NULL; idDelete(&src_id);
301  }
302  else
303#endif
304  {
305    sleftv tmpW;
306    memset(&tmpW,0,sizeof(sleftv));
307    tmpW.rtyp=POLY_CMD;
308    tmpW.data=p;
309    leftv v=(leftv)omAlloc0Bin(sleftv_bin);
310    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,nCopy))
311    {
312      WerrorS("map failed");
313      v->data=NULL;
314    }
315    res=(poly)(v->data);
316    omFreeBin((ADDRESS)v, sleftv_bin);
317  }
318  idDelete((ideal *)(&theMap));
319  return res;
320}
321
322/*2
323* substitute the n-th variable by the poly e in id
324* does not destroy id and e
325*/
326ideal  idSubstPoly(ideal id, int n, poly e)
327{
328
329#ifdef HAVE_PLURAL
330  if (rIsPluralRing(currRing))
331  {
332    int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
333    ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
334    res->rank = id->rank;
335    for(k--;k>=0;k--)
336    {
337      res->m[k]=pSubst(pCopy(id->m[k]),n,e);
338    }
339    return res;
340  }
341#endif
342  map theMap=(map)idMaxIdeal(1);
343  theMap->preimage=NULL;
344  pDelete(&(theMap->m[n-1]));
345  theMap->m[n-1]=pCopy(e);
346
347  leftv v=(leftv)omAlloc0Bin(sleftv_bin);
348  sleftv tmpW;
349  memset(&tmpW,0,sizeof(sleftv));
350  tmpW.rtyp=IDEAL_CMD;
351  tmpW.data=id;
352  if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,nCopy))
353  {
354    WerrorS("map failed");
355    v->data=NULL;
356  }
357  ideal res=(ideal)(v->data);
358  idDelete((ideal *)(&theMap));
359  omFreeBin((ADDRESS)v, sleftv_bin);
360  return res;
361}
Note: See TracBrowser for help on using the repository browser.