source: git/libpolys/polys/monomials/maps.cc @ 9d68fd

jengelh-datetimespielwiese
Last change on this file since 9d68fd was 9d68fd, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: better tests and normalization for subst
  • Property mode set to 100644
File size: 7.8 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    p_Test(p,preimage_r);
90    poly q=p_NSet(nMap(pGetCoeff(p),preimage_r->cf,dst_r->cf),dst_r);
91
92    int i;
93    for(i=1;i<=preimage_r->N; i++)
94    {
95      int pExp=p_GetExp( p,i,preimage_r);
96      if (pExp != 0)
97      {
98        if (theMap->m[i-1]!=NULL)
99        {
100          poly p1=theMap->m[i-1];
101          poly pp=maEvalVariable(p1,i,pExp,s,dst_r);
102          q = p_Mult_q(q,pp,dst_r);
103        }
104        else
105        {
106          p_Delete(&q,dst_r);
107          break;
108        }
109      }
110    }
111    int modulComp = p_GetComp( p,preimage_r);
112    if (q!=NULL) p_SetCompP(q,modulComp,dst_r);
113  return q;
114}
115
116poly maEval(map theMap, poly p,ring preimage_r,nMapFunc nMap, ideal s, const ring dst_r)
117{
118  poly result = NULL;
119  int i;
120
121//  for(i=1; i<=preimage_r->N; i++)
122//  {
123//    pTest(theMap->m[i-1]);
124//  }
125//  while (p!=NULL)
126//  {
127//    poly q=maEvalMonom(theMap,p,preimage_r,s);
128//    result = pAdd(result,q);
129//    pIter(p);
130//  }
131  if (p!=NULL)
132  {
133    int l = pLength(p)-1;
134    poly* monoms;
135    if (l>0)
136    {
137      monoms = (poly*) omAlloc(l*sizeof(poly));
138
139      for (i=0; i<l; i++)
140      {
141        monoms[i]=maEvalMonom(theMap,p,preimage_r,s, nMap, dst_r);
142        pIter(p);
143      }
144    }
145    result=maEvalMonom(theMap,p,preimage_r,s, nMap, dst_r);
146    if (l>0)
147    {
148      for(i = l-1; i>=0; i--)
149      {
150        result=p_Add_q(result, monoms[i], dst_r);
151      }
152      omFreeSize((ADDRESS)monoms,l*sizeof(poly));
153    }
154
155    assume(dst_r != NULL);
156    assume(dst_r->cf != NULL);
157   
158    if (nCoeff_is_algExt(dst_r->cf))
159      result = p_MinPolyNormalize(result, dst_r);
160  }
161  return result;
162}
163
164void maFindPerm(char const * const * const preim_names, int preim_n, char const * const * const preim_par, int preim_p,
165                char const * const * const names,       int n,       char const * const * const par,       int nop,
166                int * perm, int *par_perm, n_coeffType ch)
167{
168  int i,j;
169  /* find correspondig vars */
170  for (i=0; i<preim_n; i++)
171  {
172    for(j=0; j<n; j++)
173    {
174      if (strcmp(preim_names[i],names[j])==0)
175      {
176        if (BVERBOSE(V_IMAP))
177          Print("// var %s: nr %d -> nr %d\n",preim_names[i],i+1,j+1);
178        /* var i+1 from preimage ring is var j+1  (index j+1) from image ring */
179        perm[i+1]=j+1;
180        break;
181      }
182    }
183    if ((perm[i+1]==0)&&(par!=NULL)
184        // do not consider par of Fq
185         && (ch!=n_GF))
186    {
187      for(j=0; j<nop; j++)
188      {
189        if (strcmp(preim_names[i],par[j])==0)
190        {
191          if (BVERBOSE(V_IMAP))
192            Print("// var %s: nr %d -> par %d\n",preim_names[i],i+1,j+1);
193          /* var i+1 from preimage ring is par j+1 (index j) from image ring */
194          perm[i+1]=-(j+1);
195        }
196      }
197    }
198  }
199  if (par_perm!=NULL)
200  {
201    for (i=0; i<preim_p; i++)
202    {
203      for(j=0; j<n; j++)
204      {
205        if (strcmp(preim_par[i],names[j])==0)
206        {
207          if (BVERBOSE(V_IMAP))
208            Print("// par %s: par %d -> nr %d\n",preim_par[i],i+1,j+1);
209          /*par i+1 from preimage ring is var j+1  (index j+1) from image ring*/
210          par_perm[i]=j+1;
211          break;
212        }
213      }
214      if ((par!=NULL) && (par_perm[i]==0))
215      {
216        for(j=0; j<nop; j++)
217        {
218          if (strcmp(preim_par[i],par[j])==0)
219          {
220            if (BVERBOSE(V_IMAP))
221              Print("// par %s: nr %d -> par %d\n",preim_par[i],i+1,j+1);
222            /*par i+1 from preimage ring is par j+1 (index j) from image ring */
223            par_perm[i]=-(j+1);
224          }
225        }
226      }
227    }
228  }
229}
230
231/*2
232* embeds poly p from the subring r into the current ring
233*/
234poly maIMap(ring r, poly p, const ring dst_r)
235{
236  /* the simplest case:*/
237  if(r==dst_r) return p_Copy(p,dst_r);
238  nMapFunc nMap=n_SetMap(r->cf,dst_r->cf);
239  int *perm=(int *)omAlloc0((r->N+1)*sizeof(int));
240  //int *par_perm=(int *)omAlloc0(rPar(r)*sizeof(int));
241  maFindPerm(r->names, rVar(r), rParameter(r), rPar(r),
242             dst_r->names, rVar(dst_r),rParameter(dst_r), rPar(dst_r),
243             perm,NULL, dst_r->cf->type);
244  poly res=p_PermPoly(p,perm,r,dst_r, nMap /*,par_perm,rPar(r)*/);
245  omFreeSize((ADDRESS)perm,(r->N+1)*sizeof(int));
246  //omFreeSize((ADDRESS)par_perm,rPar(r)*sizeof(int));
247  return res;
248}
249
250/*3
251* find the max. degree in one variable, but not larger than MAX_MAP_DEG
252*/
253int maMaxDeg_Ma(ideal a,ring preimage_r)
254{
255  int i,j;
256  int N = preimage_r->N;
257  poly p;
258  int *m=(int *)omAlloc0(N*sizeof(int));
259
260  for (i=MATROWS(a)*MATCOLS(a)-1;i>=0;i--)
261  {
262    p=a->m[i];
263    //pTest(p); // cannot test p because it is from another ring
264    while(p!=NULL)
265    {
266      for(j=N-1;j>=0;j--)
267      {
268        m[j]=si_max(m[j],(int)p_GetExp( p,j+1,preimage_r));
269        if (m[j]>=MAX_MAP_DEG)
270        {
271          i=MAX_MAP_DEG;
272          goto max_deg_fertig_id;
273        }
274      }
275      pIter(p);
276    }
277  }
278  i=m[0];
279  for(j=N-1;j>0;j--)
280  {
281    i=si_max(i,m[j]);
282  }
283max_deg_fertig_id:
284  omFreeSize((ADDRESS)m,N*sizeof(int));
285  return i;
286}
287
288/*3
289* find the max. degree in one variable, but not larger than MAX_MAP_DEG
290*/
291int maMaxDeg_P(poly p,ring preimage_r)
292{
293  int i,j;
294  int N = preimage_r->N;
295  int *m=(int *)omAlloc0(N*sizeof(int));
296
297//  pTest(p);
298  while(p!=NULL)
299  {
300    for(j=N-1;j>=0;j--)
301    {
302      m[j]=si_max(m[j],(int)p_GetExp(p,j+1,preimage_r));
303      if (m[j]>=MAX_MAP_DEG)
304      {
305        i=MAX_MAP_DEG;
306        goto max_deg_fertig_p;
307      }
308    }
309    pIter(p);
310  }
311  i=m[0];
312  for(j=N-1;j>0;j--)
313  {
314    i=si_max(i,m[j]);
315  }
316max_deg_fertig_p:
317  omFreeSize((ADDRESS)m,N*sizeof(int));
318  return i;
319}
320
321// This is a very dirty way to cancel monoms whose number equals the
322// MinPoly
323poly p_MinPolyNormalize(poly p, const ring r)
324{
325  const coeffs C = r->cf;
326  number one = n_Init(1, C);
327  spolyrec rp;
328
329  poly q = &rp;
330
331  while (p != NULL)
332  {
333    // this returns 0, if p == MinPoly
334    number product = n_Mult(p_GetCoeff(p, r), one, C);
335    if ((product == NULL)||(n_IsZero(product, C)))
336    {
337      p_LmDelete(&p, r);
338    }
339    else
340    {
341      p_SetCoeff(p, product, r);
342      pNext(q) = p;
343      q = p;
344      p = pNext(p);
345    }
346  }
347  pNext(q) = NULL;
348  n_Delete(&one, C);
349  return rp.next;
350}
Note: See TracBrowser for help on using the repository browser.