source: git/Singular/maps_ip.cc @ 9d68fd

jengelh-datetimespielwiese
Last change on this file since 9d68fd was 9d68fd, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: better tests and normalization for subst
  • Property mode set to 100644
File size: 10.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT - the mapping of polynomials to other rings
6*/
7#define TRANSEXT_PRIVATES
8
9#include "config.h"
10#include <kernel/mod2.h>
11#include <omalloc/omalloc.h>
12
13#include <coeffs/numbers.h>
14#include <coeffs/coeffs.h>
15
16#include <polys/monomials/ring.h>
17#include <polys/monomials/maps.h>
18#include <polys/matpol.h>
19#include <polys/prCopy.h>
20#include <polys/ext_fields/transext.h>
21
22//#include <libpolys/polys/ext_fields/longtrans.h>
23// #include <kernel/longalg.h>
24
25#include <kernel/febase.h>
26#include <kernel/kstd1.h>
27
28#include "maps_ip.h"
29#include "ipid.h"
30
31
32#include "lists.h"
33#include "tok.h"
34
35/* debug output: Tok2Cmdname in maApplyFetch*/
36#include "ipshell.h"
37
38// define this if you want to use the fast_map routine for mapping ideals
39//#define FAST_MAP
40
41#ifdef FAST_MAP
42#include <polys/monomials/maps.h>
43#endif
44
45
46/*2
47* maps the expression w to res,
48* switch what: MAP_CMD: use theMap for mapping, N for preimage ring
49*              //FETCH_CMD: use pOrdPoly for mapping
50*              IMAP_CMD: use perm for mapping, N for preimage ring
51*              default: map only poly-structures,
52*                       use perm and par_perm, N and P,
53*/
54BOOLEAN maApplyFetch(int what,map theMap,leftv res, leftv w, ring preimage_r,
55                     int *perm, int *par_perm, int P, nMapFunc nMap)
56{
57  int i;
58  int N = preimage_r->N;
59#if 0
60  Print("N=%d what=%s ",N,Tok2Cmdname(what));
61  if (perm!=NULL) for(i=1;i<=N;i++) Print("%d -> %d ",i,perm[i]);
62  PrintS("\n");
63  Print("P=%d ",P);
64  if (par_perm!=NULL) for(i=0;i<P;i++) Print("%d -> %d ",i,par_perm[i]);
65  PrintS("\n");
66#endif
67
68  void *data=w->Data();
69  res->rtyp = w->rtyp;
70  switch (w->rtyp)
71  {
72    case NUMBER_CMD:
73      if (P!=0)
74      {
75//        WerrorS("Sorry 'napPermNumber' was lost in the refactoring process (due to Frank): needs to be fixed");
76//        return TRUE;
77#if 1
78// poly n_PermNumber(const number z, const int *par_perm, const int OldPar, const ring src, const ring dst);
79        res->data= (void *) n_PermNumber((number)data, par_perm, P, preimage_r, currRing);
80#endif
81        res->rtyp=POLY_CMD;
82        if (nCoeff_is_Extension(currRing->cf))
83          res->data=(void *)p_MinPolyNormalize((poly)res->data, currRing);
84        pTest((poly) res->data);
85      }
86      else
87      {
88        assume( nMap != NULL );
89
90        number a = nMap((number)data, preimage_r->cf, currRing->cf);
91
92
93        if (nCoeff_is_Extension(currRing->cf))
94        {
95          n_Normalize(a, currRing->cf); // ???
96/*
97          number a = (number)res->data;
98          number one = nInit(1);
99          number product = nMult(a, one );
100          nDelete(&one);
101          nDelete(&a);
102          res->data=(void *)product;
103 */
104        }
105        #ifdef LDEBUG
106        n_Test(a, currRing->cf);
107        #endif
108        res->data=(void *)a;
109
110      }
111      break;
112    case POLY_CMD:
113    case VECTOR_CMD:
114      if ((what==FETCH_CMD)&& (preimage_r->cf==currRing->cf))
115        res->data=(void *)prCopyR( (poly)data, preimage_r, currRing);
116      else
117      if ((what==IMAP_CMD) || ((what==FETCH_CMD) /* && (nMap!=nCopy)*/))
118        res->data=(void *)p_PermPoly((poly)data,perm,preimage_r,currRing, nMap,par_perm,P);
119      else /*if (what==MAP_CMD)*/
120      {
121        p_Test((poly)data,preimage_r);
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  assume( nCoeff_is_transExt(currRing->cf) ); // nCoeff_is_Extension???
238  const ring R = currRing->cf->extRing;
239
240  ideal theMapI = idInit(rPar(currRing),1);
241  nMapFunc nMap = n_SetMap(R->cf, currRing->cf);
242
243  int i;
244  for(i = rPar(currRing);i>0;i--)
245  {
246    if (i != par)
247      theMapI->m[i-1]= p_NSet(n_Param(i, currRing), currRing);
248    else
249      theMapI->m[i-1] = p_Copy(image, currRing);
250    p_Test(theMapI->m[i-1],currRing);
251  }
252  //iiWriteMatrix((matrix)theMapI,"map:",1,currRing,0);
253
254  map theMap=(map)theMapI;
255  theMap->preimage=NULL;
256
257  leftv v=(leftv)omAllocBin(sleftv_bin);
258  sleftv tmpW;
259  poly res=NULL;
260
261  p_Normalize(p,currRing);
262  while (p!=NULL)
263  {
264    memset(v,0,sizeof(sleftv));
265
266    number d = n_GetDenom(p_GetCoeff(p, currRing), currRing);
267    assume( p_Test((poly)NUM(d), R) );
268
269    if ( n_IsOne (d, currRing->cf) )
270    {
271      n_Delete(&d, currRing); d = NULL;
272    }
273    else if (!p_IsConstant((poly)NUM(d), R))
274    {
275      WarnS("ignoring denominators of coefficients...");
276      n_Delete(&d, currRing); d = NULL;
277    }
278
279    number num = n_GetNumerator(p_GetCoeff(p, currRing), currRing);
280    assume( p_Test((poly)NUM(num), R) );
281
282    memset(&tmpW,0,sizeof(sleftv));
283    tmpW.rtyp = POLY_CMD;
284    tmpW.data = NUM (num); // a copy of this poly will be used
285
286    p_Normalize(NUM(num),R);
287    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,R,NULL,NULL,0,nMap))
288    {
289      WerrorS("map failed");
290      v->data=NULL;
291    }
292    n_Delete(&num, currRing);
293
294    //TODO check for memory leaks
295    poly pp = pHead(p);
296    //PrintS("map:");pWrite(pp);
297    if( d != NULL )
298    {
299      pSetCoeff(pp, n_Invers(d, currRing->cf));
300      n_Delete(&d, currRing); // d = NULL;
301    } else
302      pSetCoeff(pp, nInit(1));
303
304    //PrintS("->");pWrite((poly)(v->data));
305    poly ppp = pMult((poly)(v->data),pp);
306    //PrintS("->");pWrite(ppp);
307    res=pAdd(res,ppp);
308    pIter(p);
309  }
310  idDelete((ideal *)(&theMap));
311  omFreeBin((ADDRESS)v, sleftv_bin);
312  return res;
313}
314
315/*2
316* substitute the n-th parameter by the poly e in id
317* does not destroy id and e
318*/
319ideal  idSubstPar(ideal id, int n, poly e)
320{
321  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
322  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
323
324  res->rank = id->rank;
325  for(k--;k>=0;k--)
326  {
327    res->m[k]=pSubstPar(id->m[k],n,e);
328  }
329  return res;
330}
331
332/*2
333* substitutes the variable var (from 1..N) by image,
334* does not destroy p and  image
335*/
336poly pSubstPoly(poly p, int var, poly image)
337{
338  if (p==NULL) return NULL;
339#ifdef HAVE_PLURAL
340  if (rIsPluralRing(currRing))
341  {
342    return pSubst(pCopy(p),var,image);
343  }
344#endif
345  map theMap=(map)idMaxIdeal(1);
346  theMap->preimage=NULL;
347  pDelete(&(theMap->m[var-1]));
348  theMap->m[var-1]=pCopy(image);
349
350  poly res=NULL;
351#ifdef FAST_MAP
352  if (pGetComp(p)==0)
353  {
354    ideal src_id=idInit(1,1);
355    src_id->m[0]=p;
356    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
357    res=res_id->m[0];
358    res_id->m[0]=NULL; idDelete(&res_id);
359    src_id->m[0]=NULL; idDelete(&src_id);
360  }
361  else
362#endif
363  {
364    sleftv tmpW;
365    memset(&tmpW,0,sizeof(sleftv));
366    tmpW.rtyp=POLY_CMD;
367    tmpW.data=p;
368    leftv v=(leftv)omAlloc0Bin(sleftv_bin);
369    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
370                            n_SetMap(currRing->cf, currRing->cf)))
371    {
372      WerrorS("map failed");
373      v->data=NULL;
374    }
375    res=(poly)(v->data);
376    omFreeBin((ADDRESS)v, sleftv_bin);
377  }
378  idDelete((ideal *)(&theMap));
379  return res;
380}
381
382/*2
383* substitute the n-th variable by the poly e in id
384* does not destroy id and e
385*/
386ideal  idSubstPoly(ideal id, int n, poly e)
387{
388
389#ifdef HAVE_PLURAL
390  if (rIsPluralRing(currRing))
391  {
392    int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
393    ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
394    res->rank = id->rank;
395    for(k--;k>=0;k--)
396    {
397      res->m[k]=pSubst(pCopy(id->m[k]),n,e);
398    }
399    return res;
400  }
401#endif
402  map theMap=(map)idMaxIdeal(1);
403  theMap->preimage=NULL;
404  pDelete(&(theMap->m[n-1]));
405  theMap->m[n-1]=pCopy(e);
406
407  leftv v=(leftv)omAlloc0Bin(sleftv_bin);
408  sleftv tmpW;
409  memset(&tmpW,0,sizeof(sleftv));
410  tmpW.rtyp=IDEAL_CMD;
411  tmpW.data=id;
412  if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
413                          n_SetMap(currRing->cf, currRing->cf)))
414  {
415    WerrorS("map failed");
416    v->data=NULL;
417  }
418  ideal res=(ideal)(v->data);
419  idDelete((ideal *)(&theMap));
420  omFreeBin((ADDRESS)v, sleftv_bin);
421  return res;
422}
Note: See TracBrowser for help on using the repository browser.