Changeset 960633 in git


Ignore:
Timestamp:
Nov 1, 2017, 3:04:36 AM (6 years ago)
Author:
Janko Boehm <boehm@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
6824e1cd64288ace28095cc51dfc7db2b75deeaf
Parents:
0b78529f4ca9d003d3c81840a34aa6b478ebb9c3
Message:
Update rootisolation.lib using numerical solver in univariate case
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/rootisolation.lib

    r0b7852 r960633  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="version rootisolation.lib 4.1.0.0 Jun_2017 ";
     2version="version rootisolation.lib 4.1.0.3 Aug_2017 "; // $Id$
    33info="
    44LIBRARY:    rootisolation.lib implements an algorithm for real root isolation
    55            using interval arithmetic
    66AUTHORS:    Dominik Bendle (bendle@rhrk.uni-kl.de)
     7            Janko Boehm (boehm@mathematik.uni-lk.de), supervisor Fachpraktikum
    78@*          Clara Petroll (petroll@rhrk.uni-kl.de)
    89
    9 OVERVIEW:   In this library the interval arithmetic from \"@code{interval.so}\"
    10             is used. The new type \"@code{ivmat}\", a matrix consiting of
    11             intervals is implemented as \"@code{newstruct}\". There are various
     10OVERVIEW:   In this library the interval arithmetic from @code{interval.so}
     11            is used. The new type @code{ivmat}, a matrix consisting of
     12            intervals, is implemented as @code{newstruct}. There are various
    1213            functions for computations with interval matrices implemented, such
    1314            as Gaussian elimination for interval matrices.
    1415
    1516@*          Interval arithmetic, the interval Newton Step and exclusion methods
    16             are used to implement the procedure \"@code{rootIsolation}\", an
     17            are used to implement the procedure @code{rootIsolation}, an
    1718            algorithm which finds boxes containing elements of the vanishing
    1819            locus of an ideal. This algorithm is specialised for
    … …  
    4142PROCEDURES:
    4243bounds(a,#);            creates a new interval with given bounds
     44length(I);              returns Euclidean length of interval
     45boxSet(B,i,I);          returns box B with B[i]==I
    4346ivmatInit(m, n);        returns m x n [0,0]-matrix
    4447ivmatSet(A,i,j,I);      returns matrix A where A[i][j]=I
    4548unitMatrix(m);          returns m x m unit matrix where 1 = [1,1]
    4649ivmatGaussian(M);       computes M^(-1) using Gaussian elimination for intervals
     50evalPolyAtBox(f,B);     returns evaluation of polynomial at a box
    4751evalJacobianAtBox(A,B); jacobian matrix of A where polynomials are evaluated at B
    4852
    … …  
    5054                        computes boxes containing unique roots of I lying in L
    5155rootIsolation(I,B,e);   slims down input box B and calls rootIsolationNoPreprocessing
     56rootIsolationPrimdec(I,B,e);   runs a primary decomposition primdecGTZ before root isoation
     57
     58
     59KEYWORDS:   real root isolation; interval arithmetic; Newton step
    5260";
    53 LIB "atkins.lib"; // for round (tmp?)
     61LIB "presolve.lib";
     62LIB "rootsur.lib";
     63LIB "primdec.lib";
     64LIB "solve.lib";
     65LIB "atkins.lib";
     66
    5467///////////////////////////////////////////////////////////////////////////////
    55 
    5668static proc mod_init()
    5769{
    58     LIB "interval.so";
    59 
    60     if (!reservedName("ivmat"))
    61     {
    62       newstruct("ivmat", "list rows");
    63       system("install", "ivmat", "print", ivmatPrint,           1);
    64       system("install", "ivmat", "[",     ivmatGet,             2);
    65       system("install", "ivmat", "nrows", ivmatNrows,           1);
    66       system("install", "ivmat", "ncols", ivmatNcols,           1);
    67       system("install", "ivmat", "*",     ivmatMultiplyGeneral, 2);
    68     }
     70  intvec opt = option(get);
     71  option(noredefine);
     72  LIB "interval.so";
     73  option(set,opt);
     74
     75  if (!reservedName("ivmat"))
     76  {
     77    newstruct("ivmat", "list rows");
     78    system("install", "ivmat", "print", ivmatPrint,           1);
     79    system("install", "ivmat", "[",     ivmatGet,             2);
     80    system("install", "ivmat", "nrows", ivmatNrows,           1);
     81    system("install", "ivmat", "ncols", ivmatNcols,           1);
     82    system("install", "ivmat", "*",     ivmatMultiplyGeneral, 2);
     83  }
    6984}
    7085
    7186///////////////////////////////////////////////////////////////////////////////
     87
     88// HELPER FUNCTIONS
     89
     90static proc floor(number a)
     91{
     92  // a = m/n
     93  bigint m, n = bigint(numerator(a)), bigint(denominator(a));
     94  // div for bigint seemingly works like floor
     95  return(number(m div n));
     96}
     97
     98static proc ceil(number a)
     99{
     100  return(-floor(-a));
     101}
    72102
    73103// INTERVAL FUNCTIONS
    … …  
    81111EXAMPLE: example bounds; create two intervals"
    82112{
    83     interval I;
    84     if (size(#) == 0)
    85     {
    86         I = a;
    87         return(I);
    88     }
    89     if (size(#) == 1 && (typeof(#[1]) == "number" || typeof(#[1]) == "int"))
    90     {
    91         I = a, #[1];
    92         return(I);
    93     }
    94     ERROR("syntax: bounds(<number>) or bounds(<number>, <number>)");
    95 }
    96 example
    97 {
    98     "EXAMPLE:"; echo = 2;
    99     ring R = 0,x,dp;
    100     interval I = bounds(1);
    101     I;
    102 
    103     interval J = bounds(2/3,3);
    104     J;
     113  interval I;
     114  if (size(#) == 0)
     115  {
     116    I = a;
     117    return(I);
     118  }
     119  if (size(#) == 1 && (typeof(#[1]) == "number" || typeof(#[1]) == "int"))
     120  {
     121    I = a, #[1];
     122    return(I);
     123  }
     124  ERROR("syntax: bounds(<number>) or bounds(<number>, <number>)");
     125}
     126example
     127{
     128  "EXAMPLE:"; echo = 2;
     129  ring R = 0,x,dp;
     130  interval I = bounds(1);
     131  I;
     132
     133  interval J = bounds(2/3,3);
     134  J;
     135}
     136
     137// function defined in interval.so
     138proc length(interval)
     139"USAGE:  @code{length(I)}; @code{I interval}
     140RETURN: number, the Euclidean length of the interval
     141EXAMPLE: example length; shows an example"
     142{
     143}
     144example
     145{
     146  "EXAMPLE:"; echo = 2;
     147  ring R = 0,x,dp;
     148  interval I = -1,3;
     149  length(I);
     150  I = 1/5,1/3;
     151  length(I);
    105152}
    106153
    107154// BOX FUNCTIONS
    108155
    109 static proc lengthBox(box B)
    110 "USAGE:  length(B), B box
     156proc boxSet(box, int, interval)
     157"USAGE: @code{boxSet(B, i, I)}; @code{B box, i int, I interval}
     158RETURN: new box @code{C} where @code{C[i]==I}
     159PURPOSE: modifies a single entry of a box
     160EXAMPLE: example boxSet; shows an example"
     161{
     162}
     163example
     164{
     165  "EXAMPLE:"; echo = 2;
     166  ring R = 0,(x,y,z),dp;
     167  box B; B;
     168  B = boxSet(B, 2, bounds(-1,1)); B;
     169  B = boxSet(B, 1, B[2]); B;
     170}
     171
     172static proc boxLength(box B)
     173"USAGE:  boxLength(B), B box
    111174RETURN: length/size in measure sense"
    112175{
    113     number maximum = 0;
    114     int n = nvars(basering);
    115 
    116     for (int i=1; i <= n; i++)
    117     {
    118         maximum = max(maximum, length(B[i]));
    119     }
    120 
    121     return(maximum);
    122 }
    123 
    124 static proc boxCenter(box M)
    125 "USAGE: boxCenter(M); M ivmat
    126 RETURN: box containing center elements of M"
    127 {
    128     int n = nvars(basering);
    129 
    130     list C;
    131     int i;
    132 
    133     for (i = 1; i <= n; i++)
    134     {
    135         C[i] = interval((M[i][1] + M[i][2])/2);
    136     }
    137 
    138     return(box(C));
     176  number maxLength, l = 0, 0;
     177  int n = nvars(basering);
     178
     179  for (int i=1; i <= n; i++)
     180  {
     181    l = length(B[i]);
     182    if (maxLength < l)
     183    {
     184      maxLength = l;
     185    }
     186  }
     187
     188  return(maxLength);
     189}
     190
     191static proc boxEmptyIntervals(box B)
     192{
     193  int N = nvars(basering);
     194  intvec z;
     195  for (int i = 1; i <= N; i++)
     196  {
     197    z[i] = length(B[i]) == 0;
     198  }
     199  return(z);
     200}
     201
     202static proc boxCenter(box B)
     203"USAGE: boxCenter(B); B box
     204RETURN: box containing center elements of B"
     205{
     206  int n = nvars(basering);
     207
     208  list C;
     209  int i;
     210
     211  for (i = 1; i <= n; i++)
     212  {
     213    C[i] = interval((B[i][1] + B[i][2])/2);
     214  }
     215
     216  return(box(C));
    139217}
    140218
    … …  
    144222        contain zeros of I
    145223NOTE:   this uses exclusion tests and Groebner bases to determine whether the
    146         intersection plane contains a root of I
    147 EXAMPLE: example splitBox; splits two-dimensional interval into two"
    148 {
    149     // at first split only at largest interval
    150     int imax = 1;
    151     int N = nvars(basering);
    152 
    153     for (int i = 2; i <= N; i++)
    154     {
    155         if (length(B[i]) > length(B[imax])) { imax = i; }
    156     }
    157 
    158     number ratio = 1/2;
    159     number mean;
    160     box intersection;
    161     ideal Inew;
    162 
    163     while(1)
    164     {
    165         mean = ratio * B[imax][1] + (1 - ratio) * B[imax][2];
    166 
    167         // note this works only for ideals with N generators or less
    168         intersection = evalIdealAtBox(I, boxSet(B, imax, interval(mean)));
    169         for (i = 1; i <= N; i++)
    170         {
    171             // check if any interval does not contain zero
    172             if (intersection[i][1]*intersection[i][2] > 0) { break; }
    173         }
    174 
    175         Inew = I + (var(imax) - mean);
    176         // check if groebner basis is trivial
    177         if (std(Inew) == 1) { break; }
    178 
    179         // else there must?/might be a zero on the intersection,
    180         // so decrease ratio slightly
    181         ratio = ratio * 15/16;
    182 
    183         // make sure algorithm terminates after taking too many steps
    184         // this may not be necessary
    185         if ( ratio < 1/100 )
    186         {
    187             print("splitBox took too long");
    188             break;
    189         }
    190     }
    191 
    192     // now split boxes
    193     box boxLeft  = boxSet(B, imax, bounds(B[imax][1], mean));
    194     box boxRight = boxSet(B, imax, bounds(mean, B[imax][2]));
    195 
    196     return(boxLeft, boxRight);
    197 }
    198 example
    199 {
    200     "EXAMPLE:"; echo = 2;
    201     ring R = 0,(x,y),dp;
    202 
    203     box B = list(bounds(0,1),
    204                  bounds(0,2));
    205 
    206     B;
    207     splitBox(B,1);
    208     splitBox(B,y-1); // contains zero on first splitting plane candidate
     224        intersection plane contains a root of I"
     225{
     226  // at first split only at largest interval
     227  int imax = 1;
     228  int N = nvars(basering);
     229
     230  for (int i = 2; i <= N; i++)
     231  {
     232    if (length(B[i]) > length(B[imax]))
     233    {
     234      imax = i;
     235    }
     236  }
     237
     238  number ratio = 1/2;
     239  number mean;
     240  box intersection;
     241  ideal Inew;
     242
     243  while(1)
     244  {
     245    mean = ratio * B[imax][1] + (1 - ratio) * B[imax][2];
     246
     247    // note this works only for ideals with N generators or less
     248    intersection = evalIdealAtBox(I, boxSet(B, imax, interval(mean)));
     249    for (i = 1; i <= N; i++)
     250    {
     251      // check if any interval does not contain zero
     252      if (intersection[i][1]*intersection[i][2] > 0)
     253      {
     254        break;
     255      }
     256    }
     257
     258    Inew = I + (var(imax) - mean);
     259    // check if groebner basis is trivial
     260    if (std(Inew) == 1)
     261    {
     262      break;
     263    }
     264
     265    // else there must?/might be a zero on the intersection, so decrease ratio
     266    // slightly
     267    ratio = ratio * 15/16;
     268
     269    // make sure algorithm terminates after taking too many steps this may not
     270    // be necessary
     271    if ( ratio < 1/100 )
     272    {
     273      print("splitBox took too long");
     274      break;
     275    }
     276  }
     277
     278  // now split boxes
     279  box boxLeft  = boxSet(B, imax, bounds(B[imax][1], mean));
     280  box boxRight = boxSet(B, imax, bounds(mean, B[imax][2]));
     281
     282  return(boxLeft, boxRight);
    209283}
    210284
    … …  
    214288EXAMPLE: example boxIsInterior; check whether A is contained in int(B)"
    215289{
    216     int N = nvars(basering);
    217     for (int i = 1; i <= N; i++)
    218     {
    219         if (A[i][1] <= B[i][1] || A[i][2] >= B[i][2]) { return(0); }
    220     }
    221     return(1);
    222 }
    223 example
    224 {
    225     "EXAMPLE:"; echo=2;
    226 
    227     ring R = 0,(x,y,z), lp;
    228 
    229     box A = list(bounds(1,2), bounds(2,3), bounds(1/2,7/2)); A;
    230     box B1 = list(bounds(0,5/2), bounds(1,4), bounds(0,9)); B1;
    231     boxIsInterior(A,B1);
    232 
    233     box B2 = list(bounds(2,4), bounds(1,4), bounds(0,9)); B2;
    234     boxIsInterior(A,B2);
     290  int N = nvars(basering);
     291  for (int i = 1; i <= N; i++)
     292  {
     293    if (A[i][1] <= B[i][1] || A[i][2] >= B[i][2])
     294    {
     295      return(0);
     296    }
     297  }
     298  return(1);
     299}
     300example
     301{
     302  "EXAMPLE:"; echo=2;
     303
     304  ring R = 0,(x,y,z), lp;
     305
     306  box A = list(bounds(1,2), bounds(2,3), bounds(1/2,7/2)); A;
     307  box B1 = list(bounds(0,5/2), bounds(1,4), bounds(0,9)); B1;
     308  boxIsInterior(A,B1);
     309
     310  box B2 = list(bounds(2,4), bounds(1,4), bounds(0,9)); B2;
     311  boxIsInterior(A,B2);
    235312}
    236313
    … …  
    239316// MATRIX FUNCTIONS
    240317
    241 proc ivmatInit(numrows, numcols)
     318proc ivmatInit(int numrows, int numcols)
    242319"USAGE: @code{ivmatInit(m, n)}; @code{m, n int}
    243320RETURN: @code{m}x@code{n} matrix of [0,0]-intervals
    … …  
    246323EXAMPLE: example ivmatInit; initialises an interval matrix"
    247324{
    248     ivmat A;
    249     A.rows = list();
    250     int i, j;
    251     interval z = 0;
    252 
    253     for (i = 1; i <= numrows; i++)
    254     {
    255         A.rows[i] = list();
    256         for (j=1; j <= numcols; j++)
    257         {
    258             A.rows[i][j] = z;
    259         }
    260     }
    261 
    262     return(A);
    263 }
    264 example
    265 {
    266     "EXAMPLE:"; echo = 2;
    267     ring R = 0,x(1..5),dp;
    268 
    269     ivmat A = ivmatInit(4, 5); A;
     325  ivmat A;
     326  A.rows = list();
     327  int i, j;
     328  interval z = 0;
     329
     330  for (i = 1; i <= numrows; i++)
     331  {
     332    A.rows[i] = list();
     333    for (j=1; j <= numcols; j++)
     334    {
     335      A.rows[i][j] = z;
     336    }
     337  }
     338
     339  return(A);
     340}
     341example
     342{
     343  "EXAMPLE:"; echo = 2;
     344  ring R = 0,x(1..5),dp;
     345
     346  ivmat A = ivmatInit(4, 5); A;
    270347}
    271348
    … …  
    275352EXAMPLE: example ivmatNrows; return number of rows"
    276353{
    277     return(size(M.rows));
    278 }
    279 example
    280 {
    281     "EXAMPLE:"; echo = 2;
    282     ring R = 0,x,dp;
    283 
    284     ivmat A = ivmatInit(2,3);
    285     nrows(A);
     354  return(size(M.rows));
     355}
     356example
     357{
     358  "EXAMPLE:"; echo = 2;
     359  ring R = 0,x,dp;
     360
     361  ivmat A = ivmatInit(2,3);
     362  nrows(A);
    286363}
    287364
    288365static proc ivmatNcols(ivmat M)
    289366{
    290     return(size(M.rows[1]));
     367  return(size(M.rows[1]));
    291368}
    292369
    293370static proc ivmatPrint(ivmat A)
    294371{
    295     int m = nrows(A);
    296     for (int i = 1; i <= m; i++)
    297     {
    298         string(A.rows[i]);
    299     }
     372  int m = nrows(A);
     373  for (int i = 1; i <= m; i++)
     374  {
     375    string(A.rows[i]);
     376  }
    300377}
    301378
    … …  
    304381RETURN: list A[i] of i-th row of A"
    305382{
    306     return(A.rows[i]);
     383  return(A.rows[i]);
    307384}
    308385
    … …  
    313390EXAMPLE: example ivmatSet; assign values to A"
    314391{
    315     A.rows[i][j] = I;
    316     return(A);
    317 }
    318 example
    319 {
    320     "EXAMPLE:"; echo = 2;
    321     ring R = 0,x,dp;
    322     ivmat A = ivmatInit(2,2);             A;
    323     A = ivmatSet(A, 1, 2, bounds(1, 2));  A;
     392  A.rows[i][j] = I;
     393  return(A);
     394}
     395example
     396{
     397  "EXAMPLE:"; echo = 2;
     398  ring R = 0,x,dp;
     399  ivmat A = ivmatInit(2,2);             A;
     400  A = ivmatSet(A, 1, 2, bounds(1, 2));  A;
    324401}
    325402
    … …  
    329406EXAMPLE: example diagMatrix; create diagonal matrix"
    330407{
    331     ivmat E = ivmatInit(n, n);
    332     for (int i = 1; i <= n; i++)
    333     {
    334         E.rows[i][i] = I;
    335     }
    336     return(E);
    337 }
    338 example
    339 {
    340     "EXAMPLE:"; echo = 2;
    341     ring R = 0,x,dp;
    342     ivmat A = diagMatrix(2, bounds(1, 2)); A;
     408  ivmat E = ivmatInit(n, n);
     409  for (int i = 1; i <= n; i++)
     410  {
     411    E.rows[i][i] = I;
     412  }
     413  return(E);
     414}
     415example
     416{
     417  "EXAMPLE:"; echo = 2;
     418  ring R = 0,x,dp;
     419  ivmat A = diagMatrix(2, bounds(1, 2)); A;
    343420}
    344421
    … …  
    348425EXAMPLE: example unitMatrix; generate a unit matrix"
    349426{
    350     return(diagMatrix(n, 1));
    351 }
    352 example
    353 {
    354     "EXAMPLE:"; echo = 2;
    355     ring R = 0,(x,y),dp;
    356     unitMatrix(2);
    357     unitMatrix(3);
     427  return(diagMatrix(n, 1));
     428}
     429example
     430{
     431  "EXAMPLE:"; echo = 2;
     432  ring R = 0,(x,y),dp;
     433  unitMatrix(2);
     434  unitMatrix(3);
    358435}
    359436
    360437static proc ivmatMultiply(ivmat A, ivmat B)
    361438{
    362     int m = nrows(A);
    363     int n = ncols(B);
    364     int p = ncols(A);
    365 
    366     if (p <> nrows(B))
    367     {
    368         ERROR("Matrices have wrong dimensions!");
    369     }
    370 
    371     ivmat C = ivmatInit(m, n);
    372     int i, j, k;
    373     interval I;
    374 
    375     for (i = 1; i <= m; i++)
    376     {
    377         for (j = 1; j <= n; j++)
    378         {
    379             I = 0;
    380             for (k = 1; k <= p; k++)
    381             {
    382                 I = I + A[i][k] * B[k][j];
    383             }
    384             C.rows[i][j] = I;
    385         }
    386     }
    387 
    388     return(C);
     439  int m = nrows(A);
     440  int n = ncols(B);
     441  int p = ncols(A);
     442
     443  if (p <> nrows(B))
     444  {
     445    ERROR("Matrices have wrong dimensions!");
     446  }
     447
     448  ivmat C = ivmatInit(m, n);
     449  int i, j, k;
     450  interval I;
     451
     452  for (i = 1; i <= m; i++)
     453  {
     454    for (j = 1; j <= n; j++)
     455    {
     456      I = 0;
     457      for (k = 1; k <= p; k++)
     458      {
     459        I = I + A[i][k] * B[k][j];
     460      }
     461      C.rows[i][j] = I;
     462    }
     463  }
     464
     465  return(C);
    389466}
    390467
    391468proc ivmatGaussian(ivmat A)
    392469"USAGE:  @code{ivmatGaussian(A)}; @code{A ivmat}
    393 RETURN: 0 if @code{A} not invertible, @code{(1, Ainv)} if @code{A} invertible
     470RETURN: 0, if @code{A} not invertible, @code{(1, Ainv)} if @code{A} invertible
    394471        where @code{Ainv} is the inverse matrix
    395472PURPOSE: Inverts an interval matrix using Gaussian elimination in the setting
    … …  
    398475EXAMPLE: example ivmatGaussian; inverts a matrix"
    399476{
    400     int n = nrows(A);
    401     if (n <> ncols(A))
    402     {
    403         ERROR("Matrix non-square");
    404     }
    405 
    406     ivmat Ainv = unitMatrix(n);
    407     list tmp;
    408     interval TMP;
    409 
    410     int i, j, pos;
    411     for (pos = 1; pos <= n; pos++)
    412     {
    413         i = pos;
    414 
    415         // get non-zero interval on diagonal
    416         while(A[i][pos][1] * A[i][pos][2] <= 0)
    417         {
    418             i++;
    419             // if no non-zero intervals exist, then matrix must be singular
    420             if (i > n)
    421             {
    422                 return(0);
    423             }
    424         }
    425         if (i <> pos)
    426         {
    427             tmp = A.rows[i];
    428             A.rows[i] = A.rows[pos];
    429             A.rows[pos] = tmp;
    430 
    431             tmp = Ainv.rows[i];
    432             Ainv.rows[i] = Ainv.rows[pos];
    433             Ainv.rows[pos] = tmp;
    434         }
    435 
    436         // pivot (pos,pos)
    437         TMP = A[pos][pos];
    438         A.rows[pos][pos] = interval(1);
    439 
     477  int n = nrows(A);
     478  if (n <> ncols(A))
     479  {
     480    ERROR("Matrix non-square");
     481  }
     482
     483  ivmat Ainv = unitMatrix(n);
     484  list tmp;
     485  interval TMP;
     486
     487  int i, j, pos;
     488  for (pos = 1; pos <= n; pos++)
     489  {
     490    i = pos;
     491
     492    // get non-zero interval on diagonal
     493    while(A[i][pos][1] * A[i][pos][2] <= 0)
     494    {
     495      i++;
     496      // if no non-zero intervals exist, then matrix must be singular
     497      if (i > n)
     498      {
     499        return(0);
     500      }
     501    }
     502    if (i <> pos)
     503    {
     504      tmp = A.rows[i];
     505      A.rows[i] = A.rows[pos];
     506      A.rows[pos] = tmp;
     507
     508      tmp = Ainv.rows[i];
     509      Ainv.rows[i] = Ainv.rows[pos];
     510      Ainv.rows[pos] = tmp;
     511    }
     512
     513    // pivot (pos,pos)
     514    TMP = A[pos][pos];
     515    A.rows[pos][pos] = interval(1);
     516
     517    for (j = 1; j <= n; j++)
     518    {
     519      if (pos <> j)
     520      {
     521        A.rows[pos][j] = A[pos][j]/TMP;
     522      }
     523      Ainv.rows[pos][j] = Ainv[pos][j]/TMP;
     524    }
     525
     526    // clear entries above and below
     527    for (i = 1; i <= n; i++)
     528    {
     529      if (i <> pos)
     530      {
     531        TMP = A[i][pos];
     532        A.rows[i][pos] = interval(0);
    440533        for (j = 1; j <= n; j++)
    441534        {
    442             if (pos <> j) { A.rows[pos][j] = A[pos][j]/TMP; }
    443             Ainv.rows[pos][j] = Ainv[pos][j]/TMP;
     535          if (j <> pos)
     536          {
     537            A.rows[i][j] = A[i][j] - A[pos][j]*TMP;
     538          }
     539          Ainv.rows[i][j] = Ainv[i][j] - Ainv[pos][j]*TMP;
    444540        }
    445 
    446         // clear entries above and below
    447         for (i = 1; i <= n; i++)
    448         {
    449             if (i <> pos)
    450             {
    451                 TMP = A[i][pos];
    452                 A.rows[i][pos] = interval(0);
    453                 for (j = 1; j <= n; j++)
    454                 {
    455                     if (j <> pos) { A.rows[i][j] = A[i][j] - A[pos][j]*TMP; }
    456                     Ainv.rows[i][j] = Ainv[i][j] - Ainv[pos][j]*TMP;
    457                 }
    458             }
    459         }
    460     }
    461     return(1, Ainv);
    462 }
    463 example
    464 {
    465     "EXAMPLE:"; echo = 2;
    466     ring R = 0,(x,y),dp;
    467 
    468     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;
    469     box B = list(bounds(7/8, 9/8), bounds(-1/10, 1/20));
    470 
    471     ivmat J = evalJacobianAtBox (I, B); J;
    472 
    473     list result = ivmatGaussian(J);
    474     ivmat Jinv = result[2];
    475     Jinv;
    476 
    477     Jinv * J;
     541      }
     542    }
     543  }
     544  return(1, Ainv);
     545}
     546example
     547{
     548  "EXAMPLE:"; echo = 2;
     549  ring R = 0,(x,y),dp;
     550
     551  ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;
     552  box B = list(bounds(7/8, 9/8), bounds(-1/10, 1/20));
     553
     554  ivmat J = evalJacobianAtBox (I, B); J;
     555
     556  list result = ivmatGaussian(J);
     557  ivmat Jinv = result[2];
     558  Jinv;
     559
     560  Jinv * J;
    478561}
    479562
    480563static proc applyMatrix(ivmat A, box b)
    481564{
    482     int n = nvars(basering);
    483 
    484     if (ncols(A) <> n || nrows(A) <> n)
    485     {
    486         ERROR("Matrix has wrong dimensions");
    487     }
    488 
    489     int i, j;
    490     list result;
    491     interval tmp;
    492 
    493     for (i = 1; i <= n; i++)
    494     {
    495         tmp = 0;
    496         for (j = 1; j <= n; j++)
    497         {
    498             tmp = tmp + A[i][j] * b[j];
    499         }
    500         result[i] = tmp;
    501     }
    502 
    503     return(box(result));
     565  int n = nvars(basering);
     566
     567  if (ncols(A) <> n || nrows(A) <> n)
     568  {
     569    ERROR("Matrix has wrong dimensions");
     570  }
     571
     572  int i, j;
     573  list result;
     574  interval tmp;
     575
     576  for (i = 1; i <= n; i++)
     577  {
     578    tmp = 0;
     579    for (j = 1; j <= n; j++)
     580    {
     581      tmp = tmp + A[i][j] * b[j];
     582    }
     583    result[i] = tmp;
     584  }
     585
     586  return(box(result));
    504587}
    505588
    506589static proc ivmatMultiplyGeneral(ivmat A, B)
    507590{
    508     if (typeof(B) == "ivmat")
    509     {
    510         return(ivmatMultiply(A, B));
    511     }
    512     if (typeof(B) == "box")
    513     {
    514         return(applyMatrix(A, B));
    515     }
    516     ERROR("Type not supported.");
     591  if (typeof(B) == "ivmat")
     592  {
     593    return(ivmatMultiply(A, B));
     594  }
     595  if (typeof(B) == "box")
     596  {
     597    return(applyMatrix(A, B));
     598  }
     599  ERROR("Type not supported.");
    517600}
    518601
    … …  
    520603
    521604// POLYNOMIAL APPLICATIONS
     605
     606proc evalPolyAtBox()
     607"USAGE: @code{evalPolyAtBox(f, B)}; @code{f poly, B box}
     608RETURN: interval, evalutaion of @code{f} at @code{B} using interval arithmetic
     609PURPOSE: computes an interval extension of the polynomial
     610EXAMPLE: example evalPolyAtBox; shows an example"
     611{
     612}
     613example
     614{
     615  "EXAMPLE:"; echo = 2;
     616  ring R = 0,(x,y),dp;
     617  poly f = x2+y-1;
     618  box B = list(bounds(-1,1), bounds(1,3)/2);
     619  interval I = evalPolyAtBox(f, B); I;
     620}
    522621
    523622proc evalJacobianAtBox(ideal I, box B)
    … …  
    526625PURPOSE: evaluates each polynomial of the Jacobian matrix of @code{I} using
    527626        interval arithmetic
    528 EXAMPLE: example evalJacobianAtBox; evalate Jacobian at box"
    529 {
    530     matrix J = jacob(I);
    531     int m = nrows(J);
    532     int n = ncols(J);
    533     ivmat M = ivmatInit(m, n);
    534 
    535     int i, j;
    536 
    537     for (i = 1; i <= m; i++)
    538     {
    539         for (j = 1; j <=n ; j++)
    540         {
    541             M.rows[i][j] = evalPolyAtBox(J[i,j], B);
    542         }
    543     }
    544     return(M);
    545 }
    546 example
    547 {
    548     "EXAMPLE:"; echo = 2;
    549     ring R = 0,(x,y),dp;
    550     ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
    551 
    552     interval J = bounds(-1,1);
    553     evalJacobianAtBox(I, list(J,J));
     627EXAMPLE: example evalJacobianAtBox; evaluate Jacobian at box"
     628{
     629  matrix J = jacob(I);
     630  int m = nrows(J);
     631  int n = ncols(J);
     632  ivmat M = ivmatInit(m, n);
     633
     634  int i, j;
     635
     636  for (i = 1; i <= m; i++)
     637  {
     638    for (j = 1; j <=n ; j++)
     639    {
     640      M.rows[i][j] = evalPolyAtBox(J[i,j], B);
     641    }
     642  }
     643  return(M);
     644}
     645example
     646{
     647  "EXAMPLE:"; echo = 2;
     648  ring R = 0,(x,y),dp;
     649  ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
     650
     651  interval J = bounds(-1,1);
     652  evalJacobianAtBox(I, list(J,J));
    554653}
    555654
    … …  
    558657RETURN: list(int, box):
    559658        -1, if ideal has no zeros in given box,
    560         1, if unique zero in given box
    561         0 if test is inconclusive;
     659        1, if unique zero in given box,
     660        0, if test is inconclusive;
    562661        box is intersection of Newton step and supplied box if applicable
    563662NOTE:   rounding is performed on fractions obtained by matrix inversion to
    … …  
    566665EXAMPLE: example testPolyBox; tests the above for intersection of ellipses."
    567666{
    568     int N = nvars(basering);
    569     int isreal = find(charstr(basering), "Float(");
    570     int i;
    571 
    572     interval tmp;
    573     number lo, up, m, n;
    574 
    575     for (i = 1; i <= ncols(I); i++)
    576     {
    577         tmp = evalPolyAtBox(I[i], B);
    578         // check if 0 contained in every interval
    579         // return -1 if not
    580         if (tmp[1]*tmp[2] > 0)
     667  int N = nvars(basering);
     668  int isreal = find(charstr(basering), "Float(");
     669  int i;
     670
     671  interval tmp;
     672  number lo, up, m, n;
     673
     674  for (i = 1; i <= ncols(I); i++)
     675  {
     676    tmp = evalPolyAtBox(I[i], B);
     677    // check if 0 contained in every interval
     678    // return -1 if not
     679    if (tmp[1]*tmp[2] > 0)
     680    {
     681      return(-1, B);
     682    }
     683  }
     684
     685  // this is always the case in our applications
     686  if (ncols(I) == N)
     687  {
     688    // calculate center as box of intervals instead of numbers so we may reuse
     689    // other procedures
     690    box Bcenter = boxCenter(B);
     691
     692    ivmat J = evalJacobianAtBox(I, B);
     693    list inverse = ivmatGaussian(J);
     694
     695    // only continue if J is invertible, i.e. J contains no singular matrix
     696    if (!inverse[1])
     697    {
     698      return(0, B);
     699    }
     700    ivmat Jinverse = inverse[2];
     701
     702    // calculate Bcenter - f(B)^(-1)f(Bcenter)
     703    box fB = evalIdealAtBox(I, Bcenter);
     704    fB = Bcenter - (Jinverse * fB);
     705
     706    // algorithm will not process box further, so do not modify
     707    int laststep = boxIsInterior(fB, B);
     708
     709    // else intersection is empty or non-trivial
     710    def bInt = intersect(B, fB);
     711
     712    // if intersection is empty bInt == -1
     713    if (typeof(bInt) == "int")
     714    {
     715      return(-1, B);
     716    }
     717
     718    // if Newton step is a single point it is a root
     719    if (boxLength(fB) == 0)
     720    {
     721      return(1,fB);
     722    }
     723
     724    if (isreal)
     725    {
     726      // fraction simplification over real basering is pointless
     727      B = bInt;
     728    }
     729    else
     730    {
     731      // attempt simplification of fractions
     732      list bb;
     733      for (i = 1; i <= N; i++)
     734      {
     735        lo = B[i][1];
     736        up = B[i][2];
     737
     738        // modify numerators of B to tighten box
     739        if (lo < bInt[i][1])
    581740        {
    582             return(-1, B);
     741          n = denominator(lo);
     742          lo = floor(bInt[i][1]*n)/n;
    583743        }
    584     }
    585 
    586     // this is always the case in our applications
    587     if (ncols(I) == N)
    588     {
    589         // calculate center as box of intervals instead of numbers
    590         // so we may reuse other procedures
    591         box Bcenter = boxCenter(B);
    592 
    593         ivmat J = evalJacobianAtBox(I, B);
    594         list inverse = ivmatGaussian(J);
    595 
    596         // only continue if J is invertible , i.e. J contains no singular matrix
    597         if (!inverse[1])
     744        if (up > bInt[i][2])
    598745        {
    599             return(0, B);
     746          n = denominator(up);
     747          up = ceil(bInt[i][2]*n)/n;
    600748        }
    601         ivmat Jinverse = inverse[2];
    602 
    603         // calculate Bcenter - f(B)^(-1)f(Bcenter)
    604         box fB = evalIdealAtBox(I, Bcenter);
    605         fB = Bcenter - (Jinverse * fB);
    606 
    607         // algorothm will not process box further, so do not modify
    608         int laststep = boxIsInterior(fB, B);
    609 
    610         // else intersection is empty or non-trivial
    611         def Bint = intersect(B, fB);
    612 
    613         // if intersection is empty Bint == -1
    614         if (typeof(Bint) == "int")
     749
     750        // make sure box does not grow
     751        if (lo >= B[i][1] && up <= B[i][2])
    615752        {
    616             return(-1, B);
    617         }
    618 
    619         if (isreal)
    620         {
    621             // fraction simplification over real basering is pointless
    622             B = Bint;
     753          bb[i] = bounds(lo, up);
    623754        }
    624755        else
    625756        {
    626             // attempt simplification of fractions
    627             // add check if we work over reals
    628             list bb;
    629             for (i = 1; i <= N; i++)
    630             {
    631                 lo = B[i][1];
    632                 up = B[i][2];
    633 
    634                 // modify numerators of B to tighten box
    635                 if (lo < Bint[i][1])
    636                 {
    637                     n = denominator(lo);
    638                     // floor
    639                     lo = round(Bint[i][1]*n - 1/2)/n;
    640                 }
    641                 if (up > Bint[i][2])
    642                 {
    643                     n = denominator(up);
    644                     // ceil
    645                     up = round(Bint[i][2]*n + 1/2)/n;
    646                 }
    647 
    648                 // make sure box does not grow
    649                 if (lo >= B[i][1] && up <= B[i][2])
    650                 {
    651                     bb[i] = bounds(lo, up);
    652                 }
    653                 else
    654                 {
    655                     bb[i] = Bint[i];
    656                 }
    657             }
    658 
    659             B = bb;
     757          bb[i] = bInt[i];
    660758        }
    661 
    662         if (laststep) { return(1, B); }
    663     }
    664 
    665     // no condition could be verified
    666     return(0, B);
    667 }
    668 example
    669 {
    670     "EXAMPLE:"; echo = 2;
    671     ring R = 0,(x,y),dp;
    672     ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
    673 
    674     interval unit = bounds(0, 1);
    675     // there may be common zeros in [0,1]x[0,1]
    676     testPolyBox(I, list(unit, unit));
    677 
    678     // there are no common zeros in [0,0.5]x[0,0.5]
    679     testPolyBox(I, list(unit/2, unit/2));
     759      }
     760
     761      B = bb;
     762    }
     763
     764    if (laststep)
     765    {
     766      return(1, B);
     767    }
     768  }
     769
     770  // no condition could be verified
     771  return(0, B);
     772}
     773example
     774{
     775  "EXAMPLE:"; echo = 2;
     776  ring R = 0,(x,y),dp;
     777  ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
     778
     779  interval unit = bounds(0, 1);
     780  // there may be common zeros in [0,1]x[0,1]
     781  testPolyBox(I, list(unit, unit));
     782
     783  // there are no common zeros in [0,0.5]x[0,0.5]
     784  testPolyBox(I, list(unit/2, unit/2));
    680785}
    681786
    … …  
    687792EXAMPLE: example evalIdealAtBox; evaluates an ideal at a box"
    688793{
    689     list resu;
    690 
    691     for (int j = 1; j <= size(I); j++) {
    692         resu[j] = evalPolyAtBox(I[j], B);
    693     }
    694 
    695     return(box(resu));
    696 }
    697 example
    698 {
    699     "EXAMPLE:"; echo = 2;
    700     ring R = 0,(x,y),dp;
    701     interval I1 = bounds(0, 1); I1;
    702     interval I2 = bounds(0, 1); I2;
    703 
    704     poly f = xy2 + 2x2 + (3/2)*y3x  + 1;
    705     poly g = 3x2 + 2y;
    706 
    707     ideal I = f,g;
    708     list intervals = I1,I2;
    709 
    710     evalIdealAtBox(I, intervals);
    711 }
    712 
    713 proc rootIsolationNoPreprocessing(ideal I, def start, number eps)
    714 "USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)}; @code{I ideal,
    715         B box/list of boxes, eps number};
     794  list resu;
     795
     796  for (int j = 1; j <= size(I); j++)
     797  {
     798    resu[j] = evalPolyAtBox(I[j], B);
     799  }
     800
     801  return(box(resu));
     802}
     803example
     804{
     805  "EXAMPLE:"; echo = 2;
     806  ring R = 0,(x,y),dp;
     807  interval I1 = bounds(0, 1); I1;
     808  interval I2 = bounds(0, 1); I2;
     809
     810  poly f = xy2 + 2x2 + (3/2)*y3x  + 1;
     811  poly g = 3x2 + 2y;
     812
     813  ideal I = f,g;
     814  list intervals = I1,I2;
     815
     816  evalIdealAtBox(I, intervals);
     817}
     818
     819// construct box in subring that contains var(i) iff u[i]==0
     820static proc boxProjection(box B, intvec v)
     821{
     822  box new;
     823  int i, j = 1, 1;
     824  for (i = 1; i <= size(v); i++)
     825  {
     826    if (v[i] == 0)
     827    {
     828      new = boxSet(new, j, B[i]);
     829      j++;
     830    }
     831  }
     832  return(new);
     833}
     834
     835// construct box from box in subring, replacing entries corresponding to
     836// variables that also exist in subring
     837static proc boxEmbedding(box B, box Bproj, intvec v)
     838{
     839  int N = nvars(basering);
     840  int i, j = 1, 1;
     841  for (i = 1; i<= N; i++)
     842  {
     843    if (v[i] == 0)
     844    {
     845      B = boxSet(B, i, Bproj[j]);
     846      j++;
     847    }
     848  }
     849  return(B);
     850}
     851
     852
     853// use laguerre solver for univariate case
     854static proc solveunivar(poly f){
     855   def R1 = basering;
     856   ring R = (real,30,i),x,dp;
     857   poly f = fetch(R1,f);
     858   int i,j,k;
     859   int epsilon=12;
     860   int tst;
     861   while (tst==0){
     862     if (epsilon>29){ERROR("Precision not sufficient");}
     863     setring R;
     864     list r = laguerre_solve(f,epsilon);
     865     list r1;
     866     for(j = 1;j<=size(r);j++)
     867     {
     868          if(impart(r[j])==0)
     869          {
     870            number rj =r[j];
     871            if (rj<>0){rj=round(rj*number(10)^epsilon);}           
     872            r1 = insert(r1,string(imap(R,rj)));
     873            kill rj;
     874          }
     875     }
     876     kill r;
     877     k = size(r1);
     878     setring R1;
     879     list r2;
     880     for(i=1; i<=k;i++)
     881     {
     882          execute("number m ="+r1[i]);
     883          m=m/(number(10)^epsilon);
     884          r2 = insert(r2,m);
     885          kill m;
     886     }
     887      number t = (number(10)^epsilon)^(-1);
     888      tst=1;
     889      for (j=1;j<=k;j++){
     890         if (sturm(f,r2[j]-t,r2[j]+t)>1){tst=0;break;};
     891      }
     892      epsilon++;
     893   }
     894   list iv;
     895   for (i=1; i<=size(r2);i++){
     896      number lo = r2[i]-t;
     897      number hi = r2[i]+t;
     898      iv[i]=bounds(lo,hi);
     899      kill lo,hi;
     900   }
     901return(iv);}
     902
     903proc rootIsolationNoPreprocessing(ideal I, list start, number eps)
     904"USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)};
     905        @code{I ideal, B list of boxes, eps number};
    716906ASSUME: @code{I} is a zero-dimensional radical ideal
    717907RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
    … …  
    727917        If the result of the Newton step is already contained in the interior
    728918        of the starting box, it contains a unique root.
     919NOTE:   While @code{rootIsolation} does, @code{rootIsolationNoPreprocessing}
     920        does not check input ideal for necessary conditions.
    729921EXAMPLE: example rootIsolationNoPreprocessing; exclusion test for intersection
    730922        of two ellipses"
    731923{
    732     //set of boxes smaller than size
    733     list B_size;
    734     //set of boxes which exactly contain one solution
    735     list B_star;
    736     //set of boxes initialised to input
    737     list B;
    738     if (typeof(start) == "box")
    739     {
    740         B = list(start);
    741     }
    742     else
    743     {
    744         if (typeof(start) == "list")
     924  int i;
     925  if (nvars(basering)==1){
     926    list iv = solveunivar(I[1]);
     927    list bo;
     928    for (i=1; i<=size(iv);i++){
     929        box B = list(iv[i]);
     930        bo[i]=B;
     931        kill B;
     932    }
     933    return(list(),bo);
     934  }
     935  // set of boxes smaller than eps
     936  list lBoxesSize;
     937  // set of boxes which exactly contain one solution
     938  list lBoxesUnique;
     939  // set of boxes initialised to input
     940  list lBoxes = start;
     941  //help set of boxes
     942  list lBoxesHelp;
     943
     944  int i, j, s, cnt;
     945  int zeroTest;
     946  int N = nvars(basering);
     947
     948  intvec empty;
     949  def rSource = basering;
     950  list lSubstBoxesSize, lSubstBoxesUnique;
     951
     952  ideal Isubst;
     953  string rSubstString;
     954
     955  while (size(lBoxes) <> 0)
     956  {
     957    // lBoxesHelp is empty set
     958    lBoxesHelp = list();
     959    s = 0;
     960
     961    for (i=1; i<=size(lBoxes); i++)
     962    {
     963      //case that maybe there is a root in the box
     964      zeroTest, lBoxes[i] = testPolyBox(I,lBoxes[i]);
     965      if (zeroTest == 1)
     966      {
     967        lBoxesUnique[size(lBoxesUnique)+1] = lBoxes[i];
     968      }
     969      if (zeroTest == 0)
     970      {
     971        // check for empty interior
     972        empty = boxEmptyIntervals(lBoxes[i]);
     973        // check if at least one interval is single number
     974        if (max(empty[1..N]))
    745975        {
    746             // add sanity check
    747             B = start;
     976          // If the box lBoxes[i] contains singleton intervals, i.e. single
     977          // numbers, we substitute the corresponding ring variables with these
     978          // numbers and pass to the subring where these variables have been
     979          // removed.
     980          // As rootIsolation now accpets ideals with arbitrary generators we
     981          // compute the zeros of the substituion ideal and reassemble the
     982          // boxes afterwards.
     983          //
     984          // This is probably overkill.
     985          Isubst = I;
     986          // name ring rSubst(voice) to avoid name conflicts if
     987          // we (unlikely) enter another level of recursion
     988          rSubstString = "ring rSubst(" + string(voice) + ") = 0,(";
     989          for (j = 1; j <= N; j++)
     990          {
     991            if (empty[j])
     992            {
     993              Isubst = subst(Isubst, var(j), lBoxes[i][j][1]);
     994            }
     995            else
     996            {
     997              rSubstString = rSubstString + varstr(j) + ",";
     998            }
     999          }
     1000          rSubstString = string(rSubstString[1..(size(rSubstString)-1)]) + "),dp;";
     1001          dbprint("passing to subring");
     1002          execute(rSubstString);
     1003          ideal Isubst = imap(rSource, Isubst);
     1004          number eps = imap(rSource, eps);
     1005          lSubstBoxesSize, lSubstBoxesUnique = rootIsolation(Isubst,
     1006              boxProjection(lBoxes[i], empty), eps);
     1007          setring rSource;
     1008          for (j = 1; j <= size(lSubstBoxesSize); j++)
     1009          {
     1010            lSubstBoxesSize[j] = boxEmbedding(lBoxes[i], lSubstBoxesSize[j], empty);
     1011          }
     1012          lBoxesSize = lBoxesSize + lSubstBoxesSize;
     1013          for (j = 1; j <= size(lSubstBoxesUnique); j++)
     1014          {
     1015            lSubstBoxesUnique[j] = boxEmbedding(lBoxes[i], lSubstBoxesUnique[j], empty);
     1016          }
     1017          lBoxesUnique = lBoxesUnique + lSubstBoxesUnique;
     1018          kill rSubst(voice);
    7481019        }
    7491020        else
    7501021        {
    751             ERROR("second arg must be box or list");
     1022          // remove boxes smaller than limit eps
     1023          if (boxLength(lBoxes[i]) < eps)
     1024          {
     1025            lBoxesSize[size(lBoxesSize)+1] = lBoxes[i];
     1026          }
     1027          else
     1028          {
     1029            // else split the box and put the smaller boxes to lBoxesHelp
     1030            lBoxesHelp[s+1..s+2] = splitBox(lBoxes[i], I);
     1031            s = s+2;
     1032          }
    7521033        }
    753     }
    754     //help set of boxes
    755     list B_prime;
    756 
    757     list split;
    758     int i, s;
    759     int zeroTest;
    760 
    761     // debug
    762     //int cnt;
    763 
    764     while (size(B) <> 0)
    765     {
    766         // B_prime is empty set
    767         B_prime = list();
    768         s = 0;
    769 
    770         for (i=1; i<=size(B); i++)
    771         {
    772             //case that maybe there is a root in the box
    773             zeroTest, B[i] = testPolyBox(I,B[i]);
    774 
    775             // maybe refine boxes in Bstar in later steps
    776             if (zeroTest == 1)
    777             {
    778                 B_star[size(B_star)+1] = B[i];
    779             }
    780             if (zeroTest == 0)
    781             {
    782                 // case that box is smaller than the input limit eps
    783                 if (lengthBox(B[i]) < eps)
    784                 {
    785                     B_size[size(B_size)+1] = B[i];
    786                 }
    787                 else
    788                 {
    789                     // else split the box and put the smaller boxes to B_prime
    790                     B_prime[s+1..s+2] = splitBox(B[i], I);
    791                     s = s+2;
    792                 }
    793             }
    794         }
    795 
    796         // debug
    797         //cnt++;
    798         //print(string(cnt, " ", s, " ", int(memory(0)/1024), "k"));
    799 
    800         B = B_prime;
    801     }
    802     return(B_size, B_star);
    803 }
    804 example
    805 {
    806     "EXAMPLE:"; echo = 2;
    807 
    808     ring R = 0,(x,y),dp;
    809     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
    810 
    811     interval i = bounds(-3/2,3/2);
    812     box B = list(i, i);
    813 
    814     list result = rootIsolationNoPreprocessing(I, B, 1/512);
    815     size(result[1]);
    816     size(result[2]);
    817 
    818     result;
    819 }
    820 
    821 static proc noRootsOnBoundary(ideal I, box B)
    822 "USAGE:  noZeroesOnBoundary(I, B), I ideal, B box
    823 RETURN: intvec iv where:
    824         iv[i] == 1 if there is no root of I on the boundary of B[i], 0 else
    825 EXAMPLE: example noRootsOnBoundary; tests boxes for roots"
    826 {
    827     int N = nvars(basering);
    828     int i, j, k;
    829     intvec noZero = 0:(2*N);
    830 
    831     box evalIntersection;
    832 
    833     for (i = 1; i <= N; i++)
    834     {
    835         for (j = 1; j <= 2; j++)
    836         {
    837             evalIntersection = evalIdealAtBox(I, boxSet(B, i, bounds(B[i][j])));
    838             for (k = 1; k <= N; k++)
    839             {
    840                 if (evalIntersection[k][1]*evalIntersection[k][2] > 0)
    841                 {
    842                     noZero[2*i+j-2] = 1;
    843                     break;
    844                 }
    845             }
    846             if (!noZero[2*i+j-2])
    847             {
    848                 // check if V(I + ...) is empty over CC[x(...)]
    849                 noZero[2*i+j-2] = std(I + (var(i) - B[i][j])) == 1;
    850             }
    851         }
    852     }
    853     return(noZero);
    854 }
    855 example
    856 {
    857     "EXAMPLE:"; echo = 2;
    858     ring R = 0,(x,y),dp;
    859 
    860     ideal I = x+1,y-2;
    861 
    862     interval I1, I2, I3 = bounds(0,1), bounds(-1,0), bounds(-2,2);
    863 
    864     noRootsOnBoundary(I, box(list(I1,I2)));
    865     noRootsOnBoundary(I, box(list(I2,I1)));
    866     noRootsOnBoundary(I, box(list(I2,I3)));
    867 }
    868 
    869 //im Moment geht das nur mit eingegebener eliminationsordnung
    870 proc rootIsolation(ideal I, box start, number eps)
    871 "USAGE:  @code{rootIsolation(I, start, eps)}; @code{I ideal, start box,
    872         eps number}
     1034      }
     1035    }
     1036
     1037    cnt++;
     1038    dbprint(string("[", cnt, "] ", s, " boxes"));
     1039
     1040    lBoxes = lBoxesHelp;
     1041  }
     1042  return(lBoxesSize, lBoxesUnique);
     1043}
     1044example
     1045{
     1046  "EXAMPLE:"; echo = 2;
     1047
     1048  ring R = 0,(x,y),dp;
     1049  ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
     1050
     1051  interval i = bounds(-3/2,3/2);
     1052  box B = list(i, i);
     1053
     1054  list result = rootIsolationNoPreprocessing(I, list(B), 1/512);
     1055  size(result[1]);
     1056  size(result[2]);
     1057
     1058  result;
     1059}
     1060
     1061// this takes a triangular decomposition of an ideal and applies rootIsolation
     1062// to all the ideals in the list and checks afterwards if the found boxes
     1063// intersect.
     1064static proc rootIsolationTriangL(list T, #)
     1065{
     1066  int n = size(T);
     1067  int i, j, k, l;
     1068  list lBoxesUnique, lBoxesSize, lRoots;
     1069  def bInt;
     1070  for (i = 1; i <= n; i++)
     1071  {
     1072    lRoots = rootIsolation(T[i], #);
     1073    // since we don't know anything about boxes in the first list just put them
     1074    // together
     1075    lBoxesSize = lBoxesSize + lRoots[1];
     1076    lBoxesUnique = lBoxesUnique + lRoots[2];
     1077  }
     1078  return(lBoxesSize, lBoxesUnique);
     1079}
     1080
     1081proc rootIsolationPrimdec(ideal I)
     1082"USAGE:  @code{rootIsolationPrimdec(I)};
     1083        @code{I ideal}
     1084ASSUME: @code{I} is a zero-dimensional radical ideal
     1085RETURN: @code{L}, where @code{L}
     1086        contains boxes which contain exactly one element of V(@code{I})
     1087PURPOSE: same as @code{rootIsolation}, but speeds up computation and improves output
     1088        by doing a primary decomposition before doing the root isolation
     1089THEORY: For the primary decomposition we use the algorithm of Gianni-Traeger-Zarcharias.
     1090NOTE:   This algorithm and some procedures used therein perform Groebner basis
     1091        computations in @code{basering}. It is thus advised to define @code{I}
     1092        w.r.t. a fast monomial ordering. @*
     1093EXAMPLE: example rootIsolationPrimdec; for intersection of two ellipses"
     1094{
     1095  def R = basering;
     1096  int n = nvars(R);
     1097  I=radical(I);
     1098  list dc = primdecGTZ(I);
     1099  list sols,va,dci,ri,RL,RLi;
     1100  int i,j,l;
     1101  ideal idi;
     1102  interval ii;
     1103  RL = ringlist(R);
     1104  for (i = 1; i<=size(dc);i++)
     1105  {
     1106    dci = elimpart(dc[i][1]);
     1107    idi = dci[1];
     1108    list keepvar,delvar;
     1109    for (j=1; j<=n;j++)
     1110    {
     1111      if ((dci[4][j]!=0) or (deg(dci[5][j])>0))
     1112      {
     1113        keepvar[size(keepvar)+1] = j;
     1114      } else {
     1115        delvar[size(delvar)+1] = j;
     1116      }
     1117      if ((dci[4][j]==0) and (deg(dci[5][j])>0))
     1118      {
     1119        idi = idi + ideal(var(j)-dci[5][j]);
     1120      }
     1121    }
     1122    RLi = RL;
     1123    if (size(delvar)>0){
     1124      RLi[2]=deleteSublist(intvec(delvar[1..size(delvar)]),RLi[2]);
     1125      RLi[3]=list(list("dp",intvec(1:size(keepvar))),list("C",0));
     1126    }
     1127    def Ri = ring(RLi);
     1128    setring Ri;
     1129    ideal idi=imap(R,idi);
     1130    ri = rootIsolation(idi);
     1131    setring R;
     1132    kill Ri;
     1133    list ivlist;
     1134    for (l=1;l<=size(ri[2]);l++)
     1135    {
     1136      for (j=1; j<=size(keepvar); j++)
     1137      {
     1138        ivlist[keepvar[j]]=ri[2][l][j];
     1139      }
     1140      for (j=1; j<=size(delvar); j++)
     1141      {
     1142        ii = bounds(number(dci[5][delvar[j]]));
     1143        ivlist[delvar[j]]=ii;
     1144      }
     1145      sols[size(sols)+1]=box(ivlist);
     1146    }
     1147    kill keepvar,delvar,ivlist;
     1148  }
     1149  return(sols);
     1150}
     1151example
     1152{
     1153  "EXAMPLE:"; echo = 2;
     1154
     1155  ring R = 0,(x,y),dp;
     1156  ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
     1157  list result = rootIsolationPrimdec(I);
     1158
     1159  result;
     1160}
     1161
     1162
     1163proc rootIsolation(ideal I, list #)
     1164"USAGE:  @code{rootIsolation(I, [start, eps])};
     1165        @code{I ideal, start box, eps number}
    8731166ASSUME: @code{I} is a zero-dimensional radical ideal
    8741167RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
    … …  
    8861179NOTE:   This algorithm and some procedures used therein perform Groebner basis
    8871180        computations in @code{basering}. It is thus advised to define @code{I}
    888         w.r.t. a fast monomial ordering.
     1181        w.r.t. a fast monomial ordering. @*
     1182        The algorithm performs checks on @code{I} to prevent errors. If
     1183        @code{I} does not have the right number of generators, we first try to
     1184        find a suitable Groebner basis. If this fails we apply the algorithm
     1185        to the triangular decomposition of @code{I}.
    8891186EXAMPLE: example rootIsolation; for intersection of two ellipses"
    8901187{
    891     int N = nvars(basering);
    892     int i, j, k, l;
    893 
    894     intvec noZeroes;
    895     // check if there are roots on the boundary of start
    896     while(1)
    897     {
    898         noZeroes = noRootsOnBoundary(I, start);
    899         // stop if all boundaries root-free
    900         if (product(noZeroes)) { break; }
    901 
    902         for (i = 1; i <= N; i++)
     1188  //voice;
     1189  //printlevel=voice;
     1190  dbprint("Start rootisolation");
     1191  int N = nvars(basering);
     1192  int i, j, k, l;
     1193  int genNumIndex = 0;
     1194
     1195  // input parameter check
     1196  int determinebox;
     1197
     1198  if (size(#)==0){
     1199    determinebox = 1;
     1200    number eps = 0;
     1201  }
     1202
     1203  if (size(#)==1){
     1204    if (typeof(#[1])=="box"){
     1205        box start = #[1];
     1206        number eps = 0;
     1207    }
     1208    if (typeof(#[1])=="number"){
     1209        number eps = #[1];
     1210        determinebox = 1;
     1211    }
     1212    if ((typeof(#[1])!="box") and (typeof(#[1])!="number")){
     1213        ERROR("second argument should be a box or a number");
     1214    }
     1215  }
     1216
     1217  if (size(#)==2)
     1218  {
     1219    if (typeof(#[1])!="box")
     1220    {
     1221      ERROR("second argument should be a box");
     1222    }
     1223    if (typeof(#[2])!="number")
     1224    {
     1225      ERROR("third argument should be a number");
     1226    }
     1227    box start = #[1];
     1228    number eps = #[2];
     1229  }
     1230
     1231  if (size(#)>2){
     1232    ERROR("too many arguments");
     1233  }
     1234
     1235  // compute reduced GB in (hopefully) fast ordering
     1236  option(redSB);
     1237  ideal fastGB = groebner(I);
     1238
     1239  if (dim(fastGB) > 0)
     1240  {
     1241    ERROR("rootIsolation: ideal not zero-dimensional");
     1242  }
     1243
     1244  ideal rad = radical(fastGB);
     1245  // since fastGB is contained in rad we only need to test one inclusion size
     1246  // counts only non-zero entries
     1247  if (size(reduce(rad, fastGB)) > 0)
     1248  {
     1249    "Warning: ideal not radical, passing to radical.";
     1250    I = rad;
     1251    fastGB = groebner(I);
     1252  }
     1253
     1254  // create lp-rings
     1255  ring rSource = basering;
     1256  list rList = ringlist(rSource);
     1257  // set ordering to lp
     1258  rList[3] = list( list( "lp", 1:N ), list( "C", 0 ) );
     1259  // save variable order for later
     1260  list varstrs = rList[2];
     1261
     1262  for (i = 1; i <= N; i++)
     1263  {
     1264    // permute variables in ringlist, create and switch to ring with
     1265    // elimination ordering
     1266    rList[2][i] = varstrs[N];
     1267    rList[2][N] = varstrs[i];
     1268    def rElimOrd(i) = ring(rList);
     1269    setring rElimOrd(i);
     1270    // get GB in elimination ordering, note that GB[1] only depends on var(i)
     1271    // of rSource, which is var(N) in rElimOrd(i)
     1272    ideal GB = fglm(rSource, fastGB);
     1273    if (size(GB) == N)
     1274    {
     1275      genNumIndex = i;
     1276    }
     1277    setring rSource;
     1278    rList[2] = varstrs;
     1279  }
     1280
     1281  if (N <> ncols(I))
     1282  {
     1283    if (genNumIndex > 0)
     1284    {
     1285      // If elements of V(I) have pairwise distinct xj coordinates for some
     1286      // variable xj, then the reduced Groebner basis w.r.t. the monomial
     1287      // ordering where xj is smaller than all other variables has exactly
     1288      // nvars(basering) generators, see [2].
     1289      // As we compute all these Groebner bases anyway we can replace the input
     1290      // ideal with a proper generator set if necessary.
     1291      I = imap(rElimOrd(genNumIndex), GB);
     1292      dbprint("Replaced ideal with suitable reduced Groebner basis.");
     1293    }
     1294    else
     1295    {
     1296      // the ideals in a triangular decomposition always satisfy the conditions
     1297      // on I, so we apply rootIsolation to them as a last resort
     1298
     1299      // note that zero-dimensional ideals can always be generated by
     1300      // nvars(basering) polynomials, however the way of computing these
     1301      // generators is currently not known to the authors
     1302      setring rElimOrd(N);
     1303      list T = triangL(GB);
     1304      setring rSource;
     1305      list T = imap(rElimOrd(N), T);
     1306      // kill elimination rings to avoid name conflicts in rootIsolation called
     1307      // later
     1308      for (i = 1; i <= N; i++)
     1309      {
     1310        kill rElimOrd(i);
     1311      }
     1312      dbprint("Applying rootIsolation to triangular decomposition.");
     1313      return(rootIsolationTriangL(T, #));
     1314    }
     1315  }
     1316
     1317  // need at least two variables
     1318  if (N < 2)
     1319  {
     1320    if (determinebox)
     1321    {
     1322      number mx = maxabs(fastGB[1]);
     1323      box start = list(bounds(-mx, mx));
     1324    }
     1325    return(rootIsolationNoPreprocessing(I, list(start), eps));
     1326  }
     1327
     1328  // need nvars(basering) generators from here on
     1329  rList = ringlist(rSource);
     1330  // first construct univariate ring
     1331  int numords = size(rList[3]);
     1332  // remove all but first variables
     1333  rList[2] = list(rList[2][1]);
     1334  // change ordering accordingly (keep last block)
     1335  rList[3] = list( list(rList[3][1][1], intvec(1)), rList[3][numords] );
     1336
     1337  // construct and switch to univariate ring
     1338  def rUnivar = ring(rList);
     1339  setring rUnivar;
     1340
     1341  // some necessary variables
     1342  ideal gbUnivar;
     1343  // maps var(N) in rElimOrd(i) to var(1) in rUnivar
     1344  intvec varmap = (0:(N-1)),1;
     1345  number eps = fetch(rSource, eps);
     1346  list univarResult, startBoxesPerDim;
     1347
     1348  intvec sizes, repCount = 0:N, 1:N;
     1349  int numBoxes = 1;
     1350
     1351  number lo, up;
     1352  number mx;
     1353
     1354  for (i = 1; i <= N; i++)
     1355  {
     1356    dbprint("variable ",i);
     1357    setring rElimOrd(i);
     1358    // throw out non-univariate polys
     1359    GB = GB[1];
     1360
     1361    setring rUnivar;
     1362    gbUnivar = fetch(rElimOrd(i), GB, varmap);
     1363    // clean up ring and its elements
     1364    kill rElimOrd(i);
     1365
     1366    // check if there are roots on the bounds of start[i]
     1367    // (very easy in univariate case)
     1368    if (determinebox)
     1369    {
     1370      dbprint("univariate GB", gbUnivar,laguerre_solve(gbUnivar[1],5));
     1371      mx = maxabs(std(gbUnivar)[1]);
     1372      dbprint("maxabs ",mx);
     1373      lo, up = -mx, mx;
     1374    }
     1375    else
     1376    {
     1377      lo, up = start[i][1], start[i][2];
     1378    }
     1379    // move bounds if they are a root of the polynomial
     1380    while(subst(gbUnivar[1], var(1), lo) == 0)
     1381    {
     1382      lo = lo - 1;
     1383    }
     1384    while(subst(gbUnivar[1], var(1), up) == 0)
     1385    {
     1386      up = up + 1;
     1387    }
     1388
     1389    // get boxes containing roots of univariate polynomial
     1390    //"rootIsolationNoPreprocessing ";
     1391    univarResult = rootIsolationNoPreprocessing(gbUnivar,
     1392      list(box(list(bounds(lo, up)))), eps);
     1393    // maybe result[1] is not empty, so take both
     1394    startBoxesPerDim[i] = univarResult[1] + univarResult[2];
     1395    // debug:
     1396    dbprint(string("Sieved variable ", varstrs[i], " to ",
     1397      size(startBoxesPerDim[i]), " intervals."));
     1398
     1399    sizes[i] = size(startBoxesPerDim[i]);
     1400    numBoxes = numBoxes * sizes[i];
     1401
     1402    // stop early if one variable already has no roots
     1403    if (numBoxes == 0)
     1404    {
     1405      dbprint("No roots exist for input. Stopping.");
     1406      return(list(), list());
     1407    }
     1408  }
     1409
     1410  setring rSource;
     1411
     1412  for (i = N-1; i >= 1; i--)
     1413  {
     1414    repCount[i] = repCount[i+1] * sizes[i+1];
     1415  }
     1416
     1417  list startBoxes, sbTemp;
     1418
     1419  // prepare a list of lists
     1420  for (i = 1; i <= numBoxes; i++)
     1421  {
     1422    sbTemp[i] = list();
     1423  }
     1424
     1425  // computes "cartesian product" of found boxes to lift to N variables
     1426  for (i = 1; i <= N; i++)
     1427  {
     1428    // full loop of elements
     1429    for (j = 0; j < numBoxes; j = j + sizes[i]*repCount[i])
     1430    {
     1431      // boxes
     1432      for (k = 0; k < sizes[i]; k++)
     1433      {
     1434        // repetions
     1435        for (l = 1; l <= repCount[i]; l++)
    9031436        {
    904             for (j = 1; j <= 2; j++)
    905             {
    906                 // change offending boundary
    907                 if (!noZeroes[2*i+j-2])
    908                 {
    909                     start = boxSet(start, i,
    910                         start[i] + (-1)^j * bounds(0, length(start[i])/10));
    911                 }
    912             }
     1437          // last index since we need interval in one-dimensional box
     1438          sbTemp[j+k*repCount[i]+l][i] = startBoxesPerDim[i][k+1][1];
    9131439        }
    914     }
    915 
    916     // need at least two variables
    917     if (N < 2)
    918     {
    919         return(rootIsolationNoPreprocessing(I, start, eps));
    920     }
    921 
    922     // compute reduced GB in (hopefully) fast ordering
    923     option(redSB);
    924     ideal fastGB = groebner(I);
    925 
    926     ring rSource = basering;
    927     list rList = ringlist(rSource);
    928     // first construct univariate ring
    929     int numords = size(rList[3]);
    930     // remove all but first variables
    931     rList[2] = list(rList[2][1]);
    932     // change ordering accordingly (keep last block)
    933     rList[3] = list( list(rList[3][1][1], intvec(1)), rList[3][numords] );
    934 
    935     // construct and switch to univariate ring
    936     ring rUnivar = ring(rList);
    937 
    938     // some necessary variables
    939     ideal gbUnivar;
    940     // maps var(N) in rElimOrd(i) to var(1) in rUnivar
    941     intvec varmap = (0:(N-1)),1;
    942     number eps = fetch(rSource, eps);
    943     list univarResult, startBoxesPerDim;
    944 
    945     // now prepare lp-ring
    946     setring rSource;
    947     rList = ringlist(rSource);
    948     // set ordering to lp
    949     rList[3] = list( list( "lp", 1:N ), list( "C", 0 ) );
    950     // save variable order for later
    951     list varstrs = rList[2];
    952 
    953     for (i = 1; i <= N; i++)
    954     {
    955         // permute variables in ringlist,
    956         // create and switch to ring with elimination ordering
    957         rList[2][i] = varstrs[N];
    958         rList[2][N] = varstrs[i];
    959         ring rElimOrd = ring(rList);
    960         // get GB in elimination ordering, note that GB[1] only depends on
    961         // var(i) of rSource, which is var(N) in rElimOrd
    962         ideal GB = fglm(rSource, fastGB);
    963         GB = GB[1];
    964 
    965         setring rUnivar;
    966         gbUnivar = fetch(rElimOrd, GB, varmap);
    967         // clean up ring and its elements
    968         kill rElimOrd;
    969 
    970         // get boxes containing roots of univariate polynomial
    971         univarResult = rootIsolationNoPreprocessing(gbUnivar,
    972             box(list(start[i])), eps);
    973         // maybe result[1] is not empty, so take both
    974         startBoxesPerDim[i] = univarResult[1] + univarResult[2];
    975         // debug:
    976         print(string("Sieved variable ", varstrs[i], " to ",
    977             size(startBoxesPerDim[i]), " intervals."));
    978 
    979         // reset ring variable order
    980         setring rSource;
    981         rList[2] = varstrs;
    982     }
    983 
    984     intvec sizes, repCount = 0:N, 1:N;
    985     int numBoxes = 1;
    986 
    987     for (i = 1; i <= N; i++)
    988     {
    989         sizes[i] = size(startBoxesPerDim[i]);
    990         numBoxes = numBoxes * sizes[i];
    991     }
    992 
    993     for (i = N-1; i >= 1; i--)
    994     {
    995         repCount[i] = repCount[i+1] * sizes[i+1];
    996     }
    997 
    998     list startBoxes, sbTemp;
    999 
    1000     // prepare a list of lists
    1001     for (i = 1; i <= numBoxes; i++)
    1002     {
    1003         sbTemp[i] = list();
    1004     }
    1005 
    1006     // computes "cartesian product" of found boxes to lift to N variables
    1007     for (i = 1; i <= N; i++)
    1008     {
    1009         // full loop of elements
    1010         for (j = 0; j < numBoxes; j = j + sizes[i]*repCount[i])
    1011         {
    1012             // boxes
    1013             for (k = 0; k < sizes[i]; k++)
    1014             {
    1015                 // repetions
    1016                 for (l = 1; l <= repCount[i]; l++)
    1017                 {
    1018                     // last index since we need interval in one-dimensional box
    1019                     sbTemp[j+k*repCount[i]+l][i] = startBoxesPerDim[i][k+1][1];
    1020                 }
    1021             }
    1022         }
    1023     }
    1024 
    1025     // since we're back in rSource box(...) will return box of proper size
    1026     for (i = 1; i <= size(sbTemp); i++)
    1027     {
    1028         startBoxes[i] = box(sbTemp[i]);
    1029     }
    1030 
    1031     return(rootIsolationNoPreprocessing(I, startBoxes, eps));
    1032 }
    1033 example
    1034 {
    1035     "EXAMPLE:"; echo = 2;
    1036 
    1037     ring R = 0,(x,y),dp;
    1038     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
    1039 
    1040     interval i = bounds(-3/2,3/2);
    1041     box B = list(i, i);
    1042 
    1043     list result = rootIsolation(I, B, 0);
    1044 
    1045     result;
     1440      }
     1441    }
     1442  }
     1443
     1444  // since we're back in rSource box(...) will return box of proper size
     1445  for (i = 1; i <= size(sbTemp); i++)
     1446  {
     1447    startBoxes[i] = box(sbTemp[i]);
     1448  }
     1449  // can be optimized
     1450  return(rootIsolationNoPreprocessing(I, startBoxes, eps));
     1451}
     1452example
     1453{
     1454  "EXAMPLE:"; echo = 2;
     1455
     1456  ring R = 0,(x,y),dp;
     1457  ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
     1458
     1459  interval i = bounds(-3/2,3/2);
     1460  box B = list(i, i);
     1461
     1462  list result = rootIsolation(I, B);
     1463
     1464  result;
    10461465}
    10471466// vim: ft=singular
Note: See TracChangeset for help on using the changeset viewer.