source: git/Singular/maps_ip.cc @ eb55f8a

spielwiese
Last change on this file since eb55f8a was eb55f8a, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
There should be no *Test in assume call (and no assume in *Test definition)
  • Property mode set to 100644
File size: 11.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
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 <libpolys/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
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        p_Test((poly)data,preimage_r);
123        matrix s=mpNew(N,maMaxDeg_P((poly)data, preimage_r));
124        res->data=(void *)maEval(theMap, (poly)data, preimage_r, nMap, (ideal)s, currRing);
125        idDelete((ideal *)&s);
126      }
127      if (nCoeff_is_Extension(currRing->cf))
128        res->data=(void *)p_MinPolyNormalize((poly)res->data, currRing);
129      pTest((poly)res->data);
130      break;
131    case MODUL_CMD:
132    case MATRIX_CMD:
133    case IDEAL_CMD:
134    case MAP_CMD:
135    {
136      int C=((matrix)data)->cols();
137      int R;
138      if (w->rtyp==MAP_CMD) R=1;
139      else R=((matrix)data)->rows();
140      matrix m=mpNew(R,C);
141      char *tmpR=NULL;
142      if(w->rtyp==MAP_CMD)
143      {
144        tmpR=((map)data)->preimage;
145        ((matrix)data)->rank=((matrix)data)->rows();
146      }
147      if ((what==FETCH_CMD)&& (preimage_r->cf == currRing->cf))
148      {
149        for (i=R*C-1;i>=0;i--)
150        {
151          m->m[i]=prCopyR(((ideal)data)->m[i], preimage_r, currRing);
152          pTest(m->m[i]);
153        }
154      }
155      else
156      if ( (what==IMAP_CMD) || /*(*/ (what==FETCH_CMD) /*)*/) /* && (nMap!=nCopy)*/
157      {
158        for (i=R*C-1;i>=0;i--)
159        {
160          m->m[i]=p_PermPoly(((ideal)data)->m[i],perm,preimage_r,currRing,
161                          nMap,par_perm,P);
162          pTest(m->m[i]);
163        }
164      }
165      else /* if(what==MAP_CMD) */
166      {
167        matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r));
168        for (i=R*C-1;i>=0;i--)
169        {
170          m->m[i]=maEval(theMap, ((ideal)data)->m[i], preimage_r, nMap, (ideal)s, currRing);
171          pTest(m->m[i]);
172        }
173        idDelete((ideal *)&s);
174      }
175      if (nCoeff_is_Extension(currRing->cf))
176      {
177        for (i=R*C-1;i>=0;i--)
178        {
179          m->m[i]=p_MinPolyNormalize(m->m[i], currRing);
180          pTest(m->m[i]);
181        }
182      }
183      if(w->rtyp==MAP_CMD)
184      {
185        ((map)data)->preimage=tmpR;
186        ((map)m)->preimage=omStrDup(tmpR);
187      }
188      else
189      {
190        m->rank=((matrix)data)->rank;
191      }
192      res->data=(char *)m;
193      idTest((ideal) m);
194      break;
195    }
196
197    case LIST_CMD:
198    {
199      lists l=(lists)data;
200      lists ml=(lists)omAllocBin(slists_bin);
201      ml->Init(l->nr+1);
202      for(i=0;i<=l->nr;i++)
203      {
204        if (((l->m[i].rtyp>BEGIN_RING)&&(l->m[i].rtyp<END_RING))
205        ||(l->m[i].rtyp==LIST_CMD))
206        {
207          if (maApplyFetch(what,theMap,&ml->m[i],&l->m[i],
208                           preimage_r,perm,par_perm,P,nMap))
209          {
210            ml->Clean();
211            omFreeBin((ADDRESS)ml, slists_bin);
212            res->rtyp=0;
213            return TRUE;
214          }
215        }
216        else
217        {
218          ml->m[i].Copy(&l->m[i]);
219        }
220      }
221      res->data=(char *)ml;
222      break;
223    }
224    default:
225    {
226      return TRUE;
227    }
228  }
229  return FALSE;
230}
231
232/*2
233* substitutes the parameter par (from 1..N) by image,
234* does not destroy p and  image
235*/
236poly pSubstPar(poly p, int par, poly image)
237{
238  const ring R = currRing->cf->extRing;
239  ideal theMapI = idInit(rPar(currRing),1);
240  nMapFunc nMap = n_SetMap(R->cf, currRing->cf);
241  int i;
242  for(i = rPar(currRing);i>0;i--)
243  {
244    if (i != par)
245      theMapI->m[i-1]= p_NSet(n_Param(i, currRing), currRing);
246    else
247      theMapI->m[i-1] = p_Copy(image, currRing);
248    p_Test(theMapI->m[i-1],currRing);
249  }
250  //iiWriteMatrix((matrix)theMapI,"map:",1,currRing,0);
251
252  map theMap=(map)theMapI;
253  theMap->preimage=NULL;
254
255  leftv v=(leftv)omAllocBin(sleftv_bin);
256  sleftv tmpW;
257  poly res=NULL;
258
259  p_Normalize(p,currRing);
260  if (currRing->cf->rep==n_rep_rat_fct )
261  {
262    while (p!=NULL)
263    {
264      memset(v,0,sizeof(sleftv));
265
266      number d = n_GetDenom(p_GetCoeff(p, currRing), currRing);
267      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      memset(&tmpW,0,sizeof(sleftv));
281      tmpW.rtyp = POLY_CMD;
282      p_Test((poly)NUM(num), R);
283
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      //TODO check for memory leaks
294      poly pp = pHead(p);
295      //PrintS("map:");pWrite(pp);
296      if( d != NULL )
297      {
298        pSetCoeff(pp, n_Invers(d, currRing->cf));
299        n_Delete(&d, currRing); // d = NULL;
300      }
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  }
311  else if (currRing->cf->rep==n_rep_poly )
312  {
313    while (p!=NULL)
314    {
315      memset(v,0,sizeof(sleftv));
316
317      number num = n_GetNumerator(p_GetCoeff(p, currRing), currRing);
318      memset(&tmpW,0,sizeof(sleftv));
319      tmpW.rtyp = POLY_CMD;
320      p_Test((poly)num, R);
321
322
323      p_Normalize((poly)num,R);
324      if (num==NULL) num=(number)R->qideal->m[0];
325      tmpW.data = num; // a copy of this poly will be used
326      if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,R,NULL,NULL,0,nMap))
327      {
328        WerrorS("map failed");
329        v->data=NULL;
330      }
331      if (num!=(number)R->qideal->m[0]) n_Delete(&num, currRing);
332      //TODO check for memory leaks
333      poly pp = pHead(p);
334      //PrintS("map:");pWrite(pp);
335      pSetCoeff(pp,n_Init(1,currRing));
336      //PrintS("cf->");pWrite((poly)(v->data));
337      poly ppp = pMult((poly)(v->data),pp);
338      //PrintS("->");pWrite(ppp);
339      res=pAdd(res,ppp);
340      pIter(p);
341    }
342  }
343  else
344  {
345    WerrorS("cannot apply subst for these coeffcients");
346  }
347  idDelete((ideal *)(&theMap));
348  omFreeBin((ADDRESS)v, sleftv_bin);
349  return res;
350}
351
352/*2
353* substitute the n-th parameter by the poly e in id
354* does not destroy id and e
355*/
356ideal  idSubstPar(ideal id, int n, poly e)
357{
358  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
359  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
360
361  res->rank = id->rank;
362  for(k--;k>=0;k--)
363  {
364    res->m[k]=pSubstPar(id->m[k],n,e);
365  }
366  return res;
367}
368
369/*2
370* substitutes the variable var (from 1..N) by image,
371* does not destroy p and  image
372*/
373poly pSubstPoly(poly p, int var, poly image)
374{
375  if (p==NULL) return NULL;
376#ifdef HAVE_PLURAL
377  if (rIsPluralRing(currRing))
378  {
379    return pSubst(pCopy(p),var,image);
380  }
381#endif
382  map theMap=(map)idMaxIdeal(1);
383  theMap->preimage=NULL;
384  pDelete(&(theMap->m[var-1]));
385  theMap->m[var-1]=pCopy(image);
386
387  poly res=NULL;
388#ifdef FAST_MAP
389  if (pGetComp(p)==0)
390  {
391    ideal src_id=idInit(1,1);
392    src_id->m[0]=p;
393    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
394    res=res_id->m[0];
395    res_id->m[0]=NULL; idDelete(&res_id);
396    src_id->m[0]=NULL; idDelete(&src_id);
397  }
398  else
399#endif
400  {
401    sleftv tmpW;
402    memset(&tmpW,0,sizeof(sleftv));
403    tmpW.rtyp=POLY_CMD;
404    tmpW.data=p;
405    leftv v=(leftv)omAlloc0Bin(sleftv_bin);
406    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
407                            n_SetMap(currRing->cf, currRing->cf)))
408    {
409      WerrorS("map failed");
410      v->data=NULL;
411    }
412    res=(poly)(v->data);
413    omFreeBin((ADDRESS)v, sleftv_bin);
414  }
415  idDelete((ideal *)(&theMap));
416  return res;
417}
418
419/*2
420* substitute the n-th variable by the poly e in id
421* does not destroy id and e
422*/
423ideal  idSubstPoly(ideal id, int n, poly e)
424{
425
426#ifdef HAVE_PLURAL
427  if (rIsPluralRing(currRing))
428  {
429    int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
430    ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
431    res->rank = id->rank;
432    for(k--;k>=0;k--)
433    {
434      res->m[k]=pSubst(pCopy(id->m[k]),n,e);
435    }
436    return res;
437  }
438#endif
439  map theMap=(map)idMaxIdeal(1);
440  theMap->preimage=NULL;
441  pDelete(&(theMap->m[n-1]));
442  theMap->m[n-1]=pCopy(e);
443
444  leftv v=(leftv)omAlloc0Bin(sleftv_bin);
445  sleftv tmpW;
446  memset(&tmpW,0,sizeof(sleftv));
447  tmpW.rtyp=IDEAL_CMD;
448  tmpW.data=id;
449  if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
450                          n_SetMap(currRing->cf, currRing->cf)))
451  {
452    WerrorS("map failed");
453    v->data=NULL;
454  }
455  ideal res=(ideal)(v->data);
456  idDelete((ideal *)(&theMap));
457  omFreeBin((ADDRESS)v, sleftv_bin);
458  return res;
459}
Note: See TracBrowser for help on using the repository browser.