source: git/libpolys/polys/monomials/maps.cc @ 01c1d0

spielwiese
Last change on this file since 01c1d0 was 01c1d0, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: moved P/minpoly/parameters from ring to coeff
  • Property mode set to 100644
File size: 8.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT - the mapping of polynomials to other rings
7*/
8
9#include <omalloc/omalloc.h>
10#include <misc/options.h>
11
12#include <coeffs/coeffs.h>
13#include <coeffs/numbers.h>
14
15#include <polys/monomials/p_polys.h>
16#include <polys/monomials/ring.h>
17#include <polys/simpleideals.h>
18#include <polys/prCopy.h>
19// #include <polys/ext_fields/longtrans.h>
20#include <polys/monomials/maps.h>
21
22#ifdef HAVE_PLURAL
23#include <polys/nc/nc.h>
24#endif
25
26// This is a very dirty way to "normalize" numbers w.r.t. a
27// MinPoly
28
29/* debug output: Tok2Cmdname in maApplyFetch*/
30//#include "ipshell.h"
31
32#define MAX_MAP_DEG 128
33
34/*2
35* copy a map
36*/
37map maCopy(map theMap, const ring r)
38{
39  int i;
40  map m=(map)idInit(IDELEMS(theMap),0);
41  for (i=IDELEMS(theMap)-1; i>=0; i--)
42      m->m[i] = p_Copy(theMap->m[i],r);
43  m->preimage=omStrDup(theMap->preimage);
44  return m;
45}
46
47
48/*2
49* return the image of var(v)^pExp, where var(v) maps to p
50*/
51poly maEvalVariable(poly p, int v,int pExp, ideal s, const ring dst_r)
52{
53  if (pExp==1)
54    return p_Copy(p,dst_r);
55
56  poly res;
57
58  if((s!=NULL)&&(pExp<MAX_MAP_DEG))
59  {
60    int j=2;
61    poly p0=p;
62    // find starting point
63    if(MATELEM(s,v,1)==NULL)
64    {
65      MATELEM(s,v,1)=p_Copy(p/*theMap->m[v-1]*/,dst_r);
66    }
67    else
68    {
69      while((j<=pExp)&&(MATELEM(s,v,j)!=NULL))
70      {
71        j++;
72      }
73      p0=MATELEM(s,v,j-1);
74    }
75    // multiply
76    for(;j<=pExp;j++)
77    {
78      p0=MATELEM(s,v,j)=pp_Mult_qq(p0, p,dst_r);
79      p_Normalize(p0, dst_r);
80    }
81    res=p_Copy(p0/*MATELEM(s,v,pExp)*/,dst_r);
82  }
83  else //if ((p->next!=NULL)&&(p->next->next==NULL))
84  {
85    res=p_Power(p_Copy(p,dst_r),pExp,dst_r);
86  }
87  return res;
88}
89
90static poly maEvalMonom(map theMap, poly p,ring preimage_r, ideal s, 
91           nMapFunc nMap, const ring dst_r)
92{
93    poly q=p_NSet(nMap(pGetCoeff(p),preimage_r->cf,dst_r->cf),dst_r);
94
95    int i;
96    for(i=1;i<=preimage_r->N; i++)
97    {
98      int pExp=p_GetExp( p,i,preimage_r);
99      if (pExp != 0)
100      {
101        if (theMap->m[i-1]!=NULL)
102        {
103          poly p1=theMap->m[i-1];
104          poly pp=maEvalVariable(p1,i,pExp,s,dst_r);
105          q = p_Mult_q(q,pp,dst_r);
106        }
107        else
108        {
109          p_Delete(&q,dst_r);
110          break;
111        }
112      }
113    }
114    int modulComp = p_GetComp( p,preimage_r);
115    if (q!=NULL) p_SetCompP(q,modulComp,dst_r);
116  return q;
117}
118
119poly maEval(map theMap, poly p,ring preimage_r,nMapFunc nMap, ideal s, const ring dst_r)
120{
121  poly result = NULL;
122  int i;
123
124//  for(i=1; i<=preimage_r->N; i++)
125//  {
126//    pTest(theMap->m[i-1]);
127//  }
128//  while (p!=NULL)
129//  {
130//    poly q=maEvalMonom(theMap,p,preimage_r,s);
131//    result = pAdd(result,q);
132//    pIter(p);
133//  }
134  if (p!=NULL)
135  {
136    int l = pLength(p)-1;
137    poly* monoms;
138    if (l>0)
139    {
140      monoms = (poly*) omAlloc(l*sizeof(poly));
141
142      for (i=0; i<l; i++)
143      {
144        monoms[i]=maEvalMonom(theMap,p,preimage_r,s, nMap, dst_r);
145        pIter(p);
146      }
147    }
148    result=maEvalMonom(theMap,p,preimage_r,s, nMap, dst_r);
149    if (l>0)
150    {
151      for(i = l-1; i>=0; i--)
152      {
153        result=p_Add_q(result, monoms[i], dst_r);
154      }
155      omFreeSize((ADDRESS)monoms,l*sizeof(poly));
156    }
157    if (dst_r->cf->minpoly!=NULL) result=p_MinPolyNormalize(result, dst_r);
158  }
159  return result;
160}
161
162/*2
163*shifts the variables between minvar and maxvar of p  \in p_ring to the
164*first maxvar-minvar+1 variables in the actual ring
165*be carefull: there is no range check for the variables of p
166*/
167static poly pChangeSizeOfPoly(ring p_ring, poly p,int minvar,int maxvar, const ring dst_r)
168{
169  int i;
170  poly result = NULL,resultWorkP;
171  number n;
172
173  if (p==NULL) return result;
174  else result = p_Init(dst_r);
175  resultWorkP = result;
176  while (p!=NULL)
177  {
178    for (i=minvar;i<=maxvar;i++)
179      p_SetExp(resultWorkP,i-minvar+1,p_GetExp(p,i,p_ring),dst_r);
180    p_SetComp(resultWorkP,p_GetComp(p,p_ring),dst_r);
181    n=n_Copy(pGetCoeff(p),dst_r->cf);
182    p_SetCoeff(resultWorkP,n,dst_r);
183    p_Setm(resultWorkP,dst_r);
184    pIter(p);
185    if (p!=NULL)
186    {
187      pNext(resultWorkP) = p_Init(dst_r);
188      pIter(resultWorkP);
189    }
190  }
191  return result;
192}
193
194
195void maFindPerm(char **preim_names, int preim_n, char **preim_par, int preim_p,
196                char **names,       int n,       char **par,       int nop,
197                int * perm, int *par_perm, int ch)
198{
199  int i,j;
200  /* find correspondig vars */
201  for (i=0; i<preim_n; i++)
202  {
203    for(j=0; j<n; j++)
204    {
205      if (strcmp(preim_names[i],names[j])==0)
206      {
207        if (BVERBOSE(V_IMAP))
208          Print("// var %s: nr %d -> nr %d\n",preim_names[i],i+1,j+1);
209        /* var i+1 from preimage ring is var j+1  (index j+1) from image ring */
210        perm[i+1]=j+1;
211        break;
212      }
213    }
214    if ((perm[i+1]==0)&&(par!=NULL)
215        // do not consider par of Fq
216         && (ch < 2))
217    {
218      for(j=0; j<nop; j++)
219      {
220        if (strcmp(preim_names[i],par[j])==0)
221        {
222          if (BVERBOSE(V_IMAP))
223            Print("// var %s: nr %d -> par %d\n",preim_names[i],i+1,j+1);
224          /* var i+1 from preimage ring is par j+1 (index j) from image ring */
225          perm[i+1]=-(j+1);
226        }
227      }
228    }
229  }
230  if (par_perm!=NULL)
231  {
232    for (i=0; i<preim_p; i++)
233    {
234      for(j=0; j<n; j++)
235      {
236        if (strcmp(preim_par[i],names[j])==0)
237        {
238          if (BVERBOSE(V_IMAP))
239            Print("// par %s: par %d -> nr %d\n",preim_par[i],i+1,j+1);
240          /*par i+1 from preimage ring is var j+1  (index j+1) from image ring*/
241          par_perm[i]=j+1;
242          break;
243        }
244      }
245      if ((par!=NULL) && (par_perm[i]==0))
246      {
247        for(j=0; j<nop; j++)
248        {
249          if (strcmp(preim_par[i],par[j])==0)
250          {
251            if (BVERBOSE(V_IMAP))
252              Print("// par %s: nr %d -> par %d\n",preim_par[i],i+1,j+1);
253            /*par i+1 from preimage ring is par j+1 (index j) from image ring */
254            par_perm[i]=-(j+1);
255          }
256        }
257      }
258    }
259  }
260}
261
262/*2
263* embeds poly p from the subring r into the current ring
264*/
265poly maIMap(ring r, poly p, const ring dst_r)
266{
267  /* the simplest case:*/
268  if(r==dst_r) return p_Copy(p,dst_r);
269  nMapFunc nMap=n_SetMap(r->cf,dst_r->cf);
270  int *perm=(int *)omAlloc0((r->N+1)*sizeof(int));
271  //int *par_perm=(int *)omAlloc0(rPar(r)*sizeof(int));
272  maFindPerm(r->names, rVar(r), r->cf->parameter, rPar(r),
273             dst_r->names, rVar(dst_r),dst_r->cf->parameter, rPar(dst_r),
274             perm,NULL, dst_r->cf->ch);
275  poly res=p_PermPoly(p,perm,r,dst_r, nMap /*,par_perm,rPar(r)*/);
276  omFreeSize((ADDRESS)perm,(r->N+1)*sizeof(int));
277  //omFreeSize((ADDRESS)par_perm,rPar(r)*sizeof(int));
278  return res;
279}
280
281/*3
282* find the max. degree in one variable, but not larger than MAX_MAP_DEG
283*/
284int maMaxDeg_Ma(ideal a,ring preimage_r)
285{
286  int i,j;
287  int N = preimage_r->N;
288  poly p;
289  int *m=(int *)omAlloc0(N*sizeof(int));
290
291  for (i=MATROWS(a)*MATCOLS(a)-1;i>=0;i--)
292  {
293    p=a->m[i];
294    //pTest(p); // cannot test p because it is from another ring
295    while(p!=NULL)
296    {
297      for(j=N-1;j>=0;j--)
298      {
299        m[j]=si_max(m[j],(int)p_GetExp( p,j+1,preimage_r));
300        if (m[j]>=MAX_MAP_DEG)
301        {
302          i=MAX_MAP_DEG;
303          goto max_deg_fertig_id;
304        }
305      }
306      pIter(p);
307    }
308  }
309  i=m[0];
310  for(j=N-1;j>0;j--)
311  {
312    i=si_max(i,m[j]);
313  }
314max_deg_fertig_id:
315  omFreeSize((ADDRESS)m,N*sizeof(int));
316  return i;
317}
318
319/*3
320* find the max. degree in one variable, but not larger than MAX_MAP_DEG
321*/
322int maMaxDeg_P(poly p,ring preimage_r)
323{
324  int i,j;
325  int N = preimage_r->N;
326  int *m=(int *)omAlloc0(N*sizeof(int));
327
328//  pTest(p);
329  while(p!=NULL)
330  {
331    for(j=N-1;j>=0;j--)
332    {
333      m[j]=si_max(m[j],(int)p_GetExp(p,j+1,preimage_r));
334      if (m[j]>=MAX_MAP_DEG)
335      {
336        i=MAX_MAP_DEG;
337        goto max_deg_fertig_p;
338      }
339    }
340    pIter(p);
341  }
342  i=m[0];
343  for(j=N-1;j>0;j--)
344  {
345    i=si_max(i,m[j]);
346  }
347max_deg_fertig_p:
348  omFreeSize((ADDRESS)m,N*sizeof(int));
349  return i;
350}
351
352// This is a very dirty way to cancel monoms whose number equals the
353// MinPoly
354poly p_MinPolyNormalize(poly p, const ring r)
355{
356  number one = n_Init(1, r->cf);
357  spolyrec rp;
358
359  poly q = &rp;
360
361  while (p != NULL)
362  {
363    // this returns 0, if p == MinPoly
364    number product = n_Mult(pGetCoeff(p), one,r->cf);
365    if ((product == NULL)||(n_IsZero(product,r->cf)))
366    {
367      p_LmDelete(&p,r);
368    }
369    else
370    {
371      p_SetCoeff(p, product,r);
372      pNext(q) = p;
373      q = p;
374      p = pNext(p);
375    }
376  }
377  pNext(q) = NULL;
378  return rp.next;
379}
Note: See TracBrowser for help on using the repository browser.