source: git/Singular/LIB/resgraph.lib @ c47090

spielwiese
Last change on this file since c47090 was c47090, checked in by jgmbenoit <quatermaster@…>, 9 years ago
correct spelling error in Singular lib resgraph.lib
  • Property mode set to 100644
File size: 25.9 KB
Line 
1////////////////////////////////////////////////////////////////////////////
2version="version resgraph.lib 4.0.0.0 Jun_2013 "; // $Id$
3category="Visualization";
4info="
5LIBRARY:  resgraph.lib            Visualization of Resolution Data
6
7AUTHOR: A. Fruehbis-Krueger,  anne@mathematik.uni-kl.de,
8
9NOTE: This library uses the external programs surf, graphviz and imagemagick.
10@*    Input data is assumed to originate from resolve.lib and reszeta.lib
11
12PROCEDURES:
13InterDiv(M[,name])   dual graph of resolution of a surface (uses graphviz,imagemagick)
14ResTree(L,M[,name])  tree of charts of resolution (uses graphviz,imagemagick)
15finalCharts(L,...)   pictures of final charts of surface (uses surf)
16";
17
18proc InterDiv(intmat M, list #)
19"USAGE:  InterDiv(M[,name]);
20@*       M = matrix
21@*       name = string
22ASSUME:  - M is first list entry of output of 'intersectionDiv'
23@*         from library reszeta.lib
24@*       - write permission in the current directory or in the
25@*         directory in which the file with name 'name' resides
26CREATE:  file 'name.jpg' containing dual graph of resolution
27@*       if filename is given
28NOTE:    only available on UNIX-type systems and programs
29@*       'display' (imagemagick package) and 'dot' (Graphviz package) need to
30@*       be in the standard search PATH
31RETURN:  nothing, only generating graphics output in separate window
32EXAMPLE: not available (for technical reasons)
33"
34{
35  "Warning: alpha testing version of this procedure";
36  int i,j;
37  string tempstr;
38//---------------------------------------------------------------------------
39//  build up temporary filename and open file for writing
40//---------------------------------------------------------------------------
41  if(size(#)>0)
42  {
43     if((typeof(#[1])=="string") && (goodFilename(#[1])))
44     {
45        string @filename=#[1];
46     }
47     else
48     {
49        ERROR("Second argument should specify a WRITABLE file");
50     }
51  }
52  if(!defined(@filename))
53  {
54     string @filename=buildFilename("InterDiv.dot");
55  }
56  link eing=":w "+@filename;
57//--------------------------------------------------------------------------
58// write input for external program dot to file
59//--------------------------------------------------------------------------
60  write(eing,"graph G{");
61  for(i=1;i<=ncols(M);i++)
62  {
63     if(M[i,i]!=-2)
64     {
65        tempstr=string(i)+"[shape=circle,label=\""+string(M[i,i])+"\"];";
66     }
67     else
68     {
69        tempstr=string(i)+"[shape=point,label=\" \"];";
70     }
71     write(eing,tempstr);
72  }
73  for(i=1;i<=nrows(M);i++)
74  {
75     for(j=i+1;j<=ncols(M);j++)
76     {
77        if(M[i,j]!=0)
78        {
79           tempstr=string(i) + "--" + string(j) + ";";
80           write(eing,tempstr);
81        }
82     }
83  }
84  write(eing,"}");
85  close(eing);
86//---------------------------------------------------------------------------
87// produce graphics output using the programs dot and display
88//---------------------------------------------------------------------------
89  string outfile=@filename + ".jpg";
90  if(!find(outfile,"/"))
91  {
92//--- display needs fully qualified path to file
93     outfile=system("getenv","PWD") + "/" + outfile;
94  }
95  j=system("sh","dot -Tjpg " + @filename + " -o "+ outfile);
96  j=system("sh","display " + outfile + " &");
97//---------------------------------------------------------------------------
98// clean up if necessary
99//---------------------------------------------------------------------------
100  "Currently showing graphics in separate window";
101  "Press <Return> to continue";
102  pause();
103  if((size(#)==0)&&(find(@filename,"/tmp/")))
104  {
105//--- do not leave any garbage in the public directories
106    j=system("sh","command rm " + outfile + " " + @filename);
107  }
108  return();
109}
110///////////////////////////////////////////////////////////////////////////////
111//static
112proc goodFilename(string datei)
113{
114//--- check whether the specified file datei is writable for us
115  if(!system("sh","touch " + datei + " > /dev/null 2>&1"))
116  {
117     return(1);
118  }
119  else
120  {
121     return(0);
122  }
123}
124//////////////////////////////////////////////////////////////////////////////
125//static
126proc buildFilename(string datei)
127{
128  if(goodFilename(datei))
129  {
130     return(datei);
131  }
132  if(find(datei,"/"))
133  {
134     ERROR("Specified directory/file is not writable");
135  }
136  datei="/tmp/" + datei;
137  if(goodFilename(datei))
138  {
139     return(datei);
140  }
141//--- not reached
142  ERROR("At least /tmp should be writable");
143}
144///////////////////////////////////////////////////////////////////////////////
145
146proc ResTree(list re, intmat DivMat, list #)
147"USAGE:  ResTree(L,M[,name][,mark]);
148@*       L = list
149@*       M = matrix
150@*       name = string
151@*       mark = intvec
152ASSUME:  - L is the output of 'resolve' from resolve.lib
153@*       - M is first entry of output of 'collectDiv(L);' from reszeta.lib
154@*       - write permission in the current directory or in the
155@*         directory in which the file with name 'name' resides
156@*       - mark intvec of size size(L[2])
157@*            mark[i]=0 (default) border of box black
158@*            mark[i]=1 border of box red
159CREATE:  file 'name.jpg' containing the tree of charts of L
160@*       if filename is given
161NOTE:    only available on UNIX-type systems and programs
162@*       'display' (imagemagick package) and 'dot' (Graphviz package) need to
163@*       be in the standard search PATH
164RETURN:  nothing, only generating graphics output in separate window
165EXAMPLE: not available (for technical reasons)
166"
167{
168//-----------------------------------------------------------------------------
169// Initialization and definition of the temporary filename
170//-----------------------------------------------------------------------------
171  int i,j,dimC,jsave;
172  string tempstr;
173  def R=basering;
174  if(size(#)>0)
175  {
176     if(typeof(#[1])=="string")
177     {
178        if(goodFilename(#[1]))
179        {
180           string @filename=#[1];
181        }
182        else
183        {
184           ERROR("optional argument of type string "+
185                  "should specify a writable file.");
186        }
187     }
188     if(typeof(#[1])=="intvec")
189     {
190        intvec @rot=#[1];
191     }
192  }
193  if(size(#)>1)
194  {
195     if((typeof(#[2])=="string")&&(!defined(@filename)))
196     {
197        if(goodFilename(#[1]))
198        {
199           string @filename=#[1];
200        }
201        else
202        {
203           ERROR("optional argument of type string "+
204                  "should specify a writable file.");
205        }
206     }
207     if((typeof(#[2])=="intvec")&&(!defined(@rot)))
208     {
209        intvec @rot=#[2];
210     }
211  }
212  if(!defined(@filename))
213  {
214     string @filename=buildFilename("ResTree.dot");
215  }
216  if(!defined(@rot))
217  {
218     intvec @rot;
219     @rot[size(re[2])]=0;
220  }
221  link eing=":w "+@filename;
222//----------------------------------------------------------------------------
223// writing the input to the program dot into a file
224//----------------------------------------------------------------------------
225  write(eing,"graph G{");
226  tempstr="1[shape=box,label=\"chart 1\"];";
227  write(eing,tempstr);
228  for(i=2;i<=size(re[2]);i++)
229  {
230     tempstr=string(i)+"[shape=box,label=\"chart " + string(i)
231            + "\\nE:"+string(simplify(ideal(DivMat[i,1..ncols(DivMat)]),2))
232            + " \"";
233     if(@rot[i]==1)
234     {
235        tempstr=tempstr + "color=\"red\"];";
236     }
237     else
238     {
239        tempstr=tempstr + "];";
240     }
241     write(eing,tempstr);
242  }
243  for(i=2;i<=size(re[2]);i++)
244  {
245     def S=re[2][i];
246     setring S;
247     j=int(leadcoef(path[1,ncols(path)]));
248     if(j!=jsave)
249     {
250        def T=re[2][j];
251        setring T;
252        dimC=dim(std(BO[1]+cent));
253        setring S;
254        kill T;
255     }
256     setring R;
257     kill S;
258     if(j!=jsave)
259     {
260        tempstr=string(j) + "--" + string(i) +"[label=\"d=" + string(dimC)
261               + "\"];";
262        jsave=j;
263     }
264     else
265     {
266       tempstr=string(j) + "--" + string(i) +";";
267     }
268     write(eing,tempstr);
269  }
270  write(eing,"}");
271  close(eing);
272//---------------------------------------------------------------------------
273// Create the graphics output using the programs dot and display
274//---------------------------------------------------------------------------
275  string outfile=@filename + ".jpg";
276  if(!find(outfile,"/"))
277  {
278//--- display needs fully qualified path to file
279     outfile=system("getenv","PWD") + "/" + outfile;
280  }
281  j=system("sh", "dot -Tjpg " + @filename + " -o "+ outfile);
282  j=system("sh","display " + outfile + "&");
283//---------------------------------------------------------------------------
284// Clean up public directories if necessary
285//---------------------------------------------------------------------------
286  "Currently showing graphics in separate window";
287  "Press <Return> to continue";
288  pause();
289  if(find(@filename,"/tmp/"))
290  {
291//--- do not leave any garbage in the public directories
292    j=system("sh","command rm " + @filename + ".jpg "+ @filename);
293  }
294  return();
295}
296/////////////////////////////////////////////////////////////////////////////
297
298proc finalCharts(list re, list inter, intvec endiv, list #)
299"USAGE:  finalCharts(L1,L2,iv[,name]);
300@*       L1 = list
301@*       L2 = list
302@*       iv = intvec
303@*       name = string
304ASSUME:  - L1 is the output of 'resolve' from resolve.lib
305@*       - L2 is the output of 'intersectionDiv(L1)' from reszeta.lib
306@*       - iv is the first entry of the output of 'abstractR(L1)'
307@*       - write permission in the current directory or in the
308@*         directory in which the file with name 'name' resides
309CREATE:  - new windows in which surf-images of the final charts are presented
310@*       - several '.ras' files in the directory in which 'name' resides
311NOTE:    only available on UNIX-type systems
312@*       external programs 'surf' and 'display' (imagemagick package) need to be
313@*       in the standard search PATH
314RETURN:  nothing, only generating graphics output in separate window
315EXAMPLE: not available (for technical reasons)
316"
317{
318//-----------------------------------------------------------------------------
319// Initialization and Sanity Checks
320//-----------------------------------------------------------------------------
321  int i,j,k,a,b,cnt,whichE,offset,haveDcE8;
322  def R=basering;
323  list endCharts;
324  string fname,tempstr;
325  for(i=1;i<=size(endiv);i++)
326  {
327    if(endiv[i]==1)
328    {
329       endCharts[size(endCharts)+1]=i;
330//--- Sanity checks for the chart i
331       if(nvars(re[2][i])!=3)
332       {
333          ERROR("This chart is not embedded in 3-dimensional space");
334       }
335       if(defined(S)){kill S;}
336       def S=re[2][i];
337       setring S;
338       if(dim(std(BO[1]+BO[2]))!=2)
339       {
340          ERROR("Strict Transform is not a surface");
341       }
342       if(size(equidim(BO[1]+BO[2]))!=1)
343       {
344          ERROR("Strict Transform has lower dimensional components.");
345       }
346//--- put the missing information into dcE, if necessary
347       for(k=1;k<=size(dcE);k++)
348       {
349          for(j=1;j<=size(dcE[k]);j++)
350          {
351             if(deg(std(dcE[k][j][1])[1])!=0)
352             {
353                if(size(dcE[k][j][1]) < 8)
354                {
355                   setring R;
356                   k=prepareDcE(re,inter,endiv);
357                   setring S;
358                }
359                haveDcE8=1;
360                break;
361             }
362          }
363          if(haveDcE8!=0) break;
364       }
365       setring R;
366    }
367  }
368  if(!defined(colorlist))
369  {
370     list colorlist;
371     for(i=1;i<=10;i++)
372     {
373       colorlist[i]="curve_red = "+string((i-1)*25)+"; curve_green = "+
374                     string(255-(i-1)*25)+"; curve_blue = "+
375                     string(((i-1) mod 5)*50)+";";
376     }
377  }
378  else
379  {
380     "Warning!";
381     "Using colors specified in variable colorlist. No syntax checks";
382     "performed on the content of this list.";
383     "If built-in colors should be used, please rename the global";
384     "object colorlist";
385  }
386  if(ncols(inter[1])>size(colorlist))
387  {
388    ERROR("Too many exceptional curves on this surface.....");
389  }
390  if((size(endCharts)>20)&&(size(#)==0))
391  {
392    ERROR("More than 20 final charts...");
393  }
394//-----------------------------------------------------------------------------
395// Determine the basename of the temporary files
396//-----------------------------------------------------------------------------
397  if(size(#)>0)
398  {
399     if((typeof(#[1])=="string") && (goodFilename(#[1])))
400     {
401        string fnamebase=#[1];
402     }
403     else
404     {
405        ERROR("Second argument should specify a WRITABLE file");
406     }
407  }
408  if(!defined(@fnamebase))
409  {
410     string fnamebase=buildFilename("Chart");
411  }
412//-----------------------------------------------------------------------------
413// Go through all final charts and write a separate surf input file for each
414//-----------------------------------------------------------------------------
415  for(i=1;i<=size(endCharts);i++)
416  {
417     fname=fnamebase + string(endCharts[i]) + ".surf";
418     if(defined(eing)){kill eing;}
419     link eing=":w "+fname;
420//--- define surf's root finding algorithm
421     tempstr="root_finder = d_chain_bisection; epsilon = 0.000000001;";
422     write(eing,tempstr);
423//--- define image size
424     tempstr="int wid = 480;width = wid; height = wid;";
425     write(eing,tempstr);
426//--- define ambient lighting
427     tempstr="ambient=50;";
428     write(eing,tempstr);
429//--- define background colour: some shade of gray
430     tempstr="background_red=200; background_green=200; background_blue=200;";
431     write(eing,tempstr);
432//--- define surface colour (one side):
433     tempstr="surface_red=22; surface_green=150; surface_blue=255;";
434     write(eing,tempstr);
435//--- define surface colour (other side):
436     tempstr="inside_red=255; inside_green=192; inside_blue=0;";
437     write(eing,tempstr);
438//--- define scaling factor
439     tempstr="scale_x=0.5; scale_y=0.5; scale_z=0.5;";
440     write(eing,tempstr);
441//--- rotate a little bit
442     tempstr="rot_y=0.9013941717697173; rot_x=-0.5556596516994916; rot_z=0.103062920202253447;";
443     write(eing,tempstr);
444//----------------------------------------------------------------------------
445// change to the chart and extract the equation of the surface
446// then draw the surface
447//----------------------------------------------------------------------------
448     if(defined(S)){kill S;}
449     def S=re[2][endCharts[i]];
450     setring S;
451//--- define equation of strict transform (hypersurface of dim 2)
452     ideal drawJ=simplify(mstd(radical(BO[1]+BO[2]))[2],2);
453     if(size(drawJ)>1)
454     {
455//!!! unnoetig! Da kann man auch surface2,... machen
456        ERROR("Did not find generator of principal ideal " + string(drawJ));
457     }
458     tempstr="poly f = " + map2Surf(string(drawJ[1])) + ";";
459     write(eing,tempstr);
460     tempstr="surface = f;draw_surface;";
461     write(eing,tempstr);
462//----------------------------------------------------------------------------
463// now consider all exceptional curves in this chart separately
464// and draw them
465//---------------------------------------------------------------------------
466     if(!defined(dcE))
467     {
468//--- Oups, exceptional curves not yet decomposed in this chart
469//--- should not happen in practice
470        ERROR("The procedure intersectionDiv has not been used
471               on this tree of charts");
472     }
473     for(j=1;j<=size(dcE);j++)
474     {
475//--- Run through all exceptional divisors (in the embedded sense)
476        for(k=1;k<=size(dcE[j]);k++)
477        {
478//--- Run through all curves corresponding to the current divisor
479           if(deg(std(dcE[j][k][1])[1])==0)
480           {
481//--- This one is empty - skip it
482              k++;
483              if(k>size(dcE[j])) {break;}
484              continue;
485           }
486           for(a=1;a<=size(inter[3]);a++)
487           {
488//--- this curve belongs to which (Q-irred.) divisor in global numbering?
489              if(inIVList(intvec(endCharts[i],j,k),inter[3][a])) break;
490           }
491           if(a>size(inter[3]))
492           {
493//--- curve not found in list
494              ERROR("Inconsistency between data of arguments 1 and 2");
495           }
496           whichE=a;
497           offset=0;
498           for(a=1;a<whichE;a++)
499           {
500//--- count C-irred. divisors up to this (Q-irred. divisor) index
501              offset=offset+inter[4][a];
502           }
503           for(a=1;a<=size(dcE[j][k][7]);a++)
504           {
505//--- run through all C-irred divisors corresponding to this Q-irred. one
506//--- run through all generators of the ideal of the curve
507//--- unfortunately the ideal is encoded in a string,
508//--- and we need to work to find the generators
509              int tempint=find(dcE[j][k][6],",");
510              b=0;
511              cnt=1;
512              while(tempint!=0)
513              {
514                if(!find(dcE[j][k][7][a],"i"))
515                {
516//--- 1) surf can only handle real numbers ==> we ignore the other curves
517//--- 2) we can only form substrings of named strings:
518                  tempstr=dcE[j][k][6];            // give it a name
519                  tempstr=tempstr[b+1..tempint-1]; // find correct substring
520                  tempstr="poly f" + string(offset+dcE[j][k][8][a])+ "_" +
521                       string(b+1) + " = " +
522                       map2Surf(plugInNumZero(tempstr,dcE[j][k][7][a])) + ";";
523                  write(eing,tempstr);
524                  tempstr="cutsurface" + string(cnt) + " = f" +
525                       string(offset+dcE[j][k][8][a]) + "_" +string(b+1) + ";";
526                  write(eing,tempstr);
527                  cnt++;
528                }
529                b=tempint;                          // save end-mark
530                tempint=find(dcE[j][k][6],",",b+1); // next one please
531              }
532              if(!find(dcE[j][k][7][a],"i"))
533              {
534                tempstr=dcE[j][k][6];                 // as before, but for last
535                tempstr=tempstr[b+1..size(tempstr)];  // fragment
536                tempstr="poly f" + string(offset+dcE[j][k][8][a])+ "_" +
537                       string(b+1) + " = " +
538                       map2Surf(plugInNumZero(tempstr,dcE[j][k][7][a])) + ";";
539                write(eing,tempstr);
540                tempstr="cutsurface" + string(cnt)+ " = f" +
541                       string(offset+dcE[j][k][8][a]) + "_" +string(b+1) + ";";
542                write(eing,tempstr);
543                cnt++;
544              }
545//--- draw the curve on the surface
546              if(cnt>1)
547              {
548                 tempstr=colorlist[offset+dcE[j][k][8][a]];
549                 write(eing,tempstr);
550                 write(eing,"cut_with_surface;");
551              }
552              kill tempint;
553           }
554        }
555     }
556     if(!find(fnamebase,"/"))
557     {
558//--- display needs fully qualified path to file
559         if(defined(outfile)) {kill outfile;}
560         string outfile=system("getenv","PWD") + "/" + fnamebase;
561     }
562     if(!defined(outfile))
563     {
564       tempstr="filename = \"" + fnamebase + string(endCharts[i]) + ".ras\";";
565     }
566     else
567     {
568       tempstr="filename = \"" + outfile + string(endCharts[i]) + ".ras\";";
569     }
570     write(eing,tempstr);
571     tempstr="color_file_format = sun;";
572     write(eing,tempstr);
573     tempstr="save_color_image;";
574     write(eing,tempstr);
575     close(eing);
576     j=system("sh", "surf -n " + fnamebase +string(endCharts[i]) + ".surf");
577     if(!defined(outfile))
578     {
579        j=system("sh","display " + fnamebase
580                       + string(endCharts[i]) + ".ras &");
581     }
582     else
583     {
584        j=system("sh","display " + outfile
585                       + string(endCharts[i]) + ".ras &");
586        kill outfile;
587     }
588     "Currently showing graphics for chart "+string(endCharts[i])
589                    +" in separate window";
590     "Press <Return> to continue";
591     pause();
592  }
593  return(0);
594}
595//////////////////////////////////////////////////////////////////////////////
596// static
597proc plugInNumZero(string f,string num)
598{
599  int i;
600  string fStr=f;
601  i=find(fStr,"t");
602  while(i!=0)
603  {
604    fStr = string(fStr[1..i-1]) + "(" + num + ")"
605         + string(fStr[i+1..size(fStr)]);
606    i=find(fStr,"t");
607  }
608  return(fStr);
609}
610//////////////////////////////////////////////////////////////////////////////
611// static
612proc replaceInStr(string alles, string alt, string neu)
613{
614  int i=find(alles,alt);
615  while(i!=0)
616  {
617    if((i-1)<1)
618    {
619       if(size(alt)==size(alles))
620       {
621          alles=neu;
622       }
623       else
624       {
625          alles=neu + string(alles[i+size(alt)..size(alles)]);
626       }
627    }
628    else
629    {
630       if(i+size(alt)>size(alles))
631       {
632          alles=string(alles[1..i-1]) + neu;
633       }
634       else
635       {
636          alles=string(alles[1..i-1]) + neu +
637                string(alles[i+size(alt)..size(alles)]);
638       }
639    }
640    i=find(alles,alt);
641  }
642  return(alles);
643}
644/////////////////////////////////////////////////////////////////////////////
645// static
646proc map2Surf(string str)
647{
648  str=replaceInStr(str,string(var(1)),"x");
649  str=replaceInStr(str,string(var(2)),"y");
650  str=replaceInStr(str,string(var(3)),"z");
651  return(str);
652}
653/////////////////////////////////////////////////////////////////////////////
654// static
655proc prepareDcE(list re, list inter, list endCharts)
656{
657   def R=basering;
658//---Test whether we are in the irreducible surface case
659   def S=re[2][1];
660   setring S;
661   BO[2]=BO[2]+BO[1];              // make sure we are living in the smooth W
662   if(dim(std(BO[2]))!=2)
663   {
664     ERROR("The given original object is not a surface");
665   }
666   if(dim(std(slocus(BO[2])))>0)
667   {
668     ERROR("The given original object has non-isolated singularities.");
669   }
670   setring R;
671//----------------------------------------------------------------------------
672// Compute a non-embedded resolution from the given embedded one by
673// dropping redundant trailing blow-ups
674//----------------------------------------------------------------------------
675//--- compute non-embedded resolution
676   int ii,j,k,a,b,comPa;
677   list abst=abstractR(re);
678   intvec endiv=abst[1];
679   intvec deleted=abst[2];
680//--- identify the divisors in the various final charts
681   list iden0=collectDiv(re,deleted)[2];
682                                    // list of final divisors
683//----------------------------------------------------------------------------
684// For each Q-irred. divisor (specified in inter[3]), run through all
685// occurrences and identify the different C-components
686//----------------------------------------------------------------------------
687   for(ii=1;ii<=size(inter[3]);ii++)
688   {
689//--- run through all Q-irred. curves,
690//--- set up the first occurrence for reference
691      if(defined(S)) {kill S;}
692      def S=re[2][inter[3][ii][1][1]];
693                                    // inter[3] = list of Q-irred div.
694                                    // inter[3][ii] = ii-th thereof
695                                    // inter[3][ii][1] = first occurrence
696                                    // inter[3][ii][1][1] = corr. chart index
697      setring S;
698      dcE[inter[3][ii][1][2]][inter[3][ii][1][3]][8]=
699              intvec(1..ncols(dcE[inter[3][ii][1][2]][inter[3][ii][1][3]][4]));
700//--- prepare the ideal of the divisor for mapping to different chart
701      if(defined(idlist1)){kill idlist1;}
702      list idlist1;
703      idlist1[1]=dcE[inter[3][ii][1][2]][inter[3][ii][1][3]][6];
704      exportto(Top,idlist1);
705      for(j=2;j<=size(inter[3][ii]);j++)
706      {
707//--- now do the comparison with the other occurrences
708//--- 1. find a common parent for inter[3][ii][1][1] and inter[3][ii][j][1]
709         if(defined(S)){kill S;}
710         def S=re[2][inter[3][ii][j][1]];
711         setring S;
712         if(defined(opath)){kill opath;}
713         def opath=imap(re[2][inter[3][ii][1][1]],path);
714         k=1;
715         while(opath[1,k]==path[1,k])
716         {
717            k++;
718            if((k>ncols(opath))||(k>ncols(path))) break;
719         }
720         comPa=int(leadcoef(opath[1,k-1]));
721//--- 2. use fetchInTree to transfer the C-components
722         if(defined(str)) {kill str;}
723         if(defined(il)) {kill il;}
724         if(defined(mpi)) {kill mpi;}
725         if(defined(nulli)) {kill nulli;}
726         string str="idlist1";
727         attrib(str,"algext",imap(re[2][inter[3][ii][1][1]],dcE)[inter[3][ii][1][2]][inter[3][ii][1][3]][5]);
728         list il=fetchInTree(re,inter[3][ii][1][1],comPa,inter[3][ii][j][1],str,iden0,1);
729         list nulli=imap(re[2][inter[3][ii][1][1]],dcE)[inter[3][ii][1][2]][inter[3][ii][1][3]][7];
730         string mpi=imap(re[2][inter[3][ii][1][1]],dcE)[inter[3][ii][1][2]][inter[3][ii][1][3]][5];
731         if(mpi!=dcE[inter[3][ii][j][2]][inter[3][ii][j][3]][5])
732         {
733            ERROR("Problem identifying the appropriate field extension.!");
734         }
735         if(defined(ringt)){kill ringt;}
736         ring ringt=0,(t),dp;
737         if(defined(St)){kill St;}
738         def St=S+ringt;
739         setring St;
740         if(defined(strId)) {kill strId;}
741         string strId="ideal id1=" + il[1] + ";";
742         execute(strId);
743         id1=std(id1);
744         strId="ideal idj="
745                         + imap(S,dcE)[inter[3][ii][j][2]][inter[3][ii][j][3]][6]
746                         +  ";";
747         execute(strId);
748         idj=std(idj);
749         if(defined(nullj)){kill nullj;}
750         list nullj=imap(S,dcE)[inter[3][ii][1][2]][inter[3][ii][1][3]][7];
751         if(defined(rcomp)){kill rcomp;}
752         strId="ring rcomp=complex,(" +varstr(basering) +"),(" +
753                           ordstr(basering) + ");";
754         execute(strId);
755         def id1=imap(St,id1);
756         def idj=imap(St,idj);
757         ideal id10,idj0;
758         if(defined(tintvec)){kill tintvec;}
759         intvec tintvec;
760         tintvec[size(nullj)]=0;
761         for(a=1;a<=size(nullj);a++)
762         {
763            if(defined(numa)) {kill numa;}
764            strId="number numa=" + string(nullj[a]) + ";";
765            execute(strId);
766            idj0=subst(idj,t,numa);
767            for(b=1;b<=size(nulli);b++)
768            {
769               if(defined(numb)) {kill numb;}
770               strId="number numb=" + string(nulli[b]) + ";";
771               execute(strId);
772               id10=subst(id1,t,numb);
773               attrib(id10,"isSB",1);
774               if(size(reduce(idj0,id10))==0)
775               {
776                  tintvec[a]=b;
777                  break;
778               }
779            }
780            if(!find(string(tintvec),string(b)))
781            {
782~;
783               ERROR("Problem identifying C-components in different charts.");
784            }
785         }
786         setring S;
787         dcE[inter[3][ii][1][2]][inter[3][ii][1][3]][8]=tintvec;
788      }
789   }
790   return(0);
791}
Note: See TracBrowser for help on using the repository browser.