Changeset 5c67581 in git


Ignore:
Timestamp:
Nov 14, 1999, 10:35:09 PM (24 years ago)
Author:
Moritz Wenk <wenk@…>
Branches:
(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
Children:
9542e4ff27140734032740c6ee3b0997bac33185
Parents:
0c85269e6cc5c8c6b47583b9c80b51d11cf82033
Message:
wenk: changes in solve.lib


git-svn-id: file:///usr/local/Singular/svn/trunk@3820 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/solve.lib

    r0c85269 r5c67581  
    11///////////////////////////////////////////////////////////////////////////////
    22
    3 version="$Id: solve.lib,v 1.12 1999-07-08 10:18:13 wenk Exp $";
     3version="$Id: solve.lib,v 1.13 1999-11-14 21:35:02 wenk Exp $";
    44info="
    55LIBRARY: solve.lib     PROCEDURES TO SOLVE POLYNOMIAL SYSTEMS
    … …  
    77
    88PROCEDURES:
    9  ures_solve(i,..);      find all roots of 0-dimensional ideal i with resultants
    10  mp_res_mat(i,..);      multipolynomial resultant matrix of ideal i
    11  laguerre_solve(p,..);  find all roots of univariate polynom p
    12  interpolate(i,j,d);    interpolate poly from evaluation points i and results j
     9ures_solve(i,..);      find all roots of 0-dimensional ideal i with resultants
     10mp_res_mat(i,..);      multipolynomial resultant matrix of ideal i
     11laguerre_solve(p,..);  find all roots of univariate polynom p
     12interpolate(i,j,d);    interpolate poly from evaluation points i and results j
     13
     14fglm_solve(i,p,...);    find roots of 0-dim. ideal using FGLM and lex_solve
     15triangL_solve(l,p,...); find roots using triangular system (Lazard)
     16triangLf_solve(l,p,..); find roots using triangular sys. (factorizing Lazard)
     17triangM_solve(l,p,...); find roots of given triangular system (Moeller)
     18
     19lex_solve(i,p,...);    find roots of reduced lexicographic standard basis
     20triang_solve(l,p,...); find roots of given triangular system
     21
     22pcheck(i,l,...);       checks if elements (numbers) of l are roots of ideal i
    1323";
     24
     25LIB "triang.lib"; // needed for triang*_solve
    1426
    1527///////////////////////////////////////////////////////////////////////////////
    … …  
    2941"
    3042{
    31   int typ=0;
    32   int polish=2;
    33   int prec=30;
    34 
    35   if ( size(#) >= 1 )
    36   {
    37     typ= #[1];
    38     if ( typ < 0 || typ > 1 )
    39     {
    40       ERROR("Valid values for second parameter k are:
    41       0: use sparse Resultant (default)
    42       1: use Macaulay Resultant");
    43     }
    44   }
    45   if ( size(#) >= 2 )
    46   {
    47     prec= #[2];
    48     if ( prec == 0 ) { prec = 30; }
    49     if ( prec < 0 )
    50     {
    51       ERROR("Third parameter l must be positive!");
    52     }
    53   }
    54   if ( size(#) >= 3 )
    55   {
    56     polish= #[3];
    57     if ( polish < 0 || polish > 2 )
    58     {
    59       ERROR("Valid values for fourth parameter m are:
    60       0,1,2: number of iterations for approximation of roots");
    61     }
    62   }
    63 
    64   if ( size(#) > 3 )
    65   {
    66     ERROR("only three parameters allowed!");
    67   }
    68 
    69   return(uressolve(gls,typ,prec,polish));
    70 
    71 }
    72 example
    73 {
    74   "EXAMPLE:";echo=2;
    75   // compute the intersection points of two curves
    76   ring rsq = 0,(x,y),lp;
    77   ideal gls=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    78   ures_solve(gls);
    79   // result is a list (x,y)-coordinates as strings
    80 
    81   // now with complex coefficient field, precision is 20 digits
    82   ring rsc= (real,20,I),(x,y),lp;
    83   ideal i = (2+3*I)*x2 + (0.35+I*45.0e-2)*y2 - 8, x2 + xy + (42.7)*y2;
    84   list l= ures_solve(i);
    85   // result is a list of (x,y)-coordinates of complex numbers
    86   l;
    87   // check the result
    88   subst(subst(i[1],x,l[1][1]),y,l[1][2]);
     43    int typ=0;
     44    int polish=2;
     45    int prec=30;
     46
     47    if ( size(#) >= 1 )
     48    {
     49        typ= #[1];
     50        if ( typ < 0 || typ > 1 )
     51        {
     52            ERROR("Valid values for second parameter k are:
     53                   0: use sparse Resultant (default)
     54                   1: use Macaulay Resultant");
     55        }
     56    }
     57    if ( size(#) >= 2 )
     58    {
     59        prec= #[2];
     60        if ( prec == 0 ) { prec = 30; }
     61        if ( prec < 0 )
     62        {
     63            ERROR("Third parameter l must be positive!");
     64        }
     65    }
     66    if ( size(#) >= 3 )
     67    {
     68        polish= #[3];
     69        if ( polish < 0 || polish > 2 )
     70        {
     71            ERROR("Valid values for fourth parameter m are:
     72                   0,1,2: number of iterations for approximation of roots");
     73        }
     74    }
     75
     76    if ( size(#) > 3 )
     77    {
     78        ERROR("only three parameters allowed!");
     79    }
     80
     81    return(uressolve(gls,typ,prec,polish));
     82
     83}
     84example
     85{
     86    "EXAMPLE:";echo=2;
     87    // compute the intersection points of two curves
     88    ring rsq = 0,(x,y),lp;
     89    ideal gls=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     90    ures_solve(gls);
     91    // result is a list (x,y)-coordinates as strings
     92
     93    // now with complex coefficient field, precision is 20 digits
     94    ring rsc= (real,20,I),(x,y),lp;
     95    ideal i = (2+3*I)*x2 + (0.35+I*45.0e-2)*y2 - 8, x2 + xy + (42.7)*y2;
     96    list l= ures_solve(i);
     97    // result is a list of (x,y)-coordinates of complex numbers
     98    l;
     99    // check the result
     100    subst(subst(i[1],x,l[1][1]),y,l[1][2]);
    89101}
    90102///////////////////////////////////////////////////////////////////////////////
    … …  
    101113"
    102114{
    103   int polish=2;
    104   int prec=30;
    105 
    106   if ( size(#) >= 1 )
    107   {
    108     prec= #[1];
    109     if ( prec == 0 ) { prec = 30; }
    110     if ( prec < 0 )
    111     {
    112       ERROR("Fisrt parameter must be positive!");
    113     }
    114   }
    115   if ( size(#) >= 2 )
    116   {
    117     polish= #[2];
    118     if ( polish < 0 || polish > 2 )
    119     {
    120       ERROR("Valid values for third parameter are:
    121       0,1,2: number of iterations for approximation of roots");
    122     }
    123   }
    124   if ( size(#) > 2 )
    125   {
    126     ERROR("only two parameters allowed!");
    127   }
    128 
    129   return(laguerre(f,prec,polish));
    130 
    131 }
    132 example
    133 {
    134   "EXAMPLE:";echo=2;
    135   // Find all roots of an univariate polynomial using Laguerre's method:
    136   ring rs1= 0,(x,y),lp;
    137   poly f = 15x5 + x3 + x2 - 10;
    138   laguerre_solve(f);
    139 
    140   // Now with 10 digits precision:
    141   laguerre_solve(f,10);
    142 
    143   // Now with complex coefficients, precision is 20 digits:
    144   ring rsc= (real,20,I),x,lp;
    145   poly f = (15.4+I*5)*x^5 + (25.0e-2+I*2)*x^3 + x2 - 10*I;
    146   list l = laguerre_solve(f);
    147   l;
    148   // check result, value of substituted poly should be near to zero
    149   subst(f,x,l[1]);
    150   subst(f,x,l[2]);
     115    int polish=2;
     116    int prec=30;
     117
     118    if ( size(#) >= 1 )
     119    {
     120        prec= #[1];
     121        if ( prec == 0 ) { prec = 30; }
     122        if ( prec < 0 )
     123        {
     124            ERROR("Fisrt parameter must be positive!");
     125        }
     126    }
     127    if ( size(#) >= 2 )
     128    {
     129        polish= #[2];
     130        if ( polish < 0 || polish > 2 )
     131        {
     132            ERROR("Valid values for third parameter are:
     133                   0,1,2: number of iterations for approximation of roots");
     134        }
     135    }
     136    if ( size(#) > 2 )
     137    {
     138        ERROR("only two parameters allowed!");
     139    }
     140
     141    return(laguerre(f,prec,polish));
     142
     143}
     144example
     145{
     146    "EXAMPLE:";echo=2;
     147    // Find all roots of an univariate polynomial using Laguerre's method:
     148    ring rs1= 0,(x,y),lp;
     149    poly f = 15x5 + x3 + x2 - 10;
     150    laguerre_solve(f);
     151
     152    // Now with 10 digits precision:
     153    laguerre_solve(f,10);
     154
     155    // Now with complex coefficients, precision is 20 digits:
     156    ring rsc= (real,20,I),x,lp;
     157    poly f = (15.4+I*5)*x^5 + (25.0e-2+I*2)*x^3 + x2 - 10*I;
     158    list l = laguerre_solve(f);
     159    l;
     160    // check result, value of substituted poly should be near to zero
     161    subst(f,x,l[1]);
     162    subst(f,x,l[2]);
    151163}
    152164///////////////////////////////////////////////////////////////////////////////
    … …  
    162174"
    163175{
    164   int typ=2;
    165 
    166   if ( size(#) == 1 )
    167   {
    168     typ= #[1];
    169     if ( typ < 0 || typ > 1 )
    170     {
    171       ERROR("Valid values for third parameter are:
    172       0: sparse resultant (default)
    173       1: Macaulay resultant");
    174     }
    175   }
    176   if ( size(#) > 1 )
    177   {
    178     ERROR("only two parameters allowed!");
    179   }
    180 
    181   return(mpresmat(i,typ));
    182 
    183 }
    184 example
    185 {
    186   "EXAMPLE:";echo=2;
    187   // compute resultant matrix in ring with parameters (sparse resultant matrix)
    188   ring rsq= (0,u0,u1,u2),(x1,x2),lp;
    189   ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
    190   module m = mp_res_mat(i);
    191   print(m);
    192   // computing sparse resultant
    193   det(m);
    194 
    195   // compute resultant matrix (Macaulay resultant matrix)
    196   ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
    197   ideal h=  homog(imap(rsq,i),x0);
    198   h;
     176    int typ=2;
     177
     178    if ( size(#) == 1 )
     179    {
     180        typ= #[1];
     181        if ( typ < 0 || typ > 1 )
     182        {
     183            ERROR("Valid values for third parameter are:
     184                   0: sparse resultant (default)
     185                   1: Macaulay resultant");
     186        }
     187    }
     188    if ( size(#) > 1 )
     189    {
     190        ERROR("only two parameters allowed!");
     191    }
     192
     193    return(mpresmat(i,typ));
     194
     195}
     196example
     197{
     198    "EXAMPLE:";echo=2;
     199    // compute resultant matrix in ring with parameters (sparse resultant matrix)
     200    ring rsq= (0,u0,u1,u2),(x1,x2),lp;
     201    ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
     202    module m = mp_res_mat(i);
     203    print(m);
     204    // computing sparse resultant
     205    det(m);
     206
     207    // compute resultant matrix (Macaulay resultant matrix)
     208    ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
     209    ideal h=  homog(imap(rsq,i),x0);
     210    h;
    199211 
    200   module m = mp_res_mat(h,1);
    201   print(m);
    202   // computing Macaulay resultant (should be the same as above!)
    203   det(m);
    204 
    205   // compute numerical sparse resultant matrix
    206   setring rsq;
    207   ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
    208   module mn = mp_res_mat(ir);
    209   print(mn);
    210   // computing sparse resultant
    211   det(mn);
     212    module m = mp_res_mat(h,1);
     213    print(m);
     214    // computing Macaulay resultant (should be the same as above!)
     215    det(m);
     216
     217    // compute numerical sparse resultant matrix
     218    setring rsq;
     219    ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
     220    module mn = mp_res_mat(ir);
     221    print(mn);
     222    // computing sparse resultant
     223    det(mn);
    212224}
    213225///////////////////////////////////////////////////////////////////////////////
    … …  
    229241"
    230242{
    231   return(vandermonde(p,w,d));
    232 }
    233 example
    234 {
    235   "EXAMPLE:";  echo=2;
    236   ring r1 = 0,(x),lp;
    237   // determine f with deg(f) = 4 and
    238   // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
    239   ideal v=16,0,11376,1046880,85949136;
    240   interpolate( 3, v, 4 );
    241 
    242   ring r2 = 0,(x,y),dp;
    243   // determine f with deg(f) = 3 and
    244   // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15)
    245   // valuation point (2,3)
    246   ideal p = 2,3;
    247   ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
    248   poly ip= interpolate( p,v,3 );
    249   ip;
    250   // compute poly at point 2,3, result must be 2
    251   subst(subst(ip,x,2),y,3);
    252   // compute poly at point 2^15,3^15, result must be 16
    253   subst(subst(ip,x,2^15),y,3^15);
    254 }
    255 ///////////////////////////////////////////////////////////////////////////////
     243    return(vandermonde(p,w,d));
     244}
     245example
     246{
     247    "EXAMPLE:";  echo=2;
     248    ring r1 = 0,(x),lp;
     249    // determine f with deg(f) = 4 and
     250    // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
     251    ideal v=16,0,11376,1046880,85949136;
     252    interpolate( 3, v, 4 );
     253
     254    ring r2 = 0,(x,y),dp;
     255    // determine f with deg(f) = 3 and
     256    // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15)
     257    // valuation point (2,3)
     258    ideal p = 2,3;
     259    ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
     260    poly ip= interpolate( p,v,3 );
     261    ip;
     262    // compute poly at point 2,3, result must be 2
     263    subst(subst(ip,x,2),y,3);
     264    // compute poly at point 2^15,3^15, result must be 16
     265    subst(subst(ip,x,2^15),y,3^15);
     266}
     267///////////////////////////////////////////////////////////////////////////////
     268
     269static proc psubst( int d, int dd, int n, list res,
     270                    ideal fi, int elem, int nv )
     271{
     272    //   nv: number of ring variables         (fixed value)
     273    // elem: number of elements in ideal fi   (fixed value)
     274    //   fi: input ideal                      (fixed value)
     275    //   rl: output list of roots
     276    //  res: actual list of roots
     277    //    n:
     278    //   dd: actual element of fi
     279    //    d: actual variable 
     280
     281    int olddd=dd;
     282
     283    if ( pdebug>=1 )
     284    {"// 0 step "+string(dd)+" of "+string(elem);}
     285 
     286    if ( dd <= elem )
     287    {
     288        int loop = 1;
     289        int k;
     290        list lsr,lh;
     291        poly ps;
     292        int thedd;
     293
     294        if ( pdebug>=1 )
     295        {"// 1 dd = "+string(dd);}
     296
     297        thedd=0;
     298        while ( (dd+1 <= elem) && loop )
     299        {
     300            ps= fi[dd+1];
     301
     302            if ( n-1 > 0 )
     303            {
     304                if ( pdebug>=2 )
     305                {
     306                    "// 2 ps=fi["+string(dd+1)+"]"+" size="
     307                        +string(size(coeffs(ps,var(n-1))))
     308                        +"  leadexp(ps)="+string(leadexp(ps));
     309                }
     310                if ( size(coeffs(ps,var(n-1))) == 1 )
     311                {
     312                    dd++;
     313                    // hier Leading-Exponent puefen???
     314                    // oder ist das Poly immer als letztes in der Liste?!?
     315                    // leadexp(ps)
     316                }
     317                else
     318                {
     319                    loop=0;
     320                }
     321            }
     322            else
     323            {
     324                if ( pdebug>=2 )
     325                {
     326                    "// 2 ps=fi["+string(dd+1)+"]"+"  leadexp(ps)="
     327                        +string(leadexp(ps));
     328                }
     329                dd++;
     330            }
     331        }
     332        thedd=dd;
     333        ps= fi[thedd];
     334
     335        if ( pdebug>=1 )
     336        {
     337            "// 3    fi["+string(thedd-1)+"]"+"  leadexp(fi[thedd-1])="
     338                +string(leadexp(fi[thedd-1]));
     339            "// 3 ps=fi["+string(thedd)+"]"+"  leadexp(ps)="
     340                +string(leadexp(ps));
     341        }
     342
     343        for ( k= nv; k > nv-d; k-- )
     344        {
     345            if ( pdebug>=2 )
     346            {
     347                "// 4 subst(fi["+string(thedd)+"],"
     348                    +string(var(k))+","+string(res[k])+");";
     349            }
     350            ps = subst(ps,var(k),res[k]);
     351        }
     352             
     353        if ( pdebug>=2 )
     354        { "// 5 substituted ps="+string(ps); }
     355
     356        if ( ps != 0 )
     357        {
     358            lsr= laguerre_solve( ps );
     359        }
     360        else
     361        {
     362            if ( pdebug>=1 )
     363            { "// 30 ps == 0, thats not cool..."; }
     364            lsr=@ln; // lsr=number(0);
     365        }
     366
     367        if ( pdebug>=1 )
     368        { "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]"; }
     369         
     370        if ( size(lsr) > 1 )
     371        {
     372            if ( pdebug>=1 )
     373            {
     374                "// 10 checking roots found before, range "
     375                    +string(dd-olddd)+" -- "+string(dd);
     376                "// 10 thedd = "+string(thedd);
     377            }
     378             
     379            int i,j,l;
     380            int ls=size(lsr);
     381            int lss;
     382            poly pss;
     383            list nares;
     384            int rroot;
     385            int nares_size;
     386
     387             
     388            for ( i = 1; i <= ls; i++ ) // lsr[1..ls]
     389            {
     390                rroot=1;
     391
     392                if ( pdebug>=2 )
     393                {"// 13 root lsr["+string(i)+"] = "+string(lsr[i]);}
     394                for ( l = 0; l <= dd-olddd; l++ )
     395                {
     396                    if ( l+olddd != thedd )
     397                    {
     398                        if ( pdebug>=2 )
     399                        {"// 11 checking ideal element "+string(l+olddd);}
     400                        ps=fi[l+olddd];
     401                        if ( pdebug>=3 )
     402                        {"// 14 ps=fi["+string(l+olddd)+"]";}
     403                        for ( k= nv; k > nv-d; k-- )
     404                        {
     405                            if ( pdebug>=3 )
     406                            {
     407                                "// 11 subst(fi["+string(olddd+l)+"],"
     408                                    +string(var(k))+","+string(res[k])+");";
     409                            }
     410                            ps = subst(ps,var(k),res[k]);
     411                         
     412                        }
     413                       
     414                        pss=subst(ps,var(k),lsr[i]); // k=nv-d
     415                        if ( pdebug>=3 )
     416                        { "// 15 0 == "+string(pss); }
     417                        if ( pss != 0 )
     418                        {
     419                            if ( system("complexNearZero",
     420                                        leadcoef(pss),
     421                                        myCompDigits) )
     422                            {
     423                                if ( pdebug>=2 )
     424                                { "// 16 root "+string(i)+" is a real root"; }
     425                            }
     426                            else
     427                            {
     428                                if ( pdebug>=2 )
     429                                { "// 17 0 == "+string(pss); }
     430                                rroot=0;
     431                            }
     432                        }
     433                     
     434                    }
     435                }
     436               
     437                if ( rroot == 1 ) // add root to list ?
     438                {
     439                    if ( size(nares) > 0 )
     440                    {
     441                        nares=nares[1..size(nares)],lsr[i];
     442                    }
     443                    else
     444                    {
     445                        nares=lsr[i];
     446                    }
     447                    if ( pdebug>=2 )
     448                    { "// 18 added root to list nares"; }
     449                }
     450            }
     451
     452            nares_size=size(nares);
     453            if ( nares_size == 0 )
     454            {
     455                "Numerical problem: No root found...";
     456                "Output may be incorrect!";
     457                nares=@ln;
     458            }
     459             
     460            if ( pdebug>=1 )
     461            { "// 20 found <"+string(size(nares))+"> roots"; }
     462
     463            for ( i= 1; i <= nares_size; i++ )
     464            {
     465                res[nv-d]= nares[i];
     466
     467                if ( dd < elem )
     468                {
     469                    if ( i > 1 )
     470                    {
     471                        rn@++;
     472                    }
     473                    psubst( d+1, dd+1, n-1, res, fi, elem, nv );
     474                }
     475                else
     476                {
     477                    if ( pdebug>=1 )
     478                    {"// 30_1 <"+string(rn@)+"> "+string(size(res))+" <-----";}
     479                    if ( pdebug>=2 )
     480                    { res; }
     481                    rlist[rn@]=res;
     482                }
     483            }
     484        }
     485        else
     486        {
     487            if ( pdebug>=2 )
     488            { "// 21 found root to be: "+string(lsr[1]); }
     489            res[nv-d]= lsr[1];
     490
     491            if ( dd < elem )
     492            {
     493                psubst( d+1, dd+1, n-1, res, fi, elem, nv );
     494            }
     495            else
     496            {
     497                if ( pdebug>=1 )
     498                { "// 30_2 <"+string(rn@)+"> "+string(size(res))+" <-----";}
     499                if ( pdebug>=2 )
     500                { res; }
     501                rlist[rn@]=res;
     502            }
     503        }
     504    }
     505}
     506
     507///////////////////////////////////////////////////////////////////////////////
     508
     509proc fglm_solve( ideal fi, list # )
     510"USAGE:   fglm_solve( i [, d ] ); i ideal, p integer.
     511         p gives precision of complex numbers in digits,
     512         uses a standard basis computed from fi to determine
     513         recursively all complex roots of fi
     514RETURN:  a list of complex roots of type number.
     515NOTE:    changes the given ring to the ring ring with complex coefficient
     516         field with precision p in digits.
     517EXAMPLE: example fglm_solve; shows an example
     518"
     519{
     520    int prec=30;
     521
     522    if ( size(#)>=1  && typeof(#[1])=="int")
     523    {
     524        prec=#[1];
     525    }
     526
     527    lex_solve(stdfglm(fi),prec);
     528    keepring basering;
     529}
     530example
     531{
     532    "EXAMPLE:";echo=2;
     533    ring r = 0,(x,y),lp;
     534    // compute the intersection points of two curves
     535    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     536    fglm_solve(s);
     537    rlist;
     538}
     539
     540///////////////////////////////////////////////////////////////////////////////
     541
     542proc lex_solve( ideal fi, int prec, list # )
     543"USAGE:   lex_solve( i, d [, l] ); i ideal, p,d,l integers.
     544         i a reduced lexicographical Groebner bases of a zero-dimensional
     545         ideal (i), sorted by increasing leading terms.
     546         p gives precision of complex numbers in digits,
     547         d gives precision (1<d<p) for near-zero-determination (default:1/2*p),
     548         l gives debug level (-1: no output,..., 3: lots of output, 0: default)
     549RETURN:  a list of complex roots of type number.
     550NOTE:    changes the given ring to the ring ring with complex coefficient
     551         field with precision p in digits.
     552EXAMPLE: example lex_solve; shows an example
     553"
     554{
     555    if ( !defined(pdebug) )
     556    {
     557        int pdebug;
     558        export pdebug;
     559    }
     560    pdebug=0;
     561
     562    string orings= nameof(basering);
     563    def oring= basering;
     564
     565    // change the ground field to complex numbers 
     566    string nrings= "ring "+orings+"C=(real,"+string(prec)
     567        +",I),("+varstr(basering)+"),lp;";
     568    execute nrings;
     569
     570    if ( pdebug>=0 )
     571    { "// name of new ring: "+string(nameof(basering));}
     572
     573    // map fi from old to new ring
     574    ideal fi= imap(oring,fi);
     575
     576    // list with entry 0 (number)
     577    number nn=0;
     578    if ( !defined(@ln) )
     579    {
     580        list @ln;
     581        export @ln;
     582    }
     583    @ln=nn;
     584
     585    // set number of digits for zero-comparison of roots
     586    if ( !defined(myCompDigits) )
     587    {
     588        int myCompDigits;
     589        export  myCompDigits;
     590    }
     591    if ( size(#)>=1  && typeof(#[1])=="int")
     592    {
     593        myCompDigits=#[1];
     594    }
     595    else
     596    {
     597        myCompDigits=(system("getPrecDigits")/2);
     598    }
     599
     600    if ( pdebug>=1 )
     601    {"// myCompDigits="+string(myCompDigits);}
     602
     603    int idelem= size(fi);
     604    int nv= nvars(basering);
     605    int i,j,k,lis;
     606    list res,li;
     607
     608    if ( !defined(rlist) )
     609    {
     610        list rlist;
     611        export rlist;
     612    }
     613   
     614    if ( !defined(rn@) )
     615    {
     616        int rn@;
     617        export rn@;
     618    }
     619    rn@=0;
     620
     621    li= laguerre_solve(fi[1]);
     622    lis= size(li);
     623
     624    if ( pdebug>=1 )
     625    {"// laguerre found roots: "+string(size(li));}
     626   
     627    for ( j= 1; j <= lis; j++ )
     628    {
     629        if ( pdebug>=1 )
     630        {"// root "+string(j);}
     631        rn@++;
     632        res[nv]= li[j];
     633        psubst( 1, 2, nv-1, res, fi, idelem, nv );
     634    }
     635   
     636    if ( pdebug>=0 )
     637    {"// list of roots: "+nameof(rlist);}
     638
     639    // keep the ring and exit
     640    keepring basering;
     641}
     642example
     643{
     644    "EXAMPLE:";echo=2;
     645    ring r = 0,(x,y),lp;
     646    // compute the intersection points of two curves
     647    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     648    lex_solve(stdfglm(s),30);
     649    rlist;
     650}
     651
     652///////////////////////////////////////////////////////////////////////////////
     653
     654proc triangLf_solve( ideal fi, list # )
     655"USAGE:   triangLf_solve( i [, d ] ); i ideal, p integer.
     656         i zero-dimensional ideal
     657         p gives precision of complex numbers in digits,
     658         uses a triangular system (Lazard's Algorithm with factorization)
     659         computed from a standard basis to determine recursively all complex
     660         roots with Laguerre's algorithm of input ideal i
     661RETURN:  a list of complex roots of type number.
     662NOTE:    changes the given ring to the ring ring with complex coefficient
     663         field with precision p in digits.
     664EXAMPLE: example triangLf_solve; shows an example
     665"
     666{
     667    int prec=30;
     668
     669    if ( size(#)>=1  && typeof(#[1])=="int")
     670    {
     671        prec=#[1];
     672    }
     673
     674    triang_solve(triangLfak(stdfglm(fi)),prec);
     675    keepring basering;
     676}
     677example
     678{
     679    "EXAMPLE:";echo=2;
     680    ring r = 0,(x,y),lp;
     681    // compute the intersection points of two curves
     682    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     683    triangLf_solve(s);
     684    rlist;
     685}
     686
     687///////////////////////////////////////////////////////////////////////////////
     688
     689proc triangM_solve( ideal fi, list # )
     690"USAGE:   triangM_solve( i [, d ] ); i ideal, p integer.
     691         i zero-dimensional ideal
     692         p gives precision of complex numbers in digits,
     693         uses a triangular system (Moellers Algorithm) computed from a standard
     694         basis to determine recursively all complex roots with Laguerre's
     695         algorithm of input ideal i
     696RETURN:  a list of complex roots of type number.
     697NOTE:    changes the given ring to the ring ring with complex coefficient
     698         field with precision p in digits.
     699EXAMPLE: example triangM_solve; shows an example
     700"
     701{
     702    int prec=30;
     703
     704    if ( size(#)>=1  && typeof(#[1])=="int")
     705    {
     706        prec=#[1];
     707    }
     708
     709    triang_solve(triangM(stdfglm(fi)),prec);
     710    keepring basering;
     711}
     712example
     713{
     714    "EXAMPLE:";echo=2;
     715    ring r = 0,(x,y),lp;
     716    // compute the intersection points of two curves
     717    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     718    triangM_solve(s);
     719    rlist;
     720}
     721
     722///////////////////////////////////////////////////////////////////////////////
     723
     724proc triangL_solve( ideal fi, list # )
     725"USAGE:   triangL_solve( i [, d ] ); i ideal, p integer.
     726         i zero-dimensional ideal
     727         p gives precision of complex numbers in digits,
     728         uses a triangular system (Lazard's Algorithm)
     729         computed from a standard basis to determine recursively all complex
     730         roots with Laguerre's algorithm of input ideal i
     731RETURN:  a list of complex roots of type number.
     732NOTE:    changes the given ring to the ring ring with complex coefficient
     733         field with precision p in digits.
     734EXAMPLE: example triangL_solve; shows an example
     735"
     736{
     737    int prec=30;
     738
     739    if ( size(#)>=1  && typeof(#[1])=="int")
     740    {
     741        prec=#[1];
     742    }
     743
     744    triang_solve(triangL(stdfglm(fi)),prec);
     745    keepring basering;
     746}
     747example
     748{
     749    "EXAMPLE:";echo=2;
     750    ring r = 0,(x,y),lp;
     751    // compute the intersection points of two curves
     752    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     753    triangL_solve(s);
     754    rlist;
     755}
     756
     757
     758///////////////////////////////////////////////////////////////////////////////
     759
     760proc triang_solve( list lfi, int prec, list # )
     761"USAGE:   triang_solve( i, d [, l] ); l list, p,d,l integers.
     762         l a list of finitely many triangular systems, such that
     763         the union of their varieties equals the variety of the
     764         initial ideal.
     765         p gives precision of complex numbers in digits,
     766         d gives precision (1<d<p) for near-zero-determination (default:1/2*p),
     767         l gives debug level (-1: no output,..., 3: lots of output, 0: default)
     768ASSUME:  l was computed using Algorithm of Lazard or
     769         Algorithm of Moeller (see triang.lib).
     770RETURN:  a list of complex roots of type number.
     771NOTE:    changes the given ring to the ring ring with complex coefficient
     772         field with precision p in digits.
     773EXAMPLE: example triang_solve; shows an example
     774"
     775{
     776    if ( !defined(pdebug) )
     777    {
     778        int pdebug;
     779        export pdebug;
     780    }
     781    pdebug=0;
     782
     783    string orings= nameof(basering);
     784    def oring= basering;
     785
     786    // change the ground field to complex numbers 
     787    string nrings= "ring "+orings+"C=(real,"+string(prec)
     788        +",I),("+varstr(basering)+"),lp;";
     789    execute nrings;
     790
     791    if ( pdebug>=0 )
     792    { "// name of new ring: "+string(nameof(basering));}
     793
     794    // list with entry 0 (number)
     795    number nn=0;
     796    if ( !defined(@ln) )
     797    {
     798        list @ln;
     799        export @ln;
     800    }
     801    @ln=nn;
     802
     803    // set number of digits for zero-comparison of roots
     804    if ( !defined(myCompDigits) )
     805    {
     806        int myCompDigits;
     807        export  myCompDigits;
     808    }
     809    if ( size(#)>=1  && typeof(#[1])=="int" )
     810    {
     811        myCompDigits=#[1];
     812    }
     813    else
     814    {
     815        myCompDigits=(system("getPrecDigits")/2);
     816    }
     817
     818    if ( pdebug>=1 )
     819    {"// myCompDigits="+string(myCompDigits);}
     820
     821    int idelem;
     822    int nv= nvars(basering);
     823    int i,j,k,lis;
     824    list res,li;
     825
     826    if ( !defined(rlist) )
     827    {
     828        list rlist;
     829        export rlist;
     830    }
     831   
     832    if ( !defined(rn@) )
     833    {
     834        int rn@;
     835        export rn@;
     836    }
     837    rn@=0;
     838
     839    // map the list
     840    list lfi= imap(oring,lfi);
     841
     842    int slfi= size(lfi);
     843    ideal fi;
     844
     845    for ( i= 1; i <= slfi; i++ )
     846    {
     847        // map fi from old to new ring
     848        fi= lfi[i]; //imap(oring,lfi[i]);
     849       
     850        idelem= size(fi);
     851
     852        // solve fi[1]
     853        li= laguerre_solve(fi[1]);
     854        lis= size(li);
     855
     856        if ( pdebug>=1 )
     857        {"// laguerre found roots: "+string(size(li));}
     858       
     859        for ( j= 1; j <= lis; j++ )
     860        {
     861            if ( pdebug>=1 )
     862            {"// root "+string(j);}
     863            rn@++;
     864            res[nv]= li[j];
     865            psubst( 1, 2, nv-1, res, fi, idelem, nv );
     866        }
     867    }
     868
     869    if ( pdebug>=0 )
     870    {"// list of roots: "+nameof(rlist);}
     871       
     872    // keep the ring and exit
     873    keepring basering;
     874}
     875example
     876{
     877    "EXAMPLE:";echo=2;
     878    ring r = 0,(x,y),lp;
     879    // compute the intersection points of two curves
     880    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     881    triang_solve(triangLfak(stdfglm(s)),30);
     882    rlist;
     883}
     884
     885
     886///////////////////////////////////////////////////////////////////////////////
     887
     888proc pcheck( ideal fi, list sols, list # )
     889"USAGE:   pcheck( i, l [, n] ); i ideal, l list, n integer
     890          i system of equations, l list of numers assumend
     891          to be roots of i.
     892RETURN:  1 iff all elements of l are roots of i, else 0
     893EXAMPLE: example pcheck; shows an example
     894"
     895{
     896    int i,j,k,nnrfound;
     897    int ss= size(sols);
     898    int nv= nvars(basering);
     899    poly ps;
     900    number nn;
     901    int cprec=0;
     902   
     903    if ( size(#)>=1  && typeof(#[1])=="int" )
     904    {
     905        cprec=#[1];
     906    }
     907    if ( cprec <= 0 )
     908    {
     909        cprec=system("getPrecDigits")/2;
     910    }
     911
     912    nnrfound=1;
     913    for ( j= 1; j <= size(fi); j++ )
     914    {
     915        for ( i= 1; i <= ss; i++ )
     916        {
     917            ps= fi[j];
     918            for ( k= 1; k <= nv; k++ )
     919            {
     920                ps = subst(ps,var(k),sols[i][k]);
     921            }
     922            //ps;
     923            nn= leadcoef(ps);
     924            if ( nn != 0 )
     925            {
     926                if ( !system("complexNearZero",nn,cprec) )
     927                {
     928                    " no root: ideal["+string(j)+"]( root "
     929                        +string(i)+"): 0 != "+string(ps);
     930                    nnrfound=0;
     931                }
     932            }
     933        }
     934    }
     935    return(nnrfound);
     936}
     937example
     938{
     939    "EXAMPLE:";echo=2;
     940    ring r = 0,(x,y),lp;
     941    // compute the intersection points of two curves
     942    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
     943    lex_solve(stdfglm(s),30);
     944    rlist;
     945    pcheck(rlist);
     946}
     947
     948
     949// local Variables: ***
     950// c-set-style: bsd ***
     951// End: ***

  • Singular/mpr_complex.cc

    r0c85269 r5c67581  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: mpr_complex.cc,v 1.17 1999-10-14 14:27:22 obachman Exp $ */
     4/* $Id: mpr_complex.cc,v 1.18 1999-11-14 21:35:09 wenk Exp $ */
    55
    66/*
    … …  
    627627}
    628628//<-
     629
     630bool complexNearZero( gmp_complex * c, int digits )
     631{
     632  gmp_float eps,epsm;
     633
     634  if ( digits < 1 ) return true;
     635
     636  eps=pow(10.0,(int)digits);
     637  //Print("eps: %s\n",floatToStr(eps,gmp_output_digits));
     638  eps=(gmp_float)1.0/eps;
     639  epsm=-eps;
     640
     641  //Print("eps: %s\n",floatToStr(eps,gmp_output_digits));
     642
     643  if ( c->real().sign() > 0 ) // +
     644    return (c->real() < eps && (c->imag() < eps && c->imag() > epsm));
     645  else // -
     646    return (c->real() > epsm && (c->imag() < eps && c->imag() > epsm));
     647}
     648
    629649//%e
    630650

  • Singular/mpr_complex.h

    r0c85269 r5c67581  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: mpr_complex.h,v 1.9 1999-10-14 14:27:23 obachman Exp $ */
     6/* $Id: mpr_complex.h,v 1.10 1999-11-14 21:35:09 wenk Exp $ */
    77
    88/*
    … …  
    392392//<-
    393393
     394bool complexNearZero( gmp_complex * c, int digits );
     395
    394396#endif MPR_COMPLEX_H
    395397
Note: See TracChangeset for help on using the changeset viewer.