source: git/Singular/LIB/surfex/Singular/surfex.lib @ 3de2ca

spielwiese
Last change on this file since 3de2ca was 3de2ca, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: surfex_0_90_00 git-svn-id: file:///usr/local/Singular/svn/trunk@11070 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 41.9 KB
Line 
1//last change: 2007/07/06 (Oliver Labs)
2///////////////////////////////////////////////////////////////////////////////
3version="$Id: surfex.lib,v 1.1 2008-09-23 11:12:53 Singular Exp $";
4category="Visualization";
5info="
6LIBRARY: surfex.lib   
7         Procedures for Visualizing Surfaces.
8         This is still an alpha version! Please, send me any comments
9         or suggestions via http://www.AlgebraicSurface.net
10         Thank you!!!
11AUTHOR: Oliver Labs
12        This library uses the program surfex.
13        surfex is written by Oliver Labs and others,
14        mainly Stephan Holzer.
15        It is based on some other tools, mainly
16        surf which was written by Stefan Endrass and others.
17        Some of the code is based on the code in surf.lib
18
19NOTE:
20 Many procedures of this library require the program surfex
21 to be installed.
22 This software is used for producing raytraced images of the surfaces.
23 You can download @code{surfex} from @uref{http://www.surfex.AlgebraicSurface.net}.
24 
25 surfex is a front-end for surf which aims to be easier to use than
26 the original tool. In the future, we plan to extend it to a tool which
27 can use other appropriate programs for producing the best visualization.
28 Please, contact me via @uref{http://www.surfex.AlgebraicSurface.net}
29 if you have any questions or suggestions,
30
31 Oliver Labs
32 
33PROCEDURES:
34proc plotRotated(poly p, list coords, list #)
35               Plot the surface given by the polynomial p
36               with the coordinates coords.
37proc plotRot(poly p, list #)
38               Similar to plotRotated, but tries to guess automatically
39               which coordinates should be used.
40proc plotRotatedList(list varieties, list coords, list #);
41               Plot the varieties given by the list varieties
42               with the coordinates coords.
43proc plotRotatedDirect(list varietiesList, list #);
44               Plot the varieties given by the list varietiesList.
45               varietiesList has a complicated format.
46               The procedure passes the data directly to surfex.
47               Thus, the coordinates have to be x,y,z.
48               Parameters are allowed; their names have to be p1, p2, ...
49proc plotRotatedListFromSpecifyList(list varietiesList, list #) ;
50               Plot the varieties given by the list varietiesList.
51               varietiesList has a complicated format; see the example.
52";
53
54LIB "solve.lib";
55LIB "primdec.lib";
56LIB "sing.lib";
57
58///////////////////////////////////////////////////////////
59//
60// the main procedures:
61//
62
63proc plotRot(poly p, list #)
64"
65USAGE: plotRot(poly p, list #)
66Similar to plotRotated, but tries to guess automatically
67which coordinates should be used.
68
69It opens the external program surfex for drawing the surface given by p,
70seen as a surface in the real affine space with coordinates coords.
71
72RETURN: the return code of the system command which executes surfex.
73
74ASSUME: The basering is characteristic zero and without parameters.
75"
76{
77    list coords = list();
78    if(num_vars_id(p)==3) {
79        execute("coords = "+string_of_vars(p)+";");
80    } else {
81        if(num_vars_id(p)<3) {
82            if(nvars(basering)==3) {
83                execute("coords = "+varstr(basering)+";");
84            } else {
85                if(nvars(basering)<3) {
86                    "Could not guess the coordinates because the number of variables in the basering is smaller than 3!";
87                    "Please use plotRotated() instead of plotRot() and specify the coordinates explicitly.";
88                    return(0);
89                } else {
90                    "Could not guess the coordinates because the number of variables in the polynomial is smaller than 3 and the number of variables in the basering is greater than three!";
91                    "Please use plotRotated() instead of plotRot() and specify the coordinates explicitly.";
92                    return(0);
93                }
94            }
95        } else {
96            "Could not guess the coordinates because the number of variables in the polynomial is greater than 3!";
97                    "Please use plotRotated() instead of plotRot() and specify the coordinates explicitly.";
98            return(0);
99        }
100    }
101    return(plotRotatedList(list(p), coords, #));
102}
103example {
104    "Example:"; echo=2;
105
106    // More variables in the basering, but only 3 variables in the polynomial:
107    ring r1 = 0, (w,x,y,z), dp;
108    poly cayley_cubic = x^3+y^3+z^3+1^3-1/4*(x+y+z+1)^3;
109    plotRot(cayley_cubic);
110
111    // Three variables in the basering, but fewer variables in the polynomial:
112    ring r2 = 0, (x,y,z), dp;
113    plotRot(x^2+y^2-1);
114    plotRot(y^2+z^2-1);
115
116    // A cubic surface with a solitary point:
117    // Use the additional parameter 3 to ask singular
118    // to compute the singular locus before calling surfex.
119    ring r3 = 0, (x,y,z), dp; 
120    poly kn_10 = x^3-3*x*y^2+z^3+3*x^2+3*y^2+z^2;
121    plotRot(kn_10, 3);
122
123    // The swallowtail:
124    // a surface with a real solitary curve sticking out of the surface.
125    // Use the additional parameter 3 to ask singular
126    // to compute the singular locus before calling surfex.
127    poly swallowtail = -4*y^2*z^3-16*x*z^4+27*y^4+144*x*y^2*z+128*x^2*z^2-256*x^3;
128    plotRot(swallowtail, 3);
129}
130
131proc plotRotated(poly p, list coords, list #)
132"
133USAGE: plotRotated(poly p, list coords, list #)
134This opens the external program surfex for drawing the surface given by p,
135seen as a surface in the real affine space with coordinates coords.
136
137RETURN: the return code of the system command which executes surfex.
138
139ASSUME: coords is a list of three variables.
140The basering is characteristic zero and without parameters.
141"
142{
143    return(plotRotatedList(list(p), coords, #));
144}
145example {
146    "Example:"; echo=2;
147
148    // An easy example: a surface with four conical nodes.
149    ring r = 0, (x,y,z), dp;
150    poly cayley_cubic = x^3+y^3+z^3+1^3-1/4*(x+y+z+1)^3;
151//     plotRotated(cayley_cubic, list(x,y,z));
152
153    // An difficult example: a surface with a one-dimensional real component!
154    poly whitney_umbrella = x^2*z-y^2;
155    // The Whitney Umbrella without its handle:
156    plotRotated(whitney_umbrella, list(x,y,z));
157    // The Whitney Umbrella together with its handle:
158    plotRotated(whitney_umbrella, list(x,y,z), 2);
159}
160
161
162proc plotRotatedList(list varieties, list coords, list #)
163"
164USAGE: plotRotatedList(list varieties, list coords, list #)
165This opens the external program surfex for drawing the surfaces given by varieties,
166seen as a surface in the real affine space with coordinates coords.
167
168RETURN: the return code of the system command which executes surfex.
169
170ASSUME: coords is a list of three variables,
171varieties is a list of ideals discribing the varieties to be shown.
172The basering is characteristic zero and without parameters.
173"
174{
175    def oring = basering;
176
177    int plotquality = 0;
178    if(size(#)>0) {
179        plotquality = #[1];
180    }
181
182    list varietiesList = list(list(), list(), list(), list());
183    list usedSurfaces = list();
184    list curveColors = list();
185   
186    // go through the list of varieties
187    // produce a list which can be used as input for plotRotatedListFromList()
188    int i;
189    int j;
190    list indList;
191    int ind;
192    ideal itmp;
193    int ncurves;
194    list pd;
195    int k;
196    int surfind;
197    list curSurfColors = list();
198
199    list listOfPoints = list();
200    string str_I = "";
201   
202    for(i=1; i<=size(varieties); i++) {
203        itmp = varieties[i];
204        if(plotquality>=3) {
205            itmp = radical(itmp);
206        }
207        itmp = simplify(itmp,1);
208        itmp = simplify(itmp,2);
209        if(size(itmp)==1) { // i.e.: a surface given by one equation
210            surfind = findInList(surfEqn(itmp[1],coords), usedSurfaces);
211            if(surfind==0) {
212                usedSurfaces = usedSurfaces + list(surfEqn(itmp[1],coords));
213                curSurfColors = list(list("insidecolor:",getInsideColorStr(size(varietiesList[1])+1)),
214                                     list("outsidecolor:",getOutsideColorStr(size(varietiesList[1])+1)));
215                varietiesList[1] = varietiesList[1] +
216                    list(list(list("eqno:",string(size(varietiesList[1])+1)),
217                              list("equation:",surfEqn(itmp[1], coords)),
218                              curSurfColors[1],
219                              curSurfColors[2],
220                              list("showcbox:","true"),
221                              list("transparency:","0")));
222                surfind = size(varietiesList[1]);
223                             
224            }
225            if(plotquality==1) {
226                varieties = varieties + list(slocus(itmp[1]));
227            }
228            if(plotquality==2 || plotquality==3) {
229                // remove doubled components and
230                // add the 1-dimensional singular components
231                // of the surface to the list of curves:
232                int dsl = dim_slocus(itmp[1]);
233                dsl;
234                if(dsl>=0) { // i.e. there is a singular locus
235                    "compute singular locus...";
236                    list eqd;
237                    //
238                    eqd = equidim(slocus(itmp[1]));
239                    ideal tmp_l;
240                    tmp_l = std(eqd[size(eqd)]);
241                    "dim:",dim(tmp_l);
242                    if(dim(tmp_l)==(nvars(basering)-3+2)) {
243                        "--- 2-dim.";
244                        // we have found a multiple component;
245                        // replace it by a simple copy of it
246                        itmp = quotient(itmp[1], tmp_l);
247                        varieties[i] = itmp[1];
248                        eqd = delete(eqd,size(eqd));
249                        if(size(eqd)>0) {
250                            tmp_l = std(eqd[size(eqd)]);
251                        }
252                    }
253                    if(dim(tmp_l)==(nvars(basering)-3+1)) {
254                        "--- 1-dim.";
255                        // we have found a 1-dimensional singular locus
256                        pd = std_primdecGTZ(tmp_l,2);
257                        for(k=1; k<=size(pd); k++) {
258                            if(pd[k][3]==(nvars(basering)-3+1)) {
259                                varieties = varieties + list(pd[k][2]);
260                                curveColors[size(varieties)] = curSurfColors;
261                            } else {
262                                "???";
263                            }
264                        }
265                        eqd = delete(eqd,size(eqd));
266                        if(size(eqd)>0) {
267                            tmp_l = std(eqd[size(eqd)]);
268                        }
269                    }
270                    if(dim(tmp_l)==(nvars(basering)-3+0)) {
271                        "--- 0-dim.";
272                        // we have found a 0-dimensional singular locus
273                        // we compute floating point approximations of the
274                        // coordinates of all singular points
275                        if(npars(oring)>0) {
276                            "str:",parstr(1),rootminpoly();
277                            list all_real_sols = allroots_minpoly();
278//                          "all sols:";all_real_sols;
279//                          sprintf("number %s = %s; ", parstr(1), rootminpoly());
280                            int minp;
281                            if((npars(basering) == 1) && (minpoly != 0)) {
282                                minp = 1;
283                            } else {
284                                minp = 0;
285                            }
286                            str_I = "";
287                            if(minp==1) {
288                                "minp=1";
289                                string str_para = parstr(1);
290                                string str_tmp_l;
291                                def cur_ring = basering;
292                                if(1) {
293                                    short=0;
294                                    str_tmp_l = "ideal eqd_tmp = "+
295//                                      string(tmp_l)+","+string(minpoly)+";";
296                                        string(tmp_l);
297                                    "str:",str_tmp_l;
298                                    string str_num_mp = "number "+parstr(1)+"="+
299                                        decstr2ratstr(rootminpoly())+";";
300                                    execute("ring Iring = 0,("
301//                                          +string(coords)+","+str_para+"),dp;");
302                                            +string(coords)+"),dp;");
303                                    basering;
304                                    execute(str_num_mp);
305                                    execute(str_tmp_l);
306                                    eqd_tmp;
307                                    list real_sols = real_solve(eqd_tmp);
308                                    real_sols;
309                                    $;
310                                    setring cur_ring;
311                                }
312                            } else {
313                                // minp==0: we do not know how to handle this
314                                "???";
315                            }
316                        } else {
317                            "no pars";
318                            ideal eqd_tmp = tmp_l;
319                            short=0;
320                            string str_tmp_l = "ideal eqd_tmp = "+string(tmp_l)+";";
321                            def cur_ring = basering;
322                            execute("ring Iring = (real,30),("+string(coords)+"),("+ordstr(oring)+");");
323//                          basering;
324                            execute(str_I);
325                            execute(str_tmp_l);
326                            list real_sols = real_solve(eqd_tmp);
327                            setring cur_ring;
328                        }
329                        "real_sols:";real_sols;
330                        for(k=1; k<=size(real_sols); k++) {
331                            "search point:";
332                            string(list(real_sols[k][1],real_sols[k][2],real_sols[k][3],string(surfind)));
333//                          listOfPoints;
334                            if(findInList(string(list(list(real_sols[k][1],real_sols[k][2],real_sols[k][3],string(surfind)))),
335                                          listOfPoints)==0) {
336                                "add pt";
337                                varietiesList[4] = varietiesList[4] +
338                                    list(list(real_sols[k][1],real_sols[k][2],real_sols[k][3],string(surfind)));
339                                listOfPoints = listOfPoints +
340                                    list(string(list(real_sols[k][1],real_sols[k][2],real_sols[k][3],string(surfind))));
341                            }
342                        }
343                    }
344                }
345            }
346        } else {
347            // i.e.: more than one equation
348            varietiesList[2] = varietiesList[2] +
349                list(list(list("surfaces:"),
350                          list("curveno:",
351                               string(size(varietiesList[2])+1)),
352                          list("showcbox:","true")));
353            if(size(curveColors) >= i) {
354                varietiesList[2][size(varietiesList[2])][4] = curveColors[i][1];
355                varietiesList[2][size(varietiesList[2])][4][1] = "color:";
356            }
357            ncurves = size(varietiesList[2]);
358            for(j=1; j<=size(itmp); j++) {
359                ind = findInList(surfEqn(itmp[j],coords), usedSurfaces);
360                usedSurfaces = usedSurfaces + list(surfEqn(itmp[1],coords));
361//              "indList:";indList;
362                if(ind == 0) {
363//                  "--------> not in list", surfEqn(itmp[j], coords);
364                    if(j==1) {
365                        varietiesList[1] = varietiesList[1] +
366                            list(list(list("eqno:",string(size(varietiesList[1])+1)),
367                                      list("equation:",surfEqn(itmp[j], coords)),
368                                      list("insidecolor:",getInsideColorStr(size(varietiesList[1])+1)),
369                                      list("outsidecolor:",getOutsideColorStr(size(varietiesList[1])+1)),
370                                      list("showcbox:","true"),
371                                      list("transparency:","100")));
372                    } else {
373                        varietiesList[1] = varietiesList[1] +
374                            list(list(list("eqno:",string(size(varietiesList[1])+1)),
375                                      list("equation:",surfEqn(itmp[j], coords)),
376                                      list("insidecolor:",getInsideColorStr(size(varietiesList[1])+1)),
377                                      list("outsidecolor:",getOutsideColorStr(size(varietiesList[1])+1)),
378                                      list("showcbox:","false"),
379                                      list("transparency:","0")));
380                    }
381                    ind = size(varietiesList[1]);
382                } else {
383                }
384                varietiesList[2][ncurves][1] = varietiesList[2][ncurves][1] + list(string(ind));
385            }
386        }
387    }
388
389//      "------------";
390//      varietiesList;
391//      "------------";
392    return(plotRotatedListFromSpecifyList(varietiesList, coords, #));
393}
394example {
395    "Example:"; echo=2;
396
397    // A cubic surface together with a tritangent plane
398    // (i.e. a plane which cuts out three lines).
399    ring r = 0, (x,y,z), dp;
400    poly cayley_cubic = x^3+y^3+z^3+1^3-1/4*(x+y+z+1)^3;
401    poly plane = 1-x-y-z;
402    plotRotatedList(list(cayley_cubic, plane), list(x,y,z));
403
404    // The same cubic and plane.
405    // The plane is not shown but only its intersection with the surface.
406    plotRotatedList(list(cayley_cubic, ideal(cayley_cubic, plane)), list(x,y,z));
407}
408
409
410proc plotRotatedListFromSpecifyList(list varietiesList, list #)
411"
412USAGE: plotRotatedListFromSpecifyList(list varietiesList, list #);
413varietiesList has a complicated format (not documented yet);
414see the example.
415
416RETURN: the return code of the system command which executes surfex.
417
418ASSUME: The basering is characteristic zero.
419
420EXAMPLE: example plotRotatedListFromSpecifyList;
421"
422{
423    // make the surfex file
424    string str = getSurfexCodeFromSpecifyList(varietiesList, #);
425
426    return(plotRotatedFromCode(str, #));
427}
428example {
429    "Example:"; echo=2;
430   
431    // A cubic surface depending on a parameter:
432    ring r = (0,p1), (x,y,z), dp;
433    poly cayley_cubic = x^3+y^3+z^3+1^3-p1*(x+y+z+1)^3;
434    poly plane = 1-x-y-z;
435    plotRotatedListFromSpecifyList(list(list(list(list("eqno:","1"),
436                                                  list("equation:",
437                                                       string(cayley_cubic))
438                                                 )
439                                            ),
440                                        list(),
441                                        list(list(1,"0.0","1.0","500","0.25+0.25*sin(PI*p1)")),
442                                        list()
443                                       ));
444}
445
446
447proc plotRotatedListFromStringList(list varieties, list #)
448"
449RETURN: the return code of the system command which executes surfex.
450
451USAGE: not documented yet.
452"
453{
454    // make the surfex file
455    getSurfexCodeFromStringList(varieties, #);
456    string str = getSurfexCodeFromStringList(varieties, #);
457   
458    return(plotRotatedFromCode(str, #));
459}
460
461
462proc plotRotatedDirect(list varieties, list #)
463"
464USAGE: plotRotatedDirect(list varieties, list #)
465This opens the external program surfex for drawing the surfaces given by varieties,
466seen as a surface in the real affine space with coordinates x,y,z.
467The format for the list varieties is not fully documented yet;
468please, try to adapt the examples to your needs.
469
470RETURN: the return code of the system command which executes surfex.
471
472ASSUME:
473Passes the equations directly to surfex,
474i.e., the variable names should be correct, i.e. x,y,z.
475The advantage is that one can use parameters p1, p2, ...;
476these will be passed to surfex.
477"
478{
479    string str = getSurfexCodeFromListDirect(varieties, #);
480   
481    return(plotRotatedFromCode(str, #));
482}
483example {
484    "Example:"; echo=2;
485   
486    // A cubic surface depending on a parameter:
487    ring r = (0,p1), (x,y,z), dp;
488    poly cayley_cubic = x^3+y^3+z^3+1^3-p1*(x+y+z+1)^3;
489    // The entries of the list of varieties can either be polynomials
490    plotRotatedDirect(list(list(list(cayley_cubic)),
491                           list(),
492                           list(list(1,"0.0","1.0","500","0.25+0.25*sin(PI*p1)"))
493                           ));
494
495    // or strings which represent surfex-readable polynomials
496    plotRotatedDirect(list(list(list("x^3+y^3+z^3+1^3-p1*(x+y+z+1)^3")),
497                           list(),
498                           list(list("1","0.0","1.0","500","0.25+0.25*sin(PI*p1)"))
499                           ));
500
501    // More complicated varieties
502    plotRotatedDirect(list(list(list("x^2+y^2-z^2-3^2"),
503                                list("x*sin(p1)+y*cos(p1)-3")),
504                           list(list(list(1,2))),
505                           list(list("1","0.0","1.0","500","2*PI*p1"))
506                           ));
507}
508
509proc plotRotatedFromCode(string str, list #)
510"
511USAGE: plotRotatedFromCode(string str, list #);
512
513This procedure is only for internal usage;
514it takes the surfex-code as a string and calls surfex.
515
516RETURN: the return code of the system command which executes surfex.
517"
518{
519     // we need a temporary .sux file for surfex
520    string tmpd = "/tmp";
521    string l="surf"+string(system("pid"))+".sux";
522    // a temporary file which stores the output of surfex
523    string erg="/tmp/surferg"+string(system("pid"));
524
525    write(":w "+tmpd+"/"+l, str);
526
527    int i=system("sh","surfex -d "+tmpd+" -i " + l +" >"+erg+" 2>/dev/null");
528
529    // delete the temporary file
530    i = system("sh","rm " + l +" 2>/dev/null");   
531    return(read(erg));
532}
533
534
535///////////////////////////////////////////////////////////
536//
537// procedures used to produce the surf-code:
538//
539
540
541proc getSurfexCodeFromListDirect(list varieties, list #)
542"
543USAGE: getSurfexCodeFromListDirect(list varieties, list #)
544
545ASSUME: varieties has four components,
546        - the first is a list of polynomials, say f_1, ..., f_k
547        - the second is a list of lists of numbers in {1, ..., k} describing the curves
548          as intersections of the corresponding f_i
549        - the third is a list of lists describing the parameters used in the polynomials f_i
550        - the fourth is a list of lists of points given by their approximate coordinates (three decimal numbers)
551
552RETURN: the surfex code (.sux)
553"
554{
555    int i;
556    int j;
557    string str = "this is surfex v0.89.07"+newline;
558
559    str = str + "TYPE:" + newline;
560    str = str + "specify"+newline;
561    str = str + "EQUATIONS:"+newline;
562    str = str + string(size(varieties[1])) + newline;
563    for(i=1; i<=size(varieties[1]); i++) {
564        str = str + "Equation:"+newline;
565        str = str + "eqno:"+newline;
566        str = str + string(i) + newline;
567        str = str + "equation:"+newline;
568        str = str + surfEqnDir(varieties[1][i][1]) + newline;
569        if(size(varieties[1][i])>=2) {
570            str = str + "showcbox:"+newline;
571            str = str + varieties[1][i][2] + newline;     // show it or not
572            if(size(varieties[1][i])>=3) {
573                str = str + "transparency:"+newline;
574                str = str + string(varieties[1][i][3]) + newline;     // transparency
575            }
576        }
577    }
578    str = str + "CURVES:"+newline;
579    str = str + string(size(varieties[2])) + newline;
580    for(i=1; i<=size(varieties[2]); i++) {
581        str = str + "Curve:"+newline;
582        str = str + "curveno:"+newline;
583        str = str + string(i) + newline;
584        str = str + "surfaces:"+newline;
585//      "curves:";varieties[2][i];
586        for(j=1; j<=size(varieties[2][i][1]); j++) {
587            str = str + string(varieties[2][i][1][j]) + newline;
588        }
589        if(size(varieties[2][i])>=2) {
590            str = str + "showcbox:"+newline;
591            str = str + varieties[2][i][2] + newline;     // show it or not
592        }
593    }
594    str = str + "PARAMETERS:"+newline;
595    str = str + string(size(varieties[3])) + newline;
596    for(i=1; i<=size(varieties[3]); i++) {
597        str = str + "Parameter:"+newline;
598        str = str + "parno:"+newline;
599        str = str + string(varieties[3][i][1]) + newline;
600        str = str + "fromtoval:"+newline;
601        str = str + varieties[3][i][2] + newline;
602        str = str + varieties[3][i][3] + newline;
603        str = str + string(varieties[3][i][4]) + newline;
604        if(size(varieties[3][i])>=5) {
605            str = str + "function:"+newline;
606            str = str + varieties[3][i][5]+newline;
607        }
608    }   
609//     str = str + "////////////////// Parameter: /////////////////////////"+newline;
610//     str = str + "1" + newline;
611//     str = str + "0.0" + newline;
612//     str = str + "1.0" + newline;
613//     str = str + "1000" + newline;
614//    str = str + string(size(varieties[3])) + newline;
615    return(str);
616}
617
618proc getSurfexCodeFromList(list varieties, list coords, list #)
619"
620ASSUME: varieties has four components,
621        - the first is a list of polynomials, say f_1, ..., f_k
622        - the second is a list of lists of numbers in {1, ..., k} describing the curves
623          as intersections of the corresponding f_i
624        - the third is a list of lists describing the parameters used in the polynomials f_i
625        - the fourth is a list of lists of points given by their approximate coordinates (three decimal numbers)
626
627RETURN: the surfex code (.sux)
628"
629{
630    int i;
631    int j;
632    string str = "this is surfex v0.89.07"+newline;
633
634    str = str + "TYPE:" + newline;
635    str = str + "specify"+newline;
636    str = str + "EQUATIONS:"+newline;
637    str = str + string(size(varieties[1])) + newline;
638    for(i=1; i<=size(varieties[1]); i++) {
639        str = str + "Equation:"+newline;
640        str = str + "eqno:"+newline;
641        str = str + string(i) + newline;
642        str = str + "equation:"+newline;
643        str = str + surfEqn(varieties[1][i][1], coords) + newline;
644        str = str + "showcbox:"+newline;
645        str = str + varieties[1][i][2] + newline;     // show it or not
646        str = str + "transparency:"+newline;
647        str = str + string(varieties[1][i][3]) + newline;     // transparency
648    }
649    str = str + "CURVES:"+newline;
650    str = str + string(size(varieties[2])) + newline;
651    for(i=1; i<=size(varieties[2]); i++) {
652        str = str + "Curve:"+newline;
653        str = str + "curveno:"+newline;
654        str = str + string(i) + newline;
655        str = str + "surfaces:"+newline;
656        for(j=1; j<=size(varieties[2][i]); j++) {
657            str = str + string(varieties[2][i][1][j]) + newline;
658        }
659        str = str + "showcbox:"+newline;
660        str = str + varieties[2][i][2] + newline;     // show it or not
661    }
662    str = str + "PARAMETERS:"+newline;
663    str = str + string(size(varieties[3])) + newline;
664    for(i=1; i<=size(varieties[3]); i++) {
665        str = str + "Parameter:"+newline;
666        str = str + "parno:"+newline;
667        str = str + string(varieties[3][i][1]) + newline;
668        str = str + "fromtoval:"+newline;
669        str = str + surfEqn(varieties[3][i][2], coords) + newline;
670        str = str + surfEqn(varieties[3][i][3], coords) + newline;
671        str = str + string(varieties[3][i][4]) + newline;
672        if(size(varieties[3][i])>=5) {
673            str = str + "function:"+newline;
674            str = str + varieties[3][i][5]+newline;
675        }
676    }   
677//     str = str + "////////////////// Parameter: /////////////////////////"+newline;
678//     str = str + "1" + newline;
679//     str = str + "0.0" + newline;
680//     str = str + "1.0" + newline;
681//     str = str + "1000" + newline;
682//    str = str + string(size(varieties[3])) + newline;
683    return(str);
684}
685
686proc getSurfexCodeFromStringList(list varieties, list #)
687"
688ASSUME: varieties has three components,
689        - the first is a list of polynomials, say f_1, ..., f_k
690        - the second is a list of lists of numbers in {1, ..., k} describing the curves
691          as intersections of the corresponding f_i
692        - the third is a list of lists describing the parameters used in the polynomials f_i
693
694RETURN: the surfex code (.sux)
695"
696{
697    int i;
698    int j;
699    string str = "this is surfex v0.89.07"+newline;
700
701    str = str + "TYPE:" + newline;
702    str = str + "specify"+newline;
703    str = str + "EQUATIONS:"+newline;
704    str = str + string(size(varieties[1])) + newline;
705    for(i=1; i<=size(varieties[1]); i++) {
706        str = str + "Equation:"+newline;
707        str = str + "eqno:"+newline;
708        str = str + string(i) + newline;
709        str = str + "equation:"+newline;
710        str = str + varieties[1][i][1] + newline;
711        str = str + "showcbox:"+newline;
712        str = str + varieties[1][i][2] + newline;     // show it or not
713        str = str + "transparency:"+newline;
714        str = str + varieties[1][i][3] + newline;     // transparency
715    }
716    str = str + "CURVES:"+newline;
717    str = str + string(size(varieties[2])) + newline;
718    for(i=1; i<=size(varieties[2]); i++) {
719        str = str + "Curve:"+newline;
720        str = str + "curveno:"+newline;
721        str = str + string(i) + newline;
722        str = str + "surfaces:"+newline;
723        for(j=1; j<=size(varieties[2][i][1]); j++) {
724            str = str + string(varieties[2][i][1][j]) + newline;
725        }
726        str = str + "showcbox:"+newline;
727        str = str + varieties[2][i][2] + newline;     // show it or not
728    }
729    str = str + "PARAMETERS:"+newline;
730    str = str + string(size(varieties[3])) + newline;
731    for(i=1; i<=size(varieties[3]); i++) {
732        str = str + "Parameter:"+newline;
733        str = str + "parno:"+newline;
734        str = str + string(varieties[3][i][1]) + newline;
735        str = str + "fromtoval:"+newline;
736        str = str + varieties[3][i][2] + newline;
737        str = str + varieties[3][i][3] + newline;
738        str = str + string(varieties[3][i][4]) + newline;
739        if(size(varieties[3][i])>=5) {
740            str = str + "function:"+newline;
741            str = str + varieties[3][i][5]+newline;
742        }
743    }   
744    return(str);
745}
746
747
748proc getSurfexCodeFromSpecifyList(list varieties, list #)
749"
750ASSUME: varieties has three components,
751        - the first is a list of polynomials, say f_1, ..., f_k
752        - the second is a list of lists of numbers in {1, ..., k} describing the curves
753          as intersections of the corresponding f_i
754        - the third is a list of lists describing the parameters used in the polynomials f_i
755        - the fourth is a list of lists describing the singular points to be shown as spheres
756
757RETURN: the surfex code (.sux)
758"
759{
760    int i;
761    int j;
762    int k;
763    string str = "this is surfex v0.89.07"+newline;
764
765    str = str + "TYPE:" + newline;
766    str = str + "specify"+newline;
767    str = str + "EQUATIONS:"+newline;
768    str = str + string(size(varieties[1])) + newline;
769    for(i=1; i<=size(varieties[1]); i++) {
770        str = str + "Equation:"+newline;
771        for(j=1; j<=size(varieties[1][i]); j++) {
772            str = str + varieties[1][i][j][1] +newline;
773            str = str + varieties[1][i][j][2] +newline;
774        }
775    }
776    str = str + "CURVES:"+newline;
777    str = str + string(size(varieties[2])) + newline;
778    for(i=1; i<=size(varieties[2]); i++) {
779        str = str + "Curve:"+newline;
780        for(j=1; j<=size(varieties[2][i]); j++) {
781            str = str + varieties[2][i][j][1] +newline;
782            if(varieties[2][i][j][1] == "surfaces:") {
783                for(k=2; k<=size(varieties[2][i][j]); k++) {
784                    str = str + string(varieties[2][i][j][k]) + newline;
785                }
786            } else {
787                str = str + varieties[2][i][j][2] +newline;
788            }
789        }
790//      str = str + "curveno:"+newline;
791//      str = str + string(i) + newline;
792//      str = str + "surfaces:"+newline;
793//      for(j=1; j<=size(varieties[2][i][1]); j++) {
794//          str = str + string(varieties[2][i][1][j]) + newline;
795//      }
796//      str = str + "showcbox:"+newline;
797//      str = str + varieties[2][i][2] + newline;     // show it or not
798    }
799    str = str + "PARAMETERS:"+newline;
800    str = str + string(size(varieties[3])) + newline;
801    for(i=1; i<=size(varieties[3]); i++) {
802        str = str + "Parameter:"+newline;
803        str = str + "parno:"+newline;
804        str = str + string(varieties[3][i][1]) + newline;
805        str = str + "fromtoval:"+newline;
806        str = str + varieties[3][i][2] + newline;
807        str = str + varieties[3][i][3] + newline;
808        str = str + string(varieties[3][i][4]) + newline;
809        if(size(varieties[3][i])>=5) {
810            str = str + "function:"+newline;
811            str = str + varieties[3][i][5]+newline;
812        }
813    }   
814    string str_from = "0.0";
815    string str_to = "5.0";
816    string str_radius = "50";
817    str = str + "SOLITARY POINTS:"+newline;
818    str = str + string(size(varieties[4])) + newline;
819    for(i=1; i<=size(varieties[4]); i++) {
820        str = str + "SolitaryPoint:"+newline;
821        str = str + "solPtNo:"+newline;
822        str = str + string(i) + newline;
823        str = str + "surface:"+newline;
824        str = str + varieties[4][i][4] + newline;
825        str = str + "fromtoval:"+newline;
826        str = str + str_from + newline;
827        str = str + str_to + newline;
828        str = str + str_radius + newline;
829        str = str + "coords:" + newline;
830        str = str + varieties[4][i][1] + newline;
831        str = str + varieties[4][i][2] + newline;
832        str = str + varieties[4][i][3] + newline;
833    }   
834    return(str);
835}
836
837///////////////////////////////////////////////////////////
838//
839// procedures for standard colors:
840//
841
842proc numBaseColors()
843"
844USAGE: numBaseColors()
845
846RETURN: the number of predefined surface colors.
847"
848{
849    return(6);
850}
851
852proc baseSurfaceColors(int no)
853"
854USAGE: baseSurfaceColors(int no)
855
856REMARK: There are currently 6=numBaseColors() basic surface colors.
857You can modify them according to your wishes
858by just redefining this procedure in your Singular-script.
859
860If you want more colors, then you also have to redefine numBaseColors() accordingly.
861
862RETURN: a list of three integers describing the RGB values of a color.
863"
864{
865    if(no%numBaseColors()==1) {
866        return(list(240,160,0));
867    }
868    if(no%numBaseColors()==2) {
869        return(list(160,240,0));
870    }
871    if(no%numBaseColors()==3) {
872        return(list(0,160,240));
873    }
874    if(no%numBaseColors()==4) {
875        return(list(240,0,160));
876    }
877    if(no%numBaseColors()==5) {
878        return(list(0,240,160));
879    }
880    if(no%numBaseColors()==0) {
881        return(list(160,0,240));
882    }
883}
884
885proc getInsideColorStr(int no)
886"
887USAGE: getInsideColorStr(int no)
888
889RETURN: a string describing inside color number no
890where the three integer RGB values are in one line each.
891"
892{
893    list bc = baseSurfaceColors(no);
894    string str = string(bc[1])+newline+string(bc[2])+newline+string(bc[3]);
895    return(str);
896}
897
898proc getOutsideColorStr(int no)
899"
900USAGE: getOutsideColorStr(int no)
901
902RETURN: a string describing outside color number no
903where the three integer RGB values are in one line each.
904"
905{
906    list bc = baseSurfaceColors(no);
907    string str = string(bc[1])+newline+string(bc[2])+newline+string(bc[3]);
908    return(str);
909}
910
911///////////////////////////////////////////////////////////
912//
913// procedures used by the plot procedures:
914//
915
916proc surfEqnDir(list #)
917"
918USAGE: surfEqnDir(list #) without any checks etc.
919
920RETURN: string(#[1]) where short=0.
921"
922{
923    int stmp = short; short = 0;
924    string str = string(#[1]);
925    short = stmp;
926    return(str);
927}
928
929proc surfEqn(poly p, list coords, list #)
930
931USAGE: surfEqn(poly p, list coords)
932       Tries to produce a string for the equation of p which is convenient for surfex.
933ASSUME: - p defines a plane curve or a surface,
934         - coords is a list of the three coordinates to use, e.g. list(x,y,z),
935           in this way, it is possible to distinguish between x^2+y^2-1 and y^2+z^2-1
936RETURN: a string, that one can use with the external program surf
937EXAMPLE: example surfEqn; shows an example
938"
939{
940    int params=0;
941    if(size(#)>0) {
942        params = #[1];
943    }
944    string err_mes; // string containing error messages
945    def base=basering;
946    int mynvars = nvars(basering);
947
948    intvec ind=num_of_vars(p);
949
950    int i,j,n;
951    int minp = 0;
952    n=0;
953    for(i=size(ind);i>0;i--)
954    {
955        if (ind[i]!=0) {
956            n++;
957        } else {
958            if(var(i)==coords[1] || var(i)==coords[2] || var(i)==coords[3]) {
959                ind[i]=1;
960                n++;
961            }
962        }
963    }
964   
965    params = params + npars(basering);
966    n = n + npars(basering);
967    if((npars(basering) == 1) && (minpoly != 0)) {
968        minp = 1;
969    } else {
970        minp = 0;
971    }
972    string str_I = "";
973    for(i=1; i<=npars(basering); i=i+1) {
974        if(!(parstr(i) == "i")) {
975            if(minp==1) {
976                str_I = str_I + sprintf("number %s = %s; ", parstr(i), rootminpoly());
977            } else {
978            }
979        }
980    }
981    int bshort = short; short = 0;
982    if(!(minp==1 || npars(basering)==0)) {
983        p=cleardenom(p);
984        err_mes="Cannot plot equations with a parameter without a specified minpoly";
985        ERROR(err_mes);
986    }
987    str_I = str_I + "poly p = " + string(p) + ";";
988
989    short = bshort;
990
991    if(params==0) {
992        if (n<=2 or n>=4)
993        {
994            err_mes="Cannot plot equations with "+string(n)+" variables";
995            ERROR(err_mes);
996//          return("0");
997        }
998        if(n==4) {
999            ring r=(real,30,30),(xx,yy,zz,ww),dp;
1000        } else {
1001            ring r=(real,30,30),(x,y,z),dp;
1002        }
1003    } else {
1004        if(n-params<=2 || n-params>=4) {
1005            err_mes="Cannot plot equations with "+string(n-params)+" variables";
1006            ERROR(err_mes);
1007//          return("0");
1008        } else {
1009            if(params == 1) {
1010                if(n-params==3) {
1011                    if(minp==1) {
1012                        // switch to a ring without minimal polynomial:
1013                        execute("ring rr = (real,30,30),("+varstr(base)+"), dp;");
1014//                      rr;
1015//                      "str_I",str_I;
1016                        execute(str_I);
1017                        def base = rr;
1018                        ring r=(real,30,30),(x,y,z),dp;
1019                    } else {
1020                        p=cleardenom(p);
1021                        ring r=(real,30,30),(x,y,z,p1),dp;
1022                    }
1023                }
1024            }
1025            if(params == 2) {
1026                if(n-params==3) {
1027                    p=cleardenom(p);
1028                    ring r=(real,30,30),(x,y,z,p1,p2),dp;
1029                }
1030            }
1031            if(params == 3) {
1032                if(n-params==3) {
1033                    p=cleardenom(p);
1034                    execute("ring rr = (real,30,30),("+varstr(base)+","+parstr(base)+"), dp;");
1035                    rr;
1036                    "str_I",str_I;
1037                    execute(str_I);
1038                    "pnew:",p;
1039                    def base = rr;
1040
1041                    ring r=(real,30,30),(x,y,z,p1,p2,p3),dp;
1042                }
1043            }
1044        }
1045    }
1046//    basering;
1047    short=0;
1048    map phi=base,0;
1049    j=1;
1050
1051    for(i=1;i<=mynvars;i++)
1052    {
1053        if (ind[i]!=0)
1054        {
1055            phi[i]=var(j);
1056            j++;
1057        }
1058    }
1059    poly p=(simplify(phi(p),1));
1060    if (leadcoef(p) <0) {
1061        if(size(#)>1) {
1062            if(#[2]!=0) {
1063                p=-p;
1064            }
1065        } else {
1066            p=-p;
1067        }
1068    }
1069    if(leadcoef(p)!=0) {
1070        p = p/leadcoef(p);
1071    }
1072    string thesurfstr = string(p);
1073    if(minp == 1) {
1074        // replace k by rootRepl
1075    }
1076
1077    return (thesurfstr);
1078} // end of surfEqn()
1079example
1080{ "EXAMPLE:"; echo =2;
1081
1082  ring rr0 = 0,(x(1..3)),dp;
1083  poly p = x(1)^3 - x(2)^2;
1084  print(surfEqn(p,list(x(1),x(2),x(3))));
1085
1086  ring rr1 = 0,(x,y,z),dp;
1087  poly I(1) = 2x2-1/2x3 +1-y+1;
1088  print(surfEqn(I(1),list(x,y,z)));
1089
1090  // Steiner surface
1091  poly J(2) = x^2*y^2+x^2*z^2+y^2*z^2-17*x*y*z;
1092  print(surfEqn(J(2),list(x,y,z)));
1093} // end of example surfEqn()
1094
1095
1096proc num_of_vars(ideal I)
1097"
1098USAGE: num_of_vars(ideal I)
1099
1100RETURN: an intvec containing one entry for each ring variable.
1101each contains the sums of all degrees in this variable of all monomials
1102occuring in the ideal.
1103An entry is zero iff the corresponding variable does not occur in the ideal.
1104"
1105{
1106  intvec v;
1107  int i;
1108  poly p;
1109  for(i=size(I);i>0;i--)
1110  {
1111    p=I[i];
1112    while(p!=0)
1113    {
1114        v=v+leadexp(p);
1115        p=p-lead(p);
1116    }
1117  }
1118  return(v);
1119}
1120example {
1121    "EXAMPLE:"; echo = 2;
1122    ring r = 0, (x,y,z),dp;
1123    ideal j0 = x^2-x*y;
1124    num_of_vars(j0);
1125    ideal j1 = x^2-x*y-y;
1126    num_of_vars(j1);
1127    ideal j2 = x^2-x*y-y, x^3-2*y;
1128    num_of_vars(j2);
1129}
1130
1131
1132proc num_vars_id(ideal I)
1133"
1134USAGE: num_vars_id(ideal I)
1135
1136RETURN: The number of ring-variables occurring in the ideal I.
1137"
1138{
1139    intvec v = num_of_vars(I);
1140    int num = 0;
1141    for(int i=size(v);i>0;i--)
1142    {
1143        if (v[i]!=0) { num++; }
1144    }
1145    return(num);
1146}
1147example {
1148    "EXAMPLE:"; echo = 2;
1149    ring r = 0, (x,y,z),dp;
1150    ideal j = x^2-y, x^3-2;
1151    num_vars_id(j);
1152}
1153
1154proc findInList(list obj, list l)
1155"
1156USAGE: findInList(list obj, list l)
1157       Tries to find the object obj in the list l.
1158
1159ASSUME: the object obj[1] can be compared to the objects in the list l
1160
1161RETURN: if obj[1]=l[i] for some i, then return the first such i,
1162        otherwise return 0
1163"
1164{
1165    for(int i=1; i<=size(l); i++) {
1166        if(l[i]==obj[1]) {
1167            return(i);
1168        }
1169    }
1170
1171    return(0);
1172}
1173example {
1174    "EXAMPLE:"; echo = 2;
1175    ring r = 0,(x,y,z), dp;
1176    list a = list(x^2+y^2+z^2+1, x^2+y^2+z^2-1, x^2+y^2-z^2+1, x^2+y^2-z^2-1);
1177    findInList(x^2+y^2+z^2-1, a);
1178    findInList(x^2+y^2+z^2, a);
1179}
1180
1181proc std_primdecGTZ(ideal I, list #)
1182"
1183USAGE: std_primdecGTZ(ideal I, list #)
1184Computes a primdary decomposition pd of I using primdecGTZ and then
1185calls std_for_pd(pd).
1186For the output and options, consult the help of std_for_pd.
1187
1188RETURN: see std_for_pd.
1189"
1190{
1191    list pd = primdecGTZ(I);
1192    return(std_for_pd(pd, #));
1193}
1194example {
1195    "EXAMPLE:"; echo = 2;
1196
1197    ring r = 0, (x,y), dp;
1198    ideal j = y-x^2,z-x^3;
1199    primdecGTZ(j);
1200    std_primdecGTZ(j);
1201    std_primdecGTZ(j,1);
1202}
1203
1204proc std_for_pd(list pd, list #)
1205"
1206USAGE: std_for_pd(list pd, list #)
1207Call std for each of the prime ideals in the list pd
1208replace the prime ideals by their standard-basis.
1209Compute dim() and mult() of each prime component using these standard bases.
1210If an additional argument is given then do the same for the primary components.
1211
1212ASSUME:
1213pd is in the format produced by primdecGTZ() or primdecSY().
1214
1215RETURN: A list, say l, of lists, similar to a list returned by primdecSY() or primdecGTZ().
1216However, each of the entries of l (which is a list l[i]) contains some additional entries:
1217l[1]: the primary ideal
1218l[2]: a standard basis of the associated prime ideal
1219l[3]: dim() of this prime ideal
1220l[4]: mult() of this prime ideal
1221
1222If an additional argument # is given then l[1] changes:
1223l[1]: a standard basis of the primary ideal
1224Morever, there are some more entries:
1225l[5]: dim() of this primary ideal
1226l[6]: mult() of this primary ideal
1227l[7]: l[6] / l[5]
1228"
1229{
1230
1231    if(typeof(pd[1])=="ideal") {
1232        // this is a Singular bug!?
1233//      "bug!";pd;"---";
1234        pd = list(list(pd[1], pd[1]));
1235//      pd;$;
1236    }
1237    list pd_neu;
1238    int i;
1239    list coords;
1240    ideal stdtmp;
1241    ideal stdtmp2;
1242    for(i=1; i<=size(pd); i++) {
1243        stdtmp = std(pd[i][2]);
1244        stdtmp2 = pd[i][1];
1245        if(size(#)>0) {
1246            stdtmp2 = std(stdtmp2);       
1247            if(mult(stdtmp)==0) {
1248                pd_neu[i] = list(stdtmp2,
1249                                 stdtmp,
1250                                 dim(stdtmp), mult(stdtmp),
1251                                 dim(stdtmp2), mult(stdtmp2),
1252                                 0);
1253            } else {
1254                pd_neu[i] = list(stdtmp2,
1255                                 stdtmp,
1256                                 dim(stdtmp), mult(stdtmp),
1257                                 dim(stdtmp2), mult(stdtmp2),
1258                                 mult(stdtmp2)/mult(stdtmp));
1259            }
1260        } else {
1261            pd_neu[i] = list(stdtmp2,
1262                             stdtmp,
1263                             dim(stdtmp), mult(stdtmp));
1264        }
1265    }
1266    return(pd_neu);
1267}
1268example {
1269    "EXAMPLE:"; echo = 2;
1270
1271    ring r = 0, (x,y,z), dp;
1272    ideal j = y-x^2,z-x^3;
1273    list pd = primdecGTZ(j);
1274    pd;
1275    std_for_pd(pd, 1);
1276}
1277
1278proc real_solve(ideal to_solve)
1279"
1280USAGE: real_solve(ideal to_solve)
1281
1282RETURN: a list of all real solutions (as strings)
1283of the zero-dimensional ideal to_solve (without multiplicities).
1284
1285REMARK: Until now, it may happen that some points appear more than once.
1286"
1287{
1288    int k;
1289    int i;
1290
1291//    def Isolring = solve(to_solve,30,0,60,"nodisplay");
1292    def Isolring = solve(to_solve,9,0,13,"nodisplay");
1293    setring Isolring;
1294//    list SOL = solve(to_solve, "oldring", "nodisplay");
1295    list real_sols = list();
1296    list tmpl;
1297    for(k=1; k<=size(SOL); k++) {
1298        if(find(string(SOL[k]),"I")==0 && find(string(SOL[k]),"i")==0) {           
1299            tmpl = list();
1300            for(i=1; i<=size(SOL[k]); i++) {
1301                tmpl = tmpl + list(string(SOL[k][i]));
1302            }
1303            real_sols = real_sols + list(tmpl);
1304        }
1305    }   
1306    return(real_sols);
1307}
1308example {
1309    "EXAMPLE:"; echo = 2;
1310    ring r = 0, (x,y), dp;
1311    number a = 2;
1312    number b = 3;
1313    ideal j = (x^2-a),(y^3-b);
1314    real_solve(j);
1315}
1316
1317proc rootminpoly(list #)
1318"
1319USAGE: rootminpoly(list #)
1320
1321RETURN: A root of the current minpoly
1322as a string representation of a complex number with
1323the given precision #[1] (default: 30).
1324E.g. ring r=(0,s),x,dp; minpoly = s^2-2; => rootminpoly() 1.41421356237309504880168872421
1325
1326ASSUME: The current minpoly is non-zero.
1327"
1328{
1329    int prec = 30;
1330    int k, done;
1331    if(size(#)>0) {
1332        prec = #[1];
1333    }
1334    short = 0;
1335    string str_lag = sprintf("list lag = laguerre_solve(%s);", minpoly);
1336    string str_ring = sprintf("ring r_sqrt = (complex,prec,I),(%s),lp;", parstr(basering));
1337    execute(str_ring);
1338    execute(str_lag);
1339//    lag;
1340    // choose a real solution, if it exists:
1341    done = 0;
1342    for(k=1; k<=size(lag) && done==0; k++) {
1343        if(find(string(lag[k]),"I")==0) {           
1344            done = k;
1345        }
1346    }     
1347    if(done==0) {
1348//      "no real solution.";
1349    }
1350
1351    if(size(lag)>2) {
1352        // return the first real solution
1353        return(sprintf("%s",lag[done]));
1354    }
1355
1356    if(sprintf("%s",lag[1])[1] == "-") {
1357        return(sprintf("%s",lag[2]));
1358    } else {
1359        if(sprintf("%s",lag[1])[1] == "(") {
1360            if(sprintf("%s",lag[1])[2] == "-") {
1361                return(sprintf("%s",lag[2]));
1362            } else {
1363                return(sprintf("%s",lag[1]));
1364            }
1365        } else {
1366            return(sprintf("%s",lag[1]));
1367        }
1368    }
1369    short = 1;
1370}
1371example
1372{
1373   "EXAMPLE:"; echo =2;
1374   ring r=(0,s),x,dp;
1375   minpoly = s^2-2;
1376   rootminpoly();
1377
1378   ring R=(0,s),x,dp;
1379   minpoly = s^2+2;
1380   rootminpoly();
1381}
1382
1383proc allroots_minpoly(list #)
1384"
1385USAGE: allroots_minpoly(list #)
1386
1387RETURN: a list of strings containing all real roots of the minimal polynomial of the active ring.
1388
1389ASSUME: The current minpoly is non-zero.
1390"
1391{
1392    int prec = 30;
1393    int k, done;
1394    if(size(#)>0) {
1395        prec = #[1];
1396    }
1397    short = 0;
1398    string str_lag = sprintf("list lag = laguerre_solve(%s);", minpoly);
1399    string str_ring = sprintf("ring r_sqrt = (complex,prec,I),(%s),lp;", parstr(basering));
1400    execute(str_ring);
1401    execute(str_lag);
1402
1403    // only take the real solutions:
1404    done = 0;
1405    list real_sols = list();
1406    for(k=1; k<=size(lag) && done==0; k++) {
1407        if(find(string(lag[k]),"I")==0) {           
1408            real_sols = real_sols + list(string(lag[k]));
1409        }
1410    }     
1411    return(real_sols);
1412}
1413example {
1414    "EXAMPLE:"; echo = 2;
1415   ring r=(0,s),x,dp;
1416   minpoly = s^3-2;
1417   allroots_minpoly();
1418
1419   ring R=(0,s),x,dp;
1420   minpoly = s^2-2;
1421   allroots_minpoly();
1422}
1423
1424proc decstr2ratstr(string str)
1425"
1426USAGE: decstr2ratstr(string str)
1427Convert a decimal number of not more than 30 digits to a rational number with 14 digits.
1428
1429REMARK: This procedure still has to be adapted to accept other precisions!
1430"
1431{
1432    ring decR = (complex,30,I),(x),lp;
1433    execute("number r="+str+";");
1434    execute("r = "+truncdec(r,14)+";");
1435    return(real2ratstr(r));
1436}
1437
1438proc real2ratstr(number r)
1439"
1440USAGE: real2ratstr(number r)
1441
1442RETURN: A string containing a rational number representing the decimal number r.
1443
1444ASSUME: The current ring has either real or complex base field.
1445"
1446{
1447    string ratstr = "number("+string(r*number(10000000000000000))+")/number(10000000000000000)";
1448    return(ratstr);
1449}
1450
1451proc truncdec(number r, int decs)
1452"
1453USAGE: truncdec(number r, int decs)
1454Truncates a decimal number r to the given number (decs) of digits.
1455
1456RETURN: A string representing the truncated number.
1457"
1458{
1459    string str = string(r);
1460    return(str[1,(decs+2)]);
1461}
1462
1463proc string_of_vars(ideal I)
1464"
1465USAGE: string_of_vars(ideal I)
1466
1467RETURN: A string of all variables contained in the ideal I, separated by commas.
1468"
1469{
1470    list listvars = list();
1471    intvec v;
1472    int i;
1473    poly p;
1474    for(i=size(I);i>0;i--)
1475    {
1476        p=I[i];
1477        while(p!=0)
1478        {
1479            v=v+leadexp(p);
1480            p=p-lead(p);
1481        }
1482    }
1483    for(i=1; i<=nvars(basering); i++) {
1484        if(v[i] > 0) {
1485            listvars = listvars + list(var(i));
1486        }
1487    }
1488    string strvars = string(listvars);
1489    return(strvars);
1490}
Note: See TracBrowser for help on using the repository browser.