source: git/Singular/maps_ip.cc @ 06c0b3

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