source: git/Singular/LIB/surfex.lib @ 334c21f

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