Annotation of OpenXM_contrib/PHC/Ada/Homotopy/reduction_of_polynomial_systems.adb, Revision 1.1
1.1 ! maekawa 1: --with text_io,integer_io; use text_io,integer_io;
! 2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 4: with Standard_Natural_Vectors;
! 5: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
! 6: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 8: with Reduction_of_Polynomials; use Reduction_of_Polynomials;
! 9:
! 10: package body Reduction_of_Polynomial_Systems is
! 11:
! 12: -- AUXALIARY DATA FOR LINEAR REDUCTION :
! 13:
! 14: mach_eps : constant double_float := 10.0**(-13);
! 15:
! 16: type Degrees_Array is array(positive range <>) of Degrees;
! 17: type Terms_Array is array(positive range <>) of Term;
! 18: type Boolean_Array is array(positive range <>) of boolean;
! 19:
! 20: -- AUXILIARY PROCEDURES FOR LINEAR REDUCTION :
! 21:
! 22: procedure Pop_First_Term ( p : in out Poly; t : in out Term ) is
! 23:
! 24: -- DESCRIPTION :
! 25: -- The term on return is the leading term of p.
! 26: -- This term is removed from p.
! 27:
! 28: procedure First_Term ( tt : in Term; continue : out boolean ) is
! 29: begin
! 30: Copy(tt,t);
! 31: continue := false;
! 32: end First_Term;
! 33: procedure Get_First_Term is new Visiting_Iterator(First_Term);
! 34:
! 35: begin
! 36: t.cf := Create(0.0);
! 37: Get_First_Term(p);
! 38: if t.cf /= Create(0.0)
! 39: then Sub(p,t);
! 40: end if;
! 41: end Pop_First_Term;
! 42:
! 43: procedure Leading_Terms ( p : in out Poly_Sys; ta : in out Terms_Array ) is
! 44:
! 45: -- DESCRIPTION :
! 46: -- Puts the leading terms of the polynomials in p in the array ta.
! 47: -- The leading terms are removed afterwards.
! 48:
! 49: begin
! 50: for i in p'range loop
! 51: Clear(ta(i));
! 52: Pop_First_Term(p(i),ta(i));
! 53: end loop;
! 54: end Leading_Terms;
! 55:
! 56: procedure Find_Max ( ta : in Terms_Array; index : in out Boolean_Array;
! 57: stop : in out boolean ) is
! 58:
! 59: res : Degrees := new Standard_Natural_Vectors.Vector'(ta(1).dg'range => 0);
! 60:
! 61: begin
! 62: stop := true;
! 63: for i in ta'range loop
! 64: if ta(i).cf /= Create(0.0)
! 65: then if ta(i).dg > res
! 66: then res.all := Standard_Natural_Vectors.Vector(ta(i).dg.all);
! 67: index(1..(i-1)) := (1..(i-1) => false);
! 68: index(i) := true;
! 69: stop := false;
! 70: elsif Equal(ta(i).dg,res)
! 71: then index(i) := true;
! 72: stop := false;
! 73: end if;
! 74: end if;
! 75: end loop;
! 76: Standard_Natural_Vectors.Clear
! 77: (Standard_Natural_Vectors.Link_to_Vector(res));
! 78: end Find_Max;
! 79:
! 80: procedure Update ( p : in out Poly_Sys; n : in natural;
! 81: ta : in out Terms_Array; da : in out Degrees_Array;
! 82: nda,cnt : in out natural; mat : in out Matrix;
! 83: stop : in out boolean ) is
! 84:
! 85: index : natural;
! 86: is_max : Boolean_Array(1..n) := (1..n => false);
! 87:
! 88: begin
! 89: Find_Max(ta,is_max,stop);
! 90: if not stop
! 91: then
! 92: -- Get the next Degrees for in the Degrees_Array
! 93: for i in is_max'range loop
! 94: if is_max(i)
! 95: then index := i; exit;
! 96: end if;
! 97: end loop;
! 98: nda := nda + 1;
! 99: Standard_Natural_Vectors.Copy
! 100: (Standard_Natural_Vectors.Link_to_Vector(ta(index).dg),
! 101: Standard_Natural_Vectors.Link_to_Vector(da(nda)));
! 102: -- Fill in the matrix and update the window :
! 103: for i in is_max'range loop
! 104: if is_max(i)
! 105: then mat(i,nda) := ta(i).cf;
! 106: Pop_First_Term(p(i),ta(i));
! 107: cnt := cnt+1;
! 108: else mat(i,nda) := Create(0.0);
! 109: end if;
! 110: end loop;
! 111: end if;
! 112: end Update;
! 113:
! 114: procedure Coefficient_Matrix
! 115: ( p : in Poly_Sys; mat : in out Matrix;
! 116: da : in out Degrees_Array; nda : in out natural;
! 117: diagonal : in out boolean ) is
! 118:
! 119: -- DESCRIPTION :
! 120: -- Constructs the coefficient matrix of the polynomial system.
! 121: -- Stops when the system is diagonal.
! 122:
! 123: -- REQUIRED :
! 124: -- mat'range(1) = p'range, mat'range(2) = 1..Sum_Number_of_Terms(p),
! 125: -- da'range = mat'range(2).
! 126:
! 127: -- ON ENTRY :
! 128: -- p a polynomial system.
! 129:
! 130: -- ON RETURN :
! 131: -- mat coefficient matrix, up to column nda filled up;
! 132: -- da da(1..nda) collects the different terms in the system;
! 133: -- nda number of different terms;
! 134: -- diagonal true if the leading terms are all different.
! 135:
! 136: work : Poly_Sys(p'range);
! 137: stop : boolean := false;
! 138: Window : Terms_Array(p'range);
! 139: cnt : natural := 0;
! 140:
! 141: begin
! 142: Copy(p,work);
! 143: Leading_Terms(work,Window);
! 144: diagonal := false;
! 145: while not stop and not diagonal loop
! 146: Update(work,p'last,Window,da,nda,cnt,mat,stop);
! 147: if (cnt = p'last) and then (cnt = nda)
! 148: then diagonal := true;
! 149: end if;
! 150: exit when diagonal;
! 151: end loop;
! 152: Clear(work);
! 153: end Coefficient_Matrix;
! 154:
! 155: procedure Coefficient_Matrix
! 156: ( p : in Poly_Sys; mat : in out Matrix;
! 157: da : in out Degrees_Array; nda : in out natural ) is
! 158:
! 159: -- DESCRIPTION :
! 160: -- Constructs the coefficient matrix of the polynomial system.
! 161:
! 162: -- REQUIRED :
! 163: -- mat'range(1) = p'range, mat'range(2) = 1..Sum_Number_of_Terms(p),
! 164: -- da'range = mat'range(2).
! 165:
! 166: -- ON ENTRY :
! 167: -- p a polynomial system.
! 168:
! 169: -- ON RETURN :
! 170: -- mat coefficient matrix, up to column nda filled up;
! 171: -- da da(1..nda) collects the different terms in the system;
! 172: -- nda number of different terms.
! 173:
! 174: work : Poly_Sys(p'range);
! 175: stop : boolean := false;
! 176: Window : Terms_Array(p'range);
! 177: cnt : natural := 0;
! 178:
! 179: begin
! 180: Copy(p,work);
! 181: Leading_Terms(work,Window);
! 182: while not stop loop
! 183: Update(work,p'last,Window,da,nda,cnt,mat,stop);
! 184: end loop;
! 185: Clear(work);
! 186: end Coefficient_Matrix;
! 187:
! 188: procedure Make_Polynomial_System
! 189: ( p : in out Poly_Sys; mat : in Matrix;
! 190: da : in Degrees_Array; nda : in natural;
! 191: inconsistent,infinite : out boolean ) is
! 192:
! 193: t : Term;
! 194: n : natural := p'length;
! 195:
! 196: begin
! 197: inconsistent := false;
! 198: infinite := false;
! 199: Clear(p);
! 200: for i in p'range loop
! 201: p(i) := Null_Poly;
! 202: for j in 1..nda loop
! 203: if AbsVal(mat(i,j)) > mach_eps
! 204: then t.dg := new Standard_Natural_Vectors.Vector'(da(j).all);
! 205: t.cf := mat(i,j);
! 206: Add(p(i),t);
! 207: Clear(t);
! 208: end if;
! 209: end loop;
! 210: if p(i) = Null_Poly
! 211: then infinite := true;
! 212: elsif Degree(p(i)) = 0
! 213: then inconsistent := true;
! 214: end if;
! 215: end loop;
! 216: end Make_Polynomial_System;
! 217:
! 218: function Sum_Number_of_Terms ( p : Poly_Sys ) return natural is
! 219:
! 220: -- DESCRIPTION :
! 221: -- Returns the sum of the number of terms of every polynomial in p.
! 222:
! 223: sum : natural := 0;
! 224:
! 225: begin
! 226: for i in p'range loop
! 227: sum := sum + Number_of_Terms(p(i));
! 228: end loop;
! 229: return sum;
! 230: end Sum_Number_of_Terms;
! 231:
! 232: -- TARGET ROUTINES FOR LINEAR REDUCTION :
! 233:
! 234: procedure Reduce ( p : in out Poly_Sys;
! 235: diagonal,inconsistent,infinite : in out boolean ) is
! 236:
! 237: n : natural := p'length;
! 238: cp : Poly_Sys(p'range);
! 239: max_terms : constant natural := Sum_Number_of_Terms(p);
! 240: columns : Degrees_Array(1..max_terms);
! 241: numb_columns : natural := 0;
! 242: mat : Matrix(p'range,1..max_terms);
! 243:
! 244: begin
! 245: Coefficient_Matrix(p,mat,columns,numb_columns,diagonal);
! 246: if diagonal
! 247: then inconsistent := false;
! 248: infinite := false;
! 249: else declare
! 250: coeffmat : Matrix(p'range,1..numb_columns);
! 251: begin
! 252: for i in coeffmat'range(1) loop
! 253: for j in coeffmat'range(2) loop
! 254: coeffmat(i,j) := mat(i,j);
! 255: end loop;
! 256: end loop;
! 257: Triangulate(coeffmat,n,numb_Columns);
! 258: Make_Polynomial_System(p,coeffmat,columns,numb_columns,
! 259: inconsistent,infinite);
! 260: for i in 1..numb_Columns loop
! 261: Standard_Natural_Vectors.Clear
! 262: (Standard_Natural_Vectors.Link_to_Vector(Columns(i)));
! 263: end loop;
! 264: end;
! 265: end if;
! 266: end Reduce;
! 267:
! 268: procedure Sparse_Reduce ( p : in out Poly_Sys;
! 269: inconsistent,infinite : in out boolean ) is
! 270:
! 271: n : natural := p'length;
! 272: max_terms : constant natural := Sum_Number_of_Terms(p);
! 273: columns : Degrees_Array(1..max_terms);
! 274: numb_columns : natural := 0;
! 275: mat : Matrix(1..n,1..max_terms);
! 276:
! 277: begin
! 278: Coefficient_Matrix(p,mat,columns,numb_columns);
! 279: declare
! 280: coeffmat : Matrix(p'range,1..numb_columns);
! 281: begin
! 282: for i in coeffmat'range(1) loop
! 283: for j in coeffmat'range(2) loop
! 284: coeffmat(i,j) := mat(i,j);
! 285: end loop;
! 286: end loop;
! 287: Diagonalize(coeffmat,n,numb_Columns);
! 288: Make_Polynomial_System(p,coeffmat,columns,numb_columns,
! 289: inconsistent,infinite);
! 290: for i in 1..numb_Columns loop
! 291: Standard_Natural_Vectors.Clear
! 292: (Standard_Natural_Vectors.Link_to_Vector(columns(i)));
! 293: end loop;
! 294: end;
! 295: end Sparse_Reduce;
! 296:
! 297: -- NONLINEAR REDUCTION :
! 298:
! 299: function Total_Degree ( p : Poly_Sys ) return natural is
! 300:
! 301: d : natural := 1;
! 302: tmp : integer;
! 303:
! 304: begin
! 305: for i in p'range loop
! 306: tmp := Degree(p(i));
! 307: if tmp >= 0
! 308: then d := d * tmp;
! 309: end if;
! 310: end loop;
! 311: return d;
! 312: end Total_Degree;
! 313:
! 314: function LEQ ( d1,d2 : Degrees ) return boolean is
! 315:
! 316: -- DESCRIPTION :
! 317: -- Returns true if all degrees of d1 are lower than
! 318: -- or equal to the degrees of d2
! 319:
! 320: begin
! 321: for i in d1'range loop
! 322: if d1(i) > d2(i)
! 323: then return false;
! 324: end if;
! 325: end loop;
! 326: return true;
! 327: end LEQ;
! 328:
! 329: function Leading_Term ( p : Poly ) return Term is
! 330:
! 331: -- DESCRIPTION :
! 332: -- Returns the leading term of the polynomial p.
! 333:
! 334: tf : Term;
! 335:
! 336: procedure First_Term (t : in Term; continue : out boolean) is
! 337: begin
! 338: Copy(t,tf);
! 339: continue := false;
! 340: end First_Term;
! 341: procedure Get_First_Term is new Visiting_Iterator (First_Term);
! 342: begin
! 343: Get_First_Term(p);
! 344: return tf;
! 345: end Leading_Term;
! 346:
! 347: function Can_Be_Eliminated ( p : Poly_Sys; j : natural ) return boolean is
! 348:
! 349: -- DESCRIPTION :
! 350: -- returns true if the degree of the j-th unknown in each equation
! 351: -- is zero.
! 352:
! 353: begin
! 354: for i in p'range loop
! 355: if Degree(p(i),j) > 0
! 356: then return false;
! 357: end if;
! 358: end loop;
! 359: return true;
! 360: end Can_Be_Eliminated;
! 361:
! 362: procedure Shift_Null_Polynomial ( p : in out Poly_Sys ) is
! 363:
! 364: -- DESCRIPTION :
! 365: -- The null polynomial in the system p will be shifted down
! 366: -- towards the end.
! 367:
! 368: begin
! 369: for i in p'range loop
! 370: if p(i) = Null_Poly
! 371: then for j in i..(p'last-1) loop
! 372: Copy(p(j+1),p(j));
! 373: Clear(p(j+1));
! 374: end loop;
! 375: end if;
! 376: end loop;
! 377: end Shift_Null_Polynomial;
! 378:
! 379: procedure Eliminate ( p : in out Poly; j : in natural ) is
! 380:
! 381: -- DESCRIPTION :
! 382: -- The j-th unknown will be eliminated out of the polynomial p
! 383:
! 384: n : natural := Number_Of_Unknowns(p);
! 385:
! 386: procedure Eliminate_Term (t : in out Term; continue : out boolean) is
! 387:
! 388: d : Degrees := new Standard_Natural_Vectors.Vector(1..(n-1));
! 389:
! 390: begin
! 391: for i in 1..(j-1) loop
! 392: d(i) := t.dg(i);
! 393: end loop;
! 394: for i in j..(n-1) loop
! 395: d(i) := t.dg(i+1);
! 396: end loop;
! 397: Clear(t);
! 398: t.dg := d;
! 399: continue := true;
! 400: end Eliminate_Term;
! 401: procedure Eliminate_Terms is new Changing_Iterator(Eliminate_Term);
! 402:
! 403: begin
! 404: Eliminate_Terms(p);
! 405: end Eliminate;
! 406:
! 407: procedure Eliminate ( p : in out Poly_Sys; j : in natural ) is
! 408:
! 409: -- DESCRIPTION :
! 410: -- The j-th unknown will be eliminated out of each equation.
! 411:
! 412: begin
! 413: for i in p'range loop
! 414: Eliminate(p(i),j);
! 415: end loop;
! 416: end Eliminate;
! 417:
! 418: procedure Replace ( p : in out Poly_Sys; pp : in Poly; i : in natural ) is
! 419:
! 420: -- DESCRIPTION :
! 421: -- This procedure replaces the i-th polynomial in the system p
! 422: -- by the polynomial pp. If pp is a null polynomial then the procedure
! 423: -- tries to eliminate an unknown, in order to have as much equations
! 424: -- as there are unknowns.
! 425:
! 426: tmp : natural;
! 427:
! 428: begin
! 429: if (pp = Null_Poly) or else (Number_Of_Unknowns(pp) = 0)
! 430: then -- try to eliminate an unknown
! 431: tmp := Number_Of_Unknowns(p(1));
! 432: Clear(p(i)); p(i) := Null_Poly;
! 433: for j in reverse 1..Number_Of_Unknowns(p(1)) loop
! 434: if Can_Be_Eliminated(p,j)
! 435: then Eliminate(p,j);
! 436: end if;
! 437: end loop;
! 438: Shift_Null_Polynomial(p);
! 439: else Clear(p(i)); Copy(pp,p(i));
! 440: end if;
! 441: end Replace;
! 442:
! 443: function red ( p,b1,b2 : Poly ) return Poly is
! 444:
! 445: Rpb1 : Poly := Rpoly(p,b1);
! 446:
! 447: begin
! 448: if Number_Of_Unknowns(Rpb1) = 0
! 449: then return Null_Poly;
! 450: else declare
! 451: Rpb2 : Poly := Rpoly(Rpb1,b2);
! 452: begin
! 453: Clear(Rpb1);
! 454: return Rpb2;
! 455: end;
! 456: end if;
! 457: end red;
! 458:
! 459: function Reduce ( p,b1,b2 : Poly ) return Poly is
! 460:
! 461: -- DESCRIPTION :
! 462: -- returns p mod < b1,b2 >
! 463:
! 464: temp : Poly := red(p,b1,b2);
! 465: begin
! 466: if Number_Of_Unknowns(temp) = 0
! 467: then return Null_Poly;
! 468: else Clear(temp);
! 469: return red(p,b2,b1);
! 470: end if;
! 471: end Reduce;
! 472:
! 473: function Simple_Criterium ( p1,p2 : Poly ) return boolean is
! 474:
! 475: -- DESCRIPTION :
! 476: -- returns true if lcm(in(p1),in(p2)) = in(p1), if in(p2) | in(p1).
! 477:
! 478: lt1,lt2 : Term;
! 479: res : boolean;
! 480:
! 481: begin
! 482: lt1 := Leading_Term(p1);
! 483: lt2 := Leading_Term(p2);
! 484: res := LEQ(lt2.dg,lt1.dg);
! 485: Clear(lt1); Clear(lt2);
! 486: return res;
! 487: end Simple_Criterium;
! 488:
! 489: procedure Rpoly_Criterium ( p,b1,b2 : in Poly; cnt : in out natural;
! 490: res : out boolean ) is
! 491:
! 492: -- DESCRIPTION :
! 493: -- Applies the R-polynomial criterium and counts the number of
! 494: -- R-polynomials computed.
! 495:
! 496: Rpb1 : Poly := Rpoly(p,b1);
! 497: Rpb2 : Poly;
! 498:
! 499: begin
! 500: cnt := cnt + 1;
! 501: if Number_Of_Unknowns(Rpb1) = 0
! 502: then res := true;
! 503: else Rpb2 := Rpoly(Rpb1,b2);
! 504: cnt := cnt + 1;
! 505: Clear(Rpb1);
! 506: if Number_of_Unknowns(Rpb2) = 0
! 507: then res := true;
! 508: else Clear(Rpb2);
! 509: Rpb2 := Rpoly(p,b2);
! 510: cnt := cnt + 1;
! 511: if Number_of_Unknowns(Rpb2) = 0
! 512: then res := true;
! 513: else Rpb1 := Rpoly(Rpb2,b1);
! 514: cnt := cnt + 1;
! 515: Clear(Rpb2);
! 516: if Number_of_Unknowns(Rpb1) = 0
! 517: then res := true;
! 518: else res := false;
! 519: Clear(Rpb1);
! 520: end if;
! 521: end if;
! 522: end if;
! 523: end if;
! 524: end Rpoly_Criterium;
! 525:
! 526: function Criterium ( p,q,s : Poly ) return boolean is
! 527:
! 528: -- DESCRIPTION :
! 529: -- returns true if p may be replaced by s.
! 530:
! 531: begin
! 532: if Simple_Criterium(p,q)
! 533: then return true;
! 534: else declare
! 535: temp : Poly := Reduce(p,q,s);
! 536: res : boolean := (Number_Of_Unknowns(temp) = 0);
! 537: begin
! 538: Clear(temp);
! 539: return res;
! 540: end;
! 541: end if;
! 542: end Criterium;
! 543:
! 544: procedure Criterium ( p,q,s : in Poly; cnt : in out natural;
! 545: res : out boolean ) is
! 546:
! 547: -- DESCRIPTION :
! 548: -- returns true if p may be replaced by s.
! 549:
! 550: begin
! 551: if Simple_Criterium(p,q)
! 552: then res := true;
! 553: else Rpoly_Criterium(p,q,s,cnt,res);
! 554: end if;
! 555: end Criterium;
! 556:
! 557: procedure Reduce ( p : in Poly_Sys; res : in out Poly_Sys;
! 558: cnt_eq : in out natural; max_eq : in natural;
! 559: cnt_sp : in out natural; max_sp : in natural;
! 560: cnt_rp : in out natural; max_rp : in natural ) is
! 561:
! 562: S : Poly;
! 563: n : natural := p'last - p'first + 1;
! 564: dS,dpi,dpj : integer;
! 565: ok : boolean;
! 566:
! 567: procedure try ( i,dpi : in natural ) is
! 568:
! 569: -- DESCRIPTION : try to replace p_i by S
! 570:
! 571: p_red : Poly_Sys(1..n);
! 572:
! 573: begin
! 574: if cnt_eq > max_eq then return; end if;
! 575: if cnt_sp > max_sp then return; end if;
! 576: Clear(p_red); Copy(p,p_red);
! 577: Replace(p_red,S,i);
! 578: -- put("replaced polynomial p("); put(i,1); put_line(").");
! 579: if dS = 0
! 580: then return;
! 581: elsif Total_Degree(p_red) < Total_Degree(res)
! 582: then Copy(p_red,res);
! 583: Reduce(p_red,res,cnt_eq,max_eq,cnt_sp,max_sp,cnt_rp,max_rp);
! 584: elsif cnt_eq <= max_eq
! 585: then cnt_eq := cnt_eq + 1;
! 586: Reduce(p_red,res,
! 587: cnt_eq,max_eq,cnt_sp,max_sp,cnt_rp,max_rp);
! 588: end if;
! 589: Clear(p_red);
! 590: end try;
! 591:
! 592: begin
! 593: if cnt_eq > max_eq then return; end if;
! 594: if cnt_sp > max_sp then return; end if;
! 595: if cnt_rp > max_rp then return; end if;
! 596: for i in 1..n loop
! 597: for j in (i+1)..n loop
! 598: if (p(i) /= Null_Poly) and (p(j) /= Null_Poly)
! 599: then
! 600: Clear(S); S := Spoly(p(i),p(j));
! 601: cnt_sp := cnt_sp + 1;
! 602: dS := Degree(S); dpi := Degree(p(i)); dpj := Degree(p(j));
! 603: -- put("S-polynomial of p("); put(i,1); put(") and p(");
! 604: -- put(j,1); put_line(") computed.");
! 605: if dS <= dpi and then dpi > dpj
! 606: and then Criterium(p(i),p(j),S)
! 607: then try(i,dpi);
! 608: elsif dS <= dpj and then dpi < dpj
! 609: and then Criterium(p(j),p(i),S)
! 610: then try(j,dpj);
! 611: else -- dpi = dpj
! 612: if dS <= dpi
! 613: then Criterium(p(i),p(j),S,cnt_rp,ok);
! 614: if ok then try(i,dpi); end if;
! 615: end if;
! 616: if dS <= dpj
! 617: then Criterium(p(j),p(i),S,cnt_rp,ok);
! 618: if ok then try(j,dpj); end if;
! 619: end if;
! 620: end if;
! 621: Clear(S);
! 622: end if;
! 623: exit when (dS = 0);
! 624: end loop;
! 625: end loop;
! 626: end Reduce;
! 627:
! 628: procedure Sparse_Reduce ( p : in Poly_Sys; res : in out Poly_Sys;
! 629: cnt_eq : in out natural; max_eq : in natural ) is
! 630:
! 631: S : Poly;
! 632: n : natural := p'last - p'first + 1;
! 633: dS,dpi,dpj : integer;
! 634:
! 635: procedure try ( i,dpi : in natural ) is
! 636:
! 637: -- DESCRIPTION : try to replace p_i by S
! 638:
! 639: p_red : Poly_Sys(1..n);
! 640: inconsistent,infinite : boolean := false;
! 641: begin
! 642: if cnt_eq > max_eq then return; end if;
! 643: Clear(p_red); Copy(p,p_red);
! 644: Replace(p_red,S,i);
! 645: if dS /= 0
! 646: then Sparse_Reduce(p_red,inconsistent,infinite);
! 647: end if;
! 648: if dS = 0 or inconsistent
! 649: then cnt_eq := max_eq + 1;
! 650: return;
! 651: elsif Total_Degree(p_red) < Total_Degree(res)
! 652: then Copy(p_red,res);
! 653: Sparse_Reduce(p_red,res,cnt_eq,max_eq);
! 654: else cnt_eq := cnt_eq + 1;
! 655: Sparse_Reduce(p_red,res,cnt_eq,max_eq);
! 656: end if;
! 657: Clear(p_red);
! 658: end try;
! 659:
! 660: begin
! 661: if cnt_eq > max_eq then return; end if;
! 662: for i in 1..n loop
! 663: for j in (i+1)..n loop
! 664: if (p(i) /= Null_Poly) and (p(j) /= Null_Poly)
! 665: then
! 666:
! 667: Clear(S); S := Spoly(p(i),p(j));
! 668:
! 669: dS := Degree(S);
! 670: dpi := Degree(p(i));
! 671: dpj := Degree(p(j));
! 672:
! 673: if dS <= dpi and then dpi > dpj
! 674: and then Criterium(p(i),p(j),S)
! 675: then try(i,dpi);
! 676: elsif dS <= dpj and then dpi < dpj
! 677: and then Criterium(p(j),p(i),S)
! 678: then try(j,dpj);
! 679: else -- dpi = dpj
! 680: if dS <= dpi
! 681: and then Criterium(p(i),p(j),S) then try(i,dpi); end if;
! 682: if dS <= dpj
! 683: and then Criterium(p(j),p(i),S) then try(j,dpj); end if;
! 684: end if;
! 685:
! 686: Clear(S);
! 687:
! 688: end if;
! 689: end loop;
! 690: end loop;
! 691: end Sparse_Reduce;
! 692:
! 693: end Reduction_of_Polynomial_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to reduction_of_polynomial_systems.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy