source: git/libpolys/polys/monomials/maps.cc @ 6ce030f

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