source: git/Singular/maps_ip.cc @ f769fd0

jengelh-datetimespielwiese
Last change on this file since f769fd0 was f769fd0, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
Fixed the handling of constant (non-1) denoms
  • Property mode set to 100644
File size: 10.2 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        matrix s=mpNew(N,maMaxDeg_P((poly)data, preimage_r));
122        res->data=(void *)maEval(theMap, (poly)data, preimage_r, nMap, (ideal)s, currRing);
123        idDelete((ideal *)&s);
124      }
125      if (nCoeff_is_Extension(currRing->cf))
126        res->data=(void *)p_MinPolyNormalize((poly)res->data, currRing);
127      pTest((poly)res->data);
128      break;
129    case MODUL_CMD:
130    case MATRIX_CMD:
131    case IDEAL_CMD:
132    case MAP_CMD:
133    {
134      int C=((matrix)data)->cols();
135      int R;
136      if (w->rtyp==MAP_CMD) R=1;
137      else R=((matrix)data)->rows();
138      matrix m=mpNew(R,C);
139      char *tmpR=NULL;
140      if(w->rtyp==MAP_CMD)
141      {
142        tmpR=((map)data)->preimage;
143        ((matrix)data)->rank=((matrix)data)->rows();
144      }
145      if ((what==FETCH_CMD)&& (preimage_r->cf == currRing->cf))
146      {
147        for (i=R*C-1;i>=0;i--)
148        {
149          m->m[i]=prCopyR(((ideal)data)->m[i], preimage_r, currRing);
150          pTest(m->m[i]);
151        }
152      }
153      else
154      if ((what==IMAP_CMD) || ((what==FETCH_CMD) /* && (nMap!=nCopy)*/))
155      {
156        for (i=R*C-1;i>=0;i--)
157        {
158          m->m[i]=p_PermPoly(((ideal)data)->m[i],perm,preimage_r,currRing,
159                          nMap,par_perm,P);
160          pTest(m->m[i]);
161        }
162      }
163      else /* if(what==MAP_CMD) */
164      {
165        matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r));
166        for (i=R*C-1;i>=0;i--)
167        {
168          m->m[i]=maEval(theMap, ((ideal)data)->m[i], preimage_r, nMap, (ideal)s, currRing);
169          pTest(m->m[i]);
170        }
171        idDelete((ideal *)&s);
172      }
173      if (nCoeff_is_Extension(currRing->cf))
174      {
175        for (i=R*C-1;i>=0;i--)
176        {
177          m->m[i]=p_MinPolyNormalize(m->m[i], currRing);
178          pTest(m->m[i]);
179        }
180      }
181      if(w->rtyp==MAP_CMD)
182      {
183        ((map)data)->preimage=tmpR;
184        ((map)m)->preimage=omStrDup(tmpR);
185      }
186      else
187      {
188        m->rank=((matrix)data)->rank;
189      }
190      res->data=(char *)m;
191      idTest((ideal) m);
192      break;
193    }
194
195    case LIST_CMD:
196    {
197      lists l=(lists)data;
198      lists ml=(lists)omAllocBin(slists_bin);
199      ml->Init(l->nr+1);
200      for(i=0;i<=l->nr;i++)
201      {
202        if (((l->m[i].rtyp>BEGIN_RING)&&(l->m[i].rtyp<END_RING))
203        ||(l->m[i].rtyp==LIST_CMD))
204        {
205          if (maApplyFetch(what,theMap,&ml->m[i],&l->m[i],
206                           preimage_r,perm,par_perm,P,nMap))
207          {
208            ml->Clean();
209            omFreeBin((ADDRESS)ml, slists_bin);
210            res->rtyp=0;
211            return TRUE;
212          }
213        }
214        else
215        {
216          ml->m[i].Copy(&l->m[i]);
217        }
218      }
219      res->data=(char *)ml;
220      break;
221    }
222    default:
223    {
224      return TRUE;
225    }
226  }
227  return FALSE;
228}
229
230/*2
231* substitutes the parameter par (from 1..N) by image,
232* does not destroy p and  image
233*/
234poly pSubstPar(poly p, int par, poly image)
235{
236  assume( nCoeff_is_transExt(currRing->cf) ); // nCoeff_is_Extension???
237  const ring R = currRing->cf->extRing;
238
239  ideal theMapI = idInit(rPar(currRing),1);
240  nMapFunc nMap = n_SetMap(R->cf, currRing->cf);
241
242  int i;
243  for(i = rPar(currRing);i>0;i--)
244  {
245    if (i != par)
246      theMapI->m[i-1]= p_NSet(n_Param(i, currRing), currRing);
247    else
248      theMapI->m[i-1] = p_Copy(image, currRing);
249  }
250
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  while (p!=NULL)
260  {
261    memset(v,0,sizeof(sleftv));
262
263    number d = n_GetDenom(p_GetCoeff(p, currRing), currRing);
264    assume( p_Test((poly)NUM(d), R) );
265
266    if ( n_IsOne (d, currRing->cf) )
267    {     
268      n_Delete(&d, currRing); d = NULL;
269    }     
270    else if (!p_IsConstant((poly)NUM(d), R))
271    {   
272      WarnS("ignoring denominators of coefficients...");
273      n_Delete(&d, currRing); d = NULL;
274    }     
275
276    number num = n_GetNumerator(p_GetCoeff(p, currRing), currRing);
277    assume( p_Test((poly)NUM(num), R) );
278
279    memset(&tmpW,0,sizeof(sleftv));
280    tmpW.rtyp = POLY_CMD;
281    tmpW.data = NUM (num); // a copy of this poly will be used
282
283    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,R,NULL,NULL,0,nMap))
284    {
285      WerrorS("map failed");
286      v->data=NULL;
287    }
288    n_Delete(&num, currRing);
289
290    //TODO check for memory leaks
291    poly pp = pHead(p);
292    //PrintS("map:");pWrite(pp);
293    if( d != NULL )
294    {
295      pSetCoeff(pp, n_Invers(d, currRing->cf));
296      n_Delete(&d, currRing); // d = NULL;
297    } else
298      pSetCoeff(pp, nInit(1));
299
300    //PrintS("->");pWrite((poly)(v->data));
301    poly ppp = pMult((poly)(v->data),pp);
302    //PrintS("->");pWrite(ppp);
303    res=pAdd(res,ppp);
304    pIter(p);
305  }
306  idDelete((ideal *)(&theMap));
307  omFreeBin((ADDRESS)v, sleftv_bin);
308  return res;
309}
310
311/*2
312* substitute the n-th parameter by the poly e in id
313* does not destroy id and e
314*/
315ideal  idSubstPar(ideal id, int n, poly e)
316{
317  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
318  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
319
320  res->rank = id->rank;
321  for(k--;k>=0;k--)
322  {
323    res->m[k]=pSubstPar(id->m[k],n,e);
324  }
325  return res;
326}
327
328/*2
329* substitutes the variable var (from 1..N) by image,
330* does not destroy p and  image
331*/
332poly pSubstPoly(poly p, int var, poly image)
333{
334  if (p==NULL) return NULL;
335#ifdef HAVE_PLURAL
336  if (rIsPluralRing(currRing))
337  {
338    return pSubst(pCopy(p),var,image);
339  }
340#endif
341  map theMap=(map)idMaxIdeal(1);
342  theMap->preimage=NULL;
343  pDelete(&(theMap->m[var-1]));
344  theMap->m[var-1]=pCopy(image);
345
346  poly res=NULL;
347#ifdef FAST_MAP
348  if (pGetComp(p)==0)
349  {
350    ideal src_id=idInit(1,1);
351    src_id->m[0]=p;
352    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
353    res=res_id->m[0];
354    res_id->m[0]=NULL; idDelete(&res_id);
355    src_id->m[0]=NULL; idDelete(&src_id);
356  }
357  else
358#endif
359  {
360    sleftv tmpW;
361    memset(&tmpW,0,sizeof(sleftv));
362    tmpW.rtyp=POLY_CMD;
363    tmpW.data=p;
364    leftv v=(leftv)omAlloc0Bin(sleftv_bin);
365    if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
366                            n_SetMap(currRing->cf, currRing->cf)))
367    {
368      WerrorS("map failed");
369      v->data=NULL;
370    }
371    res=(poly)(v->data);
372    omFreeBin((ADDRESS)v, sleftv_bin);
373  }
374  idDelete((ideal *)(&theMap));
375  return res;
376}
377
378/*2
379* substitute the n-th variable by the poly e in id
380* does not destroy id and e
381*/
382ideal  idSubstPoly(ideal id, int n, poly e)
383{
384
385#ifdef HAVE_PLURAL
386  if (rIsPluralRing(currRing))
387  {
388    int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
389    ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
390    res->rank = id->rank;
391    for(k--;k>=0;k--)
392    {
393      res->m[k]=pSubst(pCopy(id->m[k]),n,e);
394    }
395    return res;
396  }
397#endif
398  map theMap=(map)idMaxIdeal(1);
399  theMap->preimage=NULL;
400  pDelete(&(theMap->m[n-1]));
401  theMap->m[n-1]=pCopy(e);
402
403  leftv v=(leftv)omAlloc0Bin(sleftv_bin);
404  sleftv tmpW;
405  memset(&tmpW,0,sizeof(sleftv));
406  tmpW.rtyp=IDEAL_CMD;
407  tmpW.data=id;
408  if (maApplyFetch(MAP_CMD,theMap,v,&tmpW,currRing,NULL,NULL,0,
409                          n_SetMap(currRing->cf, currRing->cf)))
410  {
411    WerrorS("map failed");
412    v->data=NULL;
413  }
414  ideal res=(ideal)(v->data);
415  idDelete((ideal *)(&theMap));
416  omFreeBin((ADDRESS)v, sleftv_bin);
417  return res;
418}
Note: See TracBrowser for help on using the repository browser.