source: git/kernel/maps.cc @ 98938c

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