[eeb133] | 1 | // E. Tobis 12.Nov.2004 |
---|
| 2 | // last change 16. Apr. 2005 (G.-M. Greuel) |
---|
| 3 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 4 | category="Symbolic-numerical solving" |
---|
| 5 | info=" |
---|
| 6 | LIBRARY: signcond.lib Routines for computing realizable sign conditions |
---|
| 7 | AUTHOR: Enrique A. Tobis, etobis@dc.uba.ar |
---|
| 8 | |
---|
| 9 | PROCEDURES: |
---|
| 10 | signcnd(P,I) The sign conditions realized the polynomials of P on a V(I) |
---|
| 11 | psigncnd(P,l) Pretty prints the output of signcnd (l) |
---|
| 12 | firstoct(P) The number of elements of a V(I) with every coordinate > 0 |
---|
| 13 | |
---|
| 14 | KEYWORDS: real roots,sign conditions |
---|
| 15 | "; |
---|
| 16 | |
---|
| 17 | LIB "mrrcount.lib"; |
---|
| 18 | LIB "linalg.lib"; |
---|
| 19 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 20 | |
---|
| 21 | proc firstoct(ideal I) |
---|
| 22 | "USAGE: firstoct(i); i ideal |
---|
| 23 | RETURN: number: the number of points of V(i) lying in the first orthant |
---|
| 24 | ASSUME: i is a Groebner basis |
---|
| 25 | SEE ALSO: signcnd |
---|
| 26 | EXAMPLE: example firstoct; shows an example" |
---|
| 27 | { |
---|
| 28 | ideal firstoctant; |
---|
| 29 | int j; |
---|
| 30 | list result; |
---|
| 31 | int n; |
---|
| 32 | |
---|
| 33 | if (isparam(I)) { |
---|
| 34 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 35 | } |
---|
| 36 | |
---|
| 37 | for (j = nvars(basering);j > 0;j--) { |
---|
| 38 | firstoctant = firstoctant + var(j); |
---|
| 39 | } |
---|
| 40 | |
---|
| 41 | result = signcnd(firstoctant,I); |
---|
| 42 | |
---|
| 43 | list fst; |
---|
| 44 | for (j = nvars(basering);j > 0;j--) { |
---|
| 45 | fst[j] = 1; |
---|
| 46 | } |
---|
| 47 | |
---|
| 48 | n = isIn(fst,result[1]); |
---|
| 49 | |
---|
| 50 | if (n != -1) { |
---|
| 51 | return (result[2][n]); |
---|
| 52 | } else { |
---|
| 53 | return (0); |
---|
| 54 | } |
---|
| 55 | } |
---|
| 56 | example |
---|
| 57 | { |
---|
| 58 | echo = 2; |
---|
| 59 | ring r = 0,(x,y),dp; |
---|
| 60 | ideal i = (x-2)*(x+3)*x,y*(y-1); |
---|
| 61 | firstoct(i); |
---|
| 62 | } |
---|
| 63 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 64 | |
---|
| 65 | proc signcnd(ideal P,ideal I) |
---|
| 66 | "USAGE: signcnd(P,i); ideal P,i. i must be a grobner basis |
---|
| 67 | RETURN: list: the sign conditions realized by the polynomials of P on V(I). |
---|
| 68 | See the example for an explanation of the output. |
---|
| 69 | SEE ALSO: firstoct |
---|
| 70 | EXAMPLE: example signcnd; shows an example" |
---|
| 71 | { |
---|
| 72 | ideal B; |
---|
| 73 | |
---|
| 74 | // Cumulative stuff |
---|
| 75 | matrix M; |
---|
| 76 | matrix SQs; |
---|
| 77 | matrix C; |
---|
| 78 | list Signs; |
---|
| 79 | list Exponents; |
---|
| 80 | |
---|
| 81 | // Used to store the precalculated SQs |
---|
| 82 | list SQvalues; |
---|
| 83 | list SQpositions; |
---|
| 84 | |
---|
| 85 | int i; |
---|
| 86 | |
---|
| 87 | // Variables for each step |
---|
| 88 | matrix Mi; |
---|
| 89 | matrix M3x3[3][3]; |
---|
| 90 | matrix M3x3inv[3][3]; // Constant matrices |
---|
| 91 | matrix c[3][1]; |
---|
| 92 | matrix sq[3][1]; |
---|
| 93 | int j; |
---|
| 94 | list exponentsi; |
---|
| 95 | list signi; |
---|
| 96 | int numberOfNonZero; |
---|
| 97 | |
---|
| 98 | if (isparam(P) || isparam(I)) { |
---|
| 99 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 100 | } |
---|
| 101 | |
---|
| 102 | M3x3 = matrix(1,3,3); |
---|
| 103 | M3x3 = 1,1,1,0,1,-1,0,1,1; // The 3x3 matrix |
---|
| 104 | M3x3inv = inverse(M3x3); |
---|
| 105 | |
---|
| 106 | // First, we compute sturmquery(1,V(I)) |
---|
| 107 | I = groebner(I); |
---|
| 108 | B = qbase(I); |
---|
| 109 | sq[1,1] = sturmquery(1,B,I); // Number of real roots in V(I) |
---|
| 110 | SQvalues = SQvalues + list(sq[1,1]); |
---|
| 111 | SQpositions = SQpositions + list(1); |
---|
| 112 | |
---|
| 113 | // We initialize the cumulative variables |
---|
| 114 | M = matrix(1,1,1); |
---|
| 115 | Exponents = list(list()); |
---|
| 116 | Signs = list(list()); |
---|
| 117 | |
---|
| 118 | i = 1; |
---|
| 119 | |
---|
| 120 | while (i <= size(P)) { // for each poly in P |
---|
| 121 | |
---|
| 122 | sq[2,1] = sturmquery(P[i],B,I); |
---|
| 123 | sq[3,1] = sturmquery(P[i]^2,B,I); |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | c = M3x3inv*sq; |
---|
| 127 | |
---|
| 128 | // We have to eliminate the 0 elements in c |
---|
| 129 | exponentsi = list(); |
---|
| 130 | signi = list(); |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | // We determine the list of signs which correspond to a nonzero |
---|
| 134 | // number of roots |
---|
| 135 | numberOfNonZero = 3; |
---|
| 136 | |
---|
| 137 | if (c[1,1] != 0) { |
---|
| 138 | signi = list(0); |
---|
| 139 | } else { |
---|
| 140 | numberOfNonZero--; |
---|
| 141 | } |
---|
| 142 | |
---|
| 143 | if (c[2,1] != 0) { |
---|
| 144 | signi = signi + list(1); |
---|
| 145 | } else { |
---|
| 146 | numberOfNonZero--; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | if (c[3,1] != 0) { |
---|
| 150 | signi = signi + list(-1); |
---|
| 151 | } else { |
---|
| 152 | numberOfNonZero--; |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | // We now determine the little matrix we'll work with, |
---|
| 156 | // and the list of exponents |
---|
| 157 | if (numberOfNonZero == 3) { |
---|
| 158 | Mi = M3x3; |
---|
| 159 | exponentsi = list(0,1,2); |
---|
| 160 | } else {if (numberOfNonZero == 2) { |
---|
| 161 | Mi = matrix(1,2,2); |
---|
| 162 | Mi[1,2] = 1; |
---|
| 163 | if (c[1,1] != 0 && c[2,1] != 0) { // 0,1 |
---|
| 164 | Mi[2,1] = 0; |
---|
| 165 | Mi[2,2] = 1; |
---|
| 166 | } else {if (c[1,1] != 0 && c[3,1] != 0) { // 0,-1 |
---|
| 167 | Mi[2,1] = 0; |
---|
| 168 | Mi[2,2] = -1; |
---|
| 169 | } else { // 1,-1 |
---|
| 170 | Mi[2,1] = 1; |
---|
| 171 | Mi[2,2] = -1; |
---|
| 172 | }} |
---|
| 173 | exponentsi = list(0,1); |
---|
| 174 | } else {if (numberOfNonZero == 1) { |
---|
| 175 | Mi = matrix(1,1,1); |
---|
| 176 | exponentsi = list(0); |
---|
| 177 | }}} |
---|
| 178 | |
---|
| 179 | // We store the Sturm Queries we'll need later |
---|
| 180 | if (numberOfNonZero == 2) { |
---|
| 181 | SQvalues = SQvalues + list(sq[2,1]); |
---|
| 182 | SQpositions = SQpositions + list(size(Exponents)+1); |
---|
| 183 | } else {if (numberOfNonZero == 3) { |
---|
| 184 | SQvalues = SQvalues + list(sq[2,1],sq[3,1]); |
---|
| 185 | SQpositions = SQpositions + list(size(Exponents)+1,size(Exponents)*2+1); |
---|
| 186 | }} |
---|
| 187 | |
---|
| 188 | // Now, we accumulate information |
---|
| 189 | M = tensor(Mi,M); |
---|
| 190 | Signs = expprod(Signs,signi); |
---|
| 191 | Exponents = expprod(Exponents,exponentsi); |
---|
| 192 | |
---|
| 193 | i++; |
---|
| 194 | } |
---|
| 195 | |
---|
| 196 | // At this point, we have the cumulative matrix, |
---|
| 197 | // the vector of exponents and the matching sign conditions. |
---|
| 198 | // We have to solve the big linear system to finish. |
---|
| 199 | |
---|
| 200 | M = inverse(M); |
---|
| 201 | |
---|
| 202 | // We have to compute the constants vector (the Sturm Queries) |
---|
| 203 | |
---|
| 204 | SQs = matrix(1,size(Exponents),1); |
---|
| 205 | |
---|
| 206 | j = 1; // We'll iterate over the presaved SQs |
---|
| 207 | |
---|
| 208 | for (i = 1;i <= size(Exponents);i++) { |
---|
| 209 | if (j <= size(SQvalues)) { |
---|
| 210 | if (SQpositions[j] == i) { |
---|
| 211 | SQs[i,1] = SQvalues[j]; |
---|
| 212 | j++; |
---|
| 213 | } else { |
---|
| 214 | SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I); |
---|
| 215 | } |
---|
| 216 | } else { |
---|
| 217 | SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I); |
---|
| 218 | } |
---|
| 219 | } |
---|
| 220 | |
---|
| 221 | C = M*SQs; |
---|
| 222 | |
---|
| 223 | list result; |
---|
| 224 | result[2] = list(); |
---|
| 225 | result[1] = list(); |
---|
| 226 | |
---|
| 227 | // We have to filter the 0 elements of C |
---|
| 228 | for (i = 1;i <= size(Signs);i++) { |
---|
| 229 | if (C[i,1] != 0) { |
---|
| 230 | result[1] = result[1] + list(Signs[i]); |
---|
| 231 | result[2] = result[2] + list(C[i,1]); |
---|
| 232 | } |
---|
| 233 | } |
---|
| 234 | |
---|
| 235 | return (result); |
---|
| 236 | } |
---|
| 237 | example |
---|
| 238 | { |
---|
| 239 | echo = 2; |
---|
| 240 | ring r = 0,(x,y),dp; |
---|
| 241 | ideal i = (x-2)*(x+3)*x,y*(y-1); |
---|
| 242 | ideal P = x,y; |
---|
| 243 | list l = signcnd(P,i); |
---|
| 244 | echo = 0; |
---|
| 245 | |
---|
| 246 | print("The output of signcnd is a list of two lists. Both lists have the |
---|
| 247 | same"); |
---|
| 248 | print("length. That length is the number of sign conditions realized by the"); |
---|
| 249 | print ("polynomials of P on the set V(i). In this example, that number |
---|
| 250 | is"); |
---|
| 251 | print("print(size(l[1]));"); |
---|
| 252 | print(size(l[1])); |
---|
| 253 | print("Each element of the first list indicates a sign condition of the"); |
---|
| 254 | print("polynomials of P. For example,"); |
---|
| 255 | print("print(l[1][2]);"); |
---|
| 256 | print(l[1][2]); |
---|
| 257 | print("means P[1] > 0,P[2] = 0"); |
---|
| 258 | print("Each element of the second list indicates how many elements of V(I)"); |
---|
| 259 | print("give rise to the sign condition expressed by the same position on the"); |
---|
| 260 | print("first list. For example"); |
---|
| 261 | print("print(l[2][2]);"); |
---|
| 262 | print(l[2][2]); |
---|
| 263 | print("indicates that exactly 1 elemnt of V(I) gives rise to the condition"); |
---|
| 264 | print("P[1] > 0,P[2] = 0."); |
---|
| 265 | print("The procedure psigncnd performs some pretty printing on this output."); |
---|
| 266 | } |
---|
| 267 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 268 | |
---|
| 269 | proc psigncnd(ideal P,list l) |
---|
| 270 | "USAGE: psigncnd(P,I); ideal P, list l |
---|
| 271 | RETURN: list: a formatted version of l |
---|
| 272 | SEE ALSO: signcnd |
---|
| 273 | EXAMPLE: example psigncnd; shows an example" |
---|
| 274 | { |
---|
| 275 | string s; |
---|
| 276 | int n = size(l[1]); |
---|
| 277 | int i; |
---|
| 278 | |
---|
| 279 | for (i = 1;i <= n;i++) { |
---|
| 280 | s = s + string(l[2][i]) + " elements of V(I) satisfy " + psign(P,l[1][i]) |
---|
| 281 | + sprintf("%n",12); |
---|
| 282 | } |
---|
| 283 | return(s); |
---|
| 284 | } |
---|
| 285 | example |
---|
| 286 | { |
---|
| 287 | echo = 2; |
---|
| 288 | ring r = 0,(x,y),dp; |
---|
| 289 | ideal i = (x-2)*(x+3)*x,(y-1)*(y+2)*(y+4); |
---|
| 290 | ideal P = x,y; |
---|
| 291 | list l = signcnd(P,i); |
---|
| 292 | psigncnd(P,l); |
---|
| 293 | } |
---|
| 294 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 295 | |
---|
| 296 | static proc psign(ideal P,list s) |
---|
| 297 | { |
---|
| 298 | int i; |
---|
| 299 | int n = size(P); |
---|
| 300 | string output; |
---|
| 301 | |
---|
| 302 | output = "{P[1]"; |
---|
| 303 | |
---|
| 304 | if (s[1] == -1) { |
---|
| 305 | output = output + " < 0"; |
---|
| 306 | }; |
---|
| 307 | if (s[1] == 0) { |
---|
| 308 | output = output + " = 0"; |
---|
| 309 | }; |
---|
| 310 | if (s[1] == 1) { |
---|
| 311 | output = output + " > 0"; |
---|
| 312 | }; |
---|
| 313 | |
---|
| 314 | for (i = 2;i <= n;i++) { |
---|
| 315 | output = output + ","; |
---|
| 316 | output = output + "P[" + string(i) + "]"; |
---|
| 317 | if (s[i] == -1) { |
---|
| 318 | output = output + " < 0"; |
---|
| 319 | }; |
---|
| 320 | if (s[i] == 0) { |
---|
| 321 | output = output + " = 0"; |
---|
| 322 | }; |
---|
| 323 | if (s[i] == 1) { |
---|
| 324 | output = output + " > 0"; |
---|
| 325 | }; |
---|
| 326 | |
---|
| 327 | } |
---|
| 328 | output = output + "}"; |
---|
| 329 | return (output); |
---|
| 330 | } |
---|
| 331 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 332 | |
---|
| 333 | static proc isIn(list a,list b) //a is a list. b is a list of lists |
---|
| 334 | { |
---|
| 335 | int i,j; |
---|
| 336 | int found; |
---|
| 337 | |
---|
| 338 | found = 0; |
---|
| 339 | i = 1; |
---|
| 340 | while (i <= size(b) && !found) { |
---|
| 341 | j = 1; |
---|
| 342 | found = 1; |
---|
| 343 | if (size(a) != size(b[i])) { |
---|
| 344 | found = 0; |
---|
| 345 | } else { |
---|
| 346 | while(j <= size(a)) { |
---|
| 347 | found = found && a[j] == b[i][j]; |
---|
| 348 | j++; |
---|
| 349 | } |
---|
| 350 | } |
---|
| 351 | i++; |
---|
| 352 | } |
---|
| 353 | |
---|
| 354 | if (found) { |
---|
| 355 | return (i-1); |
---|
| 356 | } else { |
---|
| 357 | return (-1); |
---|
| 358 | } |
---|
| 359 | } |
---|
| 360 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 361 | |
---|
| 362 | static proc expprod(list A,list B) // Computes the product of the list of lists A and the list B. |
---|
| 363 | { |
---|
| 364 | int i,j; |
---|
| 365 | list result; |
---|
| 366 | int la,lb; |
---|
| 367 | |
---|
| 368 | if (size(A) == 0) { |
---|
| 369 | A = list(list()); |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | la = size(A); |
---|
| 373 | lb = size(B); |
---|
| 374 | |
---|
| 375 | result[la*lb] = 0; |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | for (i = 0;i < lb;i++) { |
---|
| 379 | for (j = 0;j < la;j++) { |
---|
| 380 | result[i*la+j+1] = A[j+1] + list(B[i+1]); |
---|
| 381 | } |
---|
| 382 | } |
---|
| 383 | |
---|
| 384 | return (result); |
---|
| 385 | } |
---|
| 386 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 387 | |
---|
| 388 | static proc initlist(int n) // Returns an n-element list of 0s. |
---|
| 389 | { |
---|
| 390 | list l; |
---|
| 391 | int i; |
---|
| 392 | l[n] = 0; |
---|
| 393 | for (i = 1;i < n;i++) { |
---|
| 394 | l[i] = 0; |
---|
| 395 | } |
---|
| 396 | return(l); |
---|
| 397 | } |
---|
| 398 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 399 | |
---|
| 400 | static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate |
---|
| 401 | { |
---|
| 402 | int i; |
---|
| 403 | int n; |
---|
| 404 | poly result; |
---|
| 405 | |
---|
| 406 | n = size(exp); |
---|
| 407 | result = 1; |
---|
| 408 | |
---|
| 409 | for (i = 1;i <= n; i++) { |
---|
| 410 | result = result * (P[i]^exp[i]); |
---|
| 411 | } |
---|
| 412 | return (result); |
---|
| 413 | } |
---|
| 414 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 415 | |
---|
| 416 | static proc incexp(list exp) |
---|
| 417 | { |
---|
| 418 | int k; |
---|
| 419 | |
---|
| 420 | k = 1; |
---|
| 421 | |
---|
| 422 | while (exp[k] == 2) { // We assume exp is not the last exponent (i.e. 2,...,2) |
---|
| 423 | exp[k] = 0; |
---|
| 424 | k++; |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | // exp[k] < 2 |
---|
| 428 | exp[k] = exp[k] + 1; |
---|
| 429 | |
---|
| 430 | return (exp); |
---|
| 431 | } |
---|
| 432 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 433 | |
---|