source: git/Singular/maps_ip.cc @ 4470e25

fieker-DuValspielwiese
Last change on this file since 4470e25 was 4470e25, checked in by Hans Schoenemann <hannes@…>, 9 years ago
Revert "fix: mapping needs normalize only for algExt" This reverts commit 7c852247b0258944ebdc6c5e200ebd31b23ff5ff.
  • Property mode set to 100644
File size: 11.7 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
10
11
12#include <kernel/mod2.h>
13#include <omalloc/omalloc.h>
14
15#include <coeffs/numbers.h>
16#include <coeffs/coeffs.h>
17
18#include <polys/monomials/ring.h>
19#include <polys/monomials/maps.h>
20#include <polys/matpol.h>
21#include <polys/prCopy.h>
22#include <polys/ext_fields/transext.h>
23
24//#include <polys/ext_fields/longtrans.h>
25// #include <kernel/longalg.h>
26
27#include <kernel/GBEngine/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/*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  const ring R = currRing->cf->extRing;
238  ideal theMapI = idInit(rPar(currRing),1);
239  nMapFunc nMap = n_SetMap(R->cf, currRing->cf);
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    p_Test(theMapI->m[i-1],currRing);
248  }
249  //iiWriteMatrix((matrix)theMapI,"map:",1,currRing,0);
250
251  map theMap=(map)theMapI;
252  theMap->preimage=NULL;
253
254  leftv v=(leftv)omAllocBin(sleftv_bin);
255  sleftv tmpW;
256  poly res=NULL;
257
258  p_Normalize(p,currRing);
259  if (currRing->cf->rep==n_rep_rat_fct )
260  {
261    while (p!=NULL)
262    {
263      memset(v,0,sizeof(sleftv));
264
265      number d = n_GetDenom(p_GetCoeff(p, currRing), currRing);
266      p_Test((poly)NUM((fraction)d), R);
267
268      if ( n_IsOne (d, currRing->cf) )
269      {
270        n_Delete(&d, currRing); d = NULL;
271      }
272       else if (!p_IsConstant((poly)NUM((fraction)d), R))
273      {
274        WarnS("ignoring denominators of coefficients...");
275        n_Delete(&d, currRing); d = NULL;
276      }
277
278      number num = n_GetNumerator(p_GetCoeff(p, currRing), currRing);
279      memset(&tmpW,0,sizeof(sleftv));
280      tmpW.rtyp = POLY_CMD;
281      p_Test((poly)NUM((fraction)num), R);
282
283      tmpW.data = NUM ((fraction)num); // a copy of this poly will be used
284
285      p_Normalize(NUM((fraction)num),R);
286      if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,R,NULL,NULL,0,nMap))
287      {
288        WerrorS("map failed");
289        v->data=NULL;
290      }
291      n_Delete(&num, currRing);
292      //TODO check for memory leaks
293      poly pp = pHead(p);
294      //PrintS("map:");pWrite(pp);
295      if( d != NULL )
296      {
297        pSetCoeff(pp, n_Invers(d, currRing->cf));
298        n_Delete(&d, currRing); // d = NULL;
299      }
300      else
301        pSetCoeff(pp, nInit(1));
302
303      //PrintS("->");pWrite((poly)(v->data));
304      poly ppp = pMult((poly)(v->data),pp);
305      //PrintS("->");pWrite(ppp);
306      res=pAdd(res,ppp);
307      pIter(p);
308    }
309  }
310  else if (currRing->cf->rep==n_rep_poly )
311  {
312    while (p!=NULL)
313    {
314      memset(v,0,sizeof(sleftv));
315
316      number num = n_GetNumerator(p_GetCoeff(p, currRing), currRing);
317      memset(&tmpW,0,sizeof(sleftv));
318      tmpW.rtyp = POLY_CMD;
319      p_Test((poly)num, R);
320
321
322      p_Normalize((poly)num,R);
323      if (num==NULL) num=(number)R->qideal->m[0];
324      tmpW.data = num; // a copy of this poly will be used
325      if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,R,NULL,NULL,0,nMap))
326      {
327        WerrorS("map failed");
328        v->data=NULL;
329      }
330      if (num!=(number)R->qideal->m[0]) n_Delete(&num, currRing);
331      //TODO check for memory leaks
332      poly pp = pHead(p);
333      //PrintS("map:");pWrite(pp);
334      pSetCoeff(pp,n_Init(1,currRing));
335      //PrintS("cf->");pWrite((poly)(v->data));
336      poly ppp = pMult((poly)(v->data),pp);
337      //PrintS("->");pWrite(ppp);
338      res=pAdd(res,ppp);
339      pIter(p);
340    }
341  }
342  else
343  {
344    WerrorS("cannot apply subst for these coeffcients");
345  }
346  idDelete((ideal *)(&theMap));
347  omFreeBin((ADDRESS)v, sleftv_bin);
348  return res;
349}
350
351/*2
352* substitute the n-th parameter by the poly e in id
353* does not destroy id and e
354*/
355ideal  idSubstPar(ideal id, int n, poly e)
356{
357  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
358  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
359
360  res->rank = id->rank;
361  for(k--;k>=0;k--)
362  {
363    res->m[k]=pSubstPar(id->m[k],n,e);
364  }
365  return res;
366}
367
368/*2
369* substitutes the variable var (from 1..N) by image,
370* does not destroy p and  image
371*/
372poly pSubstPoly(poly p, int var, poly image)
373{
374  if (p==NULL) return NULL;
375#ifdef HAVE_PLURAL
376  if (rIsPluralRing(currRing))
377  {
378    return pSubst(pCopy(p),var,image);
379  }
380#endif
381  map theMap=(map)idMaxIdeal(1);
382  theMap->preimage=NULL;
383  pDelete(&(theMap->m[var-1]));
384  theMap->m[var-1]=pCopy(image);
385
386  poly res=NULL;
387#ifdef FAST_MAP
388  if (pGetComp(p)==0)
389  {
390    ideal src_id=idInit(1,1);
391    src_id->m[0]=p;
392
393    char *tmp = theMap->preimage;
394    theMap->preimagei=(char*)1L; // map gets 1 as its rank (as an ideal)
395    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
396    theMap->preimage=tmp; // map gets its preimage back
397
398    res=res_id->m[0];
399    res_id->m[0]=NULL; idDelete(&res_id);
400    src_id->m[0]=NULL; idDelete(&src_id);
401  }
402  else
403#endif
404  {
405    sleftv tmpW;
406    memset(&tmpW,0,sizeof(sleftv));
407    tmpW.rtyp=POLY_CMD;
408    tmpW.data=p;
409    leftv v=(leftv)omAlloc0Bin(sleftv_bin);
410    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
411                            n_SetMap(currRing->cf, currRing->cf)))
412    {
413      WerrorS("map failed");
414      v->data=NULL;
415    }
416    res=(poly)(v->data);
417    omFreeBin((ADDRESS)v, sleftv_bin);
418  }
419  idDelete((ideal *)(&theMap));
420  return res;
421}
422
423/*2
424* substitute the n-th variable by the poly e in id
425* does not destroy id and e
426*/
427ideal  idSubstPoly(ideal id, int n, poly e)
428{
429
430#ifdef HAVE_PLURAL
431  if (rIsPluralRing(currRing))
432  {
433    int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
434    ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
435    res->rank = id->rank;
436    for(k--;k>=0;k--)
437    {
438      res->m[k]=pSubst(pCopy(id->m[k]),n,e);
439    }
440    return res;
441  }
442#endif
443  map theMap=(map)idMaxIdeal(1);
444  theMap->preimage=NULL;
445  pDelete(&(theMap->m[n-1]));
446  theMap->m[n-1]=pCopy(e);
447
448  leftv v=(leftv)omAlloc0Bin(sleftv_bin);
449  sleftv tmpW;
450  memset(&tmpW,0,sizeof(sleftv));
451  tmpW.rtyp=IDEAL_CMD;
452  tmpW.data=id;
453  if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
454                          n_SetMap(currRing->cf, currRing->cf)))
455  {
456    WerrorS("map failed");
457    v->data=NULL;
458  }
459  ideal res=(ideal)(v->data);
460  idDelete((ideal *)(&theMap));
461  omFreeBin((ADDRESS)v, sleftv_bin);
462  return res;
463}
Note: See TracBrowser for help on using the repository browser.