source: git/kernel/maps.cc @ 599326

spielwiese
Last change on this file since 599326 was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.6 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 <kernel/mod2.h>
10#include <kernel/options.h>
11#include <kernel/febase.h>
12#include <kernel/polys.h>
13#include <kernel/numbers.h>
14#include <kernel/ring.h>
15#include <kernel/ideals.h>
16#include <kernel/matpol.h>
17#include <omalloc.h>
18#include <kernel/kstd1.h>
19#include <kernel/longalg.h>
20#include <kernel/maps.h>
21#include <kernel/prCopy.h>
22
23#ifdef HAVE_PLURAL
24#include <kernel/gring.h>
25#endif
26
27// This is a very dirty way to "normalize" numbers w.r.t. a
28// MinPoly
29
30/* debug output: Tok2Cmdname in maApplyFetch*/
31//#include "ipshell.h"
32
33#define MAX_MAP_DEG 128
34
35/*2
36* copy a map
37*/
38map maCopy(map theMap)
39{
40  int i;
41  map m=(map)idInit(IDELEMS(theMap),0);
42  for (i=IDELEMS(theMap)-1; i>=0; i--)
43      m->m[i] = pCopy(theMap->m[i]);
44  m->preimage=omStrDup(theMap->preimage);
45  return m;
46}
47
48
49/*2
50* return the image of var(v)^pExp, where var(v) maps to p
51*/
52poly maEvalVariable(poly p, int v,int pExp,matrix s)
53{
54  if (pExp==1)
55    return pCopy(p);
56
57  poly res;
58
59  if((s!=NULL)&&(pExp<MAX_MAP_DEG))
60  {
61    int j=2;
62    poly p0=p;
63    // find starting point
64    if(MATELEM(s,v,1)==NULL)
65    {
66      MATELEM(s,v,1)=pCopy(p/*theMap->m[v-1]*/);
67    }
68    else
69    {
70      while((j<=pExp)&&(MATELEM(s,v,j)!=NULL))
71      {
72        j++;
73      }
74      p0=MATELEM(s,v,j-1);
75    }
76    // multiply
77    for(;j<=pExp;j++)
78    {
79      p0=MATELEM(s,v,j)=ppMult_qq(p0, p);
80      pNormalize(p0);
81    }
82    res=pCopy(p0/*MATELEM(s,v,pExp)*/);
83  }
84  else //if ((p->next!=NULL)&&(p->next->next==NULL))
85  {
86    res=pPower(pCopy(p),pExp);
87  }
88  return res;
89}
90
91static poly maEvalMonom(map theMap, poly p,ring preimage_r,matrix s, nMapFunc nMap)
92{
93    poly q=pNSet(nMap(pGetCoeff(p)));
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);
105          q = pMult(q,pp);
106        }
107        else
108        {
109          pDelete(&q);
110          break;
111        }
112      }
113    }
114    int modulComp = p_GetComp( p,preimage_r);
115    if (q!=NULL) pSetCompP(q,modulComp);
116  return q;
117}
118
119poly maEval(map theMap, poly p,ring preimage_r,nMapFunc nMap,matrix s)
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);
145        pIter(p);
146      }
147    }
148    result=maEvalMonom(theMap,p,preimage_r,s, nMap);
149    if (l>0)
150    {
151      for(i = l-1; i>=0; i--)
152      {
153        result=pAdd(result, monoms[i]);
154      }
155      omFreeSize((ADDRESS)monoms,l*sizeof(poly));
156    }
157    if (currRing->minpoly!=NULL) result=pMinPolyNormalize(result);
158  }
159  pTest(result);
160  return result;
161}
162
163/*2
164*shifts the variables between minvar and maxvar of p  \in p_ring to the
165*first maxvar-minvar+1 variables in the actual ring
166*be carefull: there is no range check for the variables of p
167*/
168static poly pChangeSizeOfPoly(ring p_ring, poly p,int minvar,int maxvar)
169{
170  int i;
171  poly result = NULL,resultWorkP;
172  number n;
173
174  if (p==NULL) return result;
175  else result = pInit();
176  resultWorkP = result;
177  while (p!=NULL)
178  {
179    for (i=minvar;i<=maxvar;i++)
180      pSetExp(resultWorkP,i-minvar+1,p_GetExp(p,i,p_ring));
181    pSetComp(resultWorkP,p_GetComp(p,p_ring));
182    n=nCopy(pGetCoeff(p));
183    pSetCoeff(resultWorkP,n);
184    pSetm(resultWorkP);
185    pIter(p);
186    if (p!=NULL)
187    {
188      pNext(resultWorkP) = pInit();
189      pIter(resultWorkP);
190    }
191  }
192  return result;
193}
194
195
196/*2
197*returns the preimage of id under theMap,
198*if id is empty or zero the kernel is computed
199* (assumes) that both ring have the same coeff.field
200*/
201ideal maGetPreimage(ring theImageRing, map theMap, ideal id)
202{
203  ring sourcering = currRing;
204
205#ifdef HAVE_PLURAL
206  if (rIsPluralRing(theImageRing))
207  {
208    if ((rIsPluralRing(sourcering)) && (ncRingType(sourcering)!=nc_comm)) 
209    {
210      Werror("Sorry, not yet implemented for noncomm. rings");
211      return NULL;
212    }
213  }
214#endif
215 
216  int i,j;
217  poly p,pp,q;
218  ideal temp1;
219  ideal temp2;
220
221  int imagepvariables = theImageRing->N;
222  int N = pVariables+imagepvariables;
223
224  ring tmpR;
225  if (rSumInternal(theImageRing,sourcering,tmpR,FALSE,TRUE)!=1)
226  {
227     WerrorS("error in rSumInternal");
228     return NULL;
229  }
230
231  if (nSetMap(theImageRing) != nCopy)
232  {
233    Werror("Coefficient fields must be equal");
234    return NULL;
235  }
236
237  // change to new ring
238  rChangeCurrRing(tmpR);
239  if (id==NULL)
240    j = 0;
241  else
242    j = IDELEMS(id);
243  int j0=j;
244  if (theImageRing->qideal!=NULL) j+=IDELEMS(theImageRing->qideal);
245  temp1 = idInit(sourcering->N+j,1);
246  for (i=0;i<sourcering->N;i++)
247  {
248    q = pISet(-1);
249    pSetExp(q,i+1+imagepvariables,1);
250    pSetm(q);
251    if ((i<IDELEMS(theMap)) && (theMap->m[i]!=NULL))
252    {
253      p = pSort(pChangeSizeOfPoly(theImageRing,theMap->m[i],1,imagepvariables));
254      p=pAdd(p,q);
255    }
256    else
257    {
258      p = q;
259    }
260    temp1->m[i] = p;
261  }
262  idTest(temp1);
263  for (i=sourcering->N;i<sourcering->N+j0;i++)
264  {
265    temp1->m[i] = pSort(pChangeSizeOfPoly(theImageRing,
266                                    id->m[i-sourcering->N],1,imagepvariables));
267  }
268  for (i=sourcering->N+j0;i<sourcering->N+j;i++)
269  {
270    temp1->m[i] = pSort(pChangeSizeOfPoly(theImageRing,
271                                    theImageRing->qideal->m[i-sourcering->N-j0],
272                                    1,imagepvariables));
273  }
274  // we ignore here homogenity - may be changed later:
275  temp2 = kStd(temp1,NULL,isNotHomog,NULL);
276  idDelete(&temp1);
277  for (i=0;i<IDELEMS(temp2);i++)
278  {
279    if (pLowVar(temp2->m[i])<imagepvariables) pDelete(&(temp2->m[i]));
280  }
281
282  // let's get back to the original ring
283  rChangeCurrRing(sourcering);
284  temp1 = idInit(5,1);
285  j = 0;
286  for (i=0;i<IDELEMS(temp2);i++)
287  {
288    p = temp2->m[i];
289    if (p!=NULL)
290    {
291      q = pSort(pChangeSizeOfPoly(tmpR, p,imagepvariables+1,N));
292      if (j>=IDELEMS(temp1))
293      {
294        pEnlargeSet(&(temp1->m),IDELEMS(temp1),5);
295        IDELEMS(temp1)+=5;
296      }
297      temp1->m[j] = q;
298      j++;
299    }
300  }
301  id_Delete(&temp2, tmpR);
302  idSkipZeroes(temp1);
303  rKill(tmpR);
304  return temp1;
305}
306
307void maFindPerm(char **preim_names, int preim_n, char **preim_par, int preim_p,
308                char **names,       int n,       char **par,       int nop,
309                int * perm, int *par_perm, int ch)
310{
311  int i,j;
312  /* find correspondig vars */
313  for (i=0; i<preim_n; i++)
314  {
315    for(j=0; j<n; j++)
316    {
317      if (strcmp(preim_names[i],names[j])==0)
318      {
319        if (BVERBOSE(V_IMAP))
320          Print("// var %s: nr %d -> nr %d\n",preim_names[i],i+1,j+1);
321        /* var i+1 from preimage ring is var j+1  (index j+1) from image ring */
322        perm[i+1]=j+1;
323        break;
324      }
325    }
326    if ((perm[i+1]==0)&&(par!=NULL)
327        // do not consider par of Fq
328         && (ch < 2))
329    {
330      for(j=0; j<nop; j++)
331      {
332        if (strcmp(preim_names[i],par[j])==0)
333        {
334          if (BVERBOSE(V_IMAP))
335            Print("// var %s: nr %d -> par %d\n",preim_names[i],i+1,j+1);
336          /* var i+1 from preimage ring is par j+1 (index j) from image ring */
337          perm[i+1]=-(j+1);
338        }
339      }
340    }
341  }
342  if (par_perm!=NULL)
343  {
344    for (i=0; i<preim_p; i++)
345    {
346      for(j=0; j<n; j++)
347      {
348        if (strcmp(preim_par[i],names[j])==0)
349        {
350          if (BVERBOSE(V_IMAP))
351            Print("// par %s: par %d -> nr %d\n",preim_par[i],i+1,j+1);
352          /*par i+1 from preimage ring is var j+1  (index j+1) from image ring*/
353          par_perm[i]=j+1;
354          break;
355        }
356      }
357      if ((par!=NULL) && (par_perm[i]==0))
358      {
359        for(j=0; j<nop; j++)
360        {
361          if (strcmp(preim_par[i],par[j])==0)
362          {
363            if (BVERBOSE(V_IMAP))
364              Print("// par %s: nr %d -> par %d\n",preim_par[i],i+1,j+1);
365            /*par i+1 from preimage ring is par j+1 (index j) from image ring */
366            par_perm[i]=-(j+1);
367          }
368        }
369      }
370    }
371  }
372}
373
374/*2
375* embeds poly p from the subring r into the current ring
376*/
377poly maIMap(ring r, poly p)
378{
379  /* the simplest case:*/
380  if(r==currRing) return pCopy(p);
381  nMapFunc nMap=nSetMap(r);
382  int *perm=(int *)omAlloc0((r->N+1)*sizeof(int));
383  //int *par_perm=(int *)omAlloc0(rPar(r)*sizeof(int));
384  maFindPerm(r->names,r->N, r->parameter, r->P,
385             currRing->names,currRing->N,currRing->parameter, currRing->P,
386             perm,NULL, currRing->ch);
387  poly res=pPermPoly(p,perm,r, nMap /*,par_perm,rPar(r)*/);
388  omFreeSize((ADDRESS)perm,(r->N+1)*sizeof(int));
389  //omFreeSize((ADDRESS)par_perm,rPar(r)*sizeof(int));
390  return res;
391}
392
393/*3
394* find the max. degree in one variable, but not larger than MAX_MAP_DEG
395*/
396int maMaxDeg_Ma(ideal a,ring preimage_r)
397{
398  int i,j;
399  int N = preimage_r->N;
400  poly p;
401  int *m=(int *)omAlloc0(N*sizeof(int));
402
403  for (i=MATROWS(a)*MATCOLS(a)-1;i>=0;i--)
404  {
405    p=a->m[i];
406    //pTest(p); // cannot test p because it is from another ring
407    while(p!=NULL)
408    {
409      for(j=N-1;j>=0;j--)
410      {
411        m[j]=si_max(m[j],(int)p_GetExp( p,j+1,preimage_r));
412        if (m[j]>=MAX_MAP_DEG)
413        {
414          i=MAX_MAP_DEG;
415          goto max_deg_fertig_id;
416        }
417      }
418      pIter(p);
419    }
420  }
421  i=m[0];
422  for(j=N-1;j>0;j--)
423  {
424    i=si_max(i,m[j]);
425  }
426max_deg_fertig_id:
427  omFreeSize((ADDRESS)m,N*sizeof(int));
428  return i;
429}
430
431/*3
432* find the max. degree in one variable, but not larger than MAX_MAP_DEG
433*/
434int maMaxDeg_P(poly p,ring preimage_r)
435{
436  int i,j;
437  int N = preimage_r->N;
438  int *m=(int *)omAlloc0(N*sizeof(int));
439
440//  pTest(p);
441  while(p!=NULL)
442  {
443    for(j=N-1;j>=0;j--)
444    {
445      m[j]=si_max(m[j],(int)p_GetExp(p,j+1,preimage_r));
446      if (m[j]>=MAX_MAP_DEG)
447      {
448        i=MAX_MAP_DEG;
449        goto max_deg_fertig_p;
450      }
451    }
452    pIter(p);
453  }
454  i=m[0];
455  for(j=N-1;j>0;j--)
456  {
457    i=si_max(i,m[j]);
458  }
459max_deg_fertig_p:
460  omFreeSize((ADDRESS)m,N*sizeof(int));
461  return i;
462}
463
464// This is a very dirty way to cancel monoms whose number equals the
465// MinPoly
466poly pMinPolyNormalize(poly p)
467{
468  number one = nInit(1);
469  spolyrec rp;
470
471  poly q = &rp;
472
473  while (p != NULL)
474  {
475    // this returns 0, if p == MinPoly
476    number product = nMult(pGetCoeff(p), one);
477    if (product == NULL)
478    {
479      pLmDelete(&p);
480    }
481    else
482    {
483      pSetCoeff(p, product);
484      pNext(q) = p;
485      q = p;
486      p = pNext(p);
487    }
488  }
489  pNext(q) = NULL;
490  return rp.next;
491}
Note: See TracBrowser for help on using the repository browser.