Annotation of OpenXM_contrib/PHC/Ada/Homotopy/reduction_of_polynomials.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 3: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
! 4:
! 5: package body Reduction_of_Polynomials is
! 6:
! 7: -- AUXILIARIES :
! 8:
! 9: function Leading_Term ( p : Poly ) return Term is
! 10:
! 11: -- DESCRIPTION :
! 12: -- Returns the leading term of the polynomial p.
! 13:
! 14: res : Term;
! 15:
! 16: procedure First_Term ( t : in Term; continue : out boolean ) is
! 17: begin
! 18: Standard_Natural_Vectors.Copy
! 19: (Link_to_Vector(t.dg),Link_to_Vector(res.dg));
! 20: res.cf := t.cf;
! 21: continue := false;
! 22: end First_Term;
! 23: procedure Get_First_Term is new Visiting_Iterator(First_Term);
! 24:
! 25: begin
! 26: Get_First_Term(p);
! 27: return res;
! 28: end Leading_Term;
! 29:
! 30: function LEQ ( d1,d2 : Degrees ) return boolean is
! 31:
! 32: -- DESCRIPTION :
! 33: -- Returns true if all degrees of d1 are lower than
! 34: -- or equal to the degrees of d2.
! 35:
! 36: begin
! 37: for i in d1'range loop
! 38: if d1(i) > d2(i)
! 39: then return false;
! 40: end if;
! 41: end loop;
! 42: return true;
! 43: end LEQ;
! 44:
! 45: function Search (p : Poly; pt : Term) return Term is
! 46:
! 47: -- DESCRIPTION :
! 48: -- Returns the first term of p for which the following holds :
! 49: -- LEQ(result,pt).
! 50:
! 51: result : Term;
! 52:
! 53: procedure Search_Term (t : in Term; continue : out boolean) is
! 54: begin
! 55: if LEQ(t.dg,pt.dg)
! 56: then result.cf := t.cf;
! 57: result.dg := new Standard_Natural_Vectors.Vector(t.dg'range);
! 58: for i in t.dg'range loop
! 59: result.dg(i) := t.dg(i);
! 60: end loop;
! 61: continue := false;
! 62: else continue := true;
! 63: end if;
! 64: end Search_Term;
! 65: procedure Search_Terms is new Visiting_Iterator(Search_Term);
! 66:
! 67: begin
! 68: result.cf := Create(0.0);
! 69: Search_Terms(p);
! 70: return result;
! 71: end Search;
! 72:
! 73: procedure Purify ( p : in out Poly; tol : in double_float ) is
! 74:
! 75: -- DESCRIPTION :
! 76: -- All terms of which the coefficient are in AbsVal smaller
! 77: -- than tol are deleted.
! 78:
! 79: procedure Purify_Term ( t : in out Term; continue : out boolean ) is
! 80: begin
! 81: if AbsVal(t.cf) < tol
! 82: then t.cf := Create(0.0);
! 83: end if;
! 84: continue := true;
! 85: end Purify_Term;
! 86: procedure Purify_Terms is new Changing_Iterator(Purify_Term);
! 87:
! 88: begin
! 89: Purify_Terms(p);
! 90: if Number_Of_Unknowns(p) = 0
! 91: then Clear(p);
! 92: p := Null_Poly;
! 93: end if;
! 94: end Purify;
! 95:
! 96: -- TARGET ROUTINES :
! 97:
! 98: function Spoly ( p,q : poly ) return Poly is
! 99:
! 100: S,pp,qq : Poly;
! 101: tfp,tfq,facq,facp : Term;
! 102: abstmp : double_float;
! 103: tol : constant double_float := 10.0**(-13);
! 104:
! 105: begin
! 106: if (p = Null_Poly) or else (Number_Of_Unknowns(p) = 0)
! 107: or else (q = Null_Poly) or else (Number_Of_Unknowns(q) = 0)
! 108: then return Null_Poly;
! 109: end if;
! 110:
! 111: tfp := Leading_Term(p);
! 112: tfq := Leading_Term(q);
! 113:
! 114: facp.dg := new Standard_Natural_Vectors.Vector'(tfp.dg'range => 0);
! 115: facq.dg := new Standard_Natural_Vectors.Vector'(tfq.dg'range => 0);
! 116: for i in tfp.dg'range loop
! 117: if tfp.dg(i) > tfq.dg(i)
! 118: then facq.dg(i) := tfp.dg(i) - tfq.dg(i);
! 119: elsif tfp.dg(i) < tfq.dg(i)
! 120: then facp.dg(i) := tfq.dg(i) - tfp.dg(i);
! 121: end if;
! 122: end loop;
! 123: abstmp := AbsVal(tfp.cf);
! 124: if abstmp > AbsVal(tfq.cf)
! 125: then facp.cf := Create(1.0);
! 126: facq.cf := - tfp.cf / tfq.cf;
! 127: else facp.cf := tfq.cf / tfp.cf;
! 128: facq.cf := Create(-1.0);
! 129: end if;
! 130: pp := facp * p;
! 131: qq := facq * q;
! 132: S := pp + qq;
! 133: Clear(pp); Clear(qq);
! 134: Standard_Natural_Vectors.Clear(Link_to_Vector(facp.dg));
! 135: Standard_Natural_Vectors.Clear(Link_to_Vector(tfp.dg));
! 136: Standard_Natural_Vectors.Clear(Link_to_Vector(facq.dg));
! 137: Standard_Natural_Vectors.Clear(Link_to_Vector(tfq.dg));
! 138: Purify(S,tol);
! 139: return S;
! 140: end Spoly;
! 141:
! 142: function Rpoly ( p,q : Poly ) return Poly is
! 143:
! 144: tol : constant double_float := 10.0**(-13);
! 145:
! 146: begin
! 147: if (p = Null_Poly) or else (Number_Of_Unknowns(p) = 0)
! 148: then return Null_Poly;
! 149: elsif (q = Null_Poly) or else (Number_Of_Unknowns(q) = 0)
! 150: then null;
! 151: else declare
! 152: ltp,tq : Term;
! 153: begin
! 154: ltp := Leading_Term(p);
! 155: tq := Search(q,ltp);
! 156: if tq.cf /= Create(0.0)
! 157: then
! 158: declare
! 159: R,qq : Poly;
! 160: fac : Term;
! 161: begin
! 162: fac.dg
! 163: := new Standard_Natural_Vectors.Vector'(ltp.dg'range => 0);
! 164: for i in ltp.dg'range loop
! 165: if ltp.dg(i) > tq.dg(i)
! 166: then fac.dg(i) := ltp.dg(i) - tq.dg(i);
! 167: end if;
! 168: end loop;
! 169: fac.cf := -ltp.cf / tq.cf;
! 170: qq := fac * q;
! 171: R := p + qq;
! 172: Clear(qq);
! 173: Standard_Natural_Vectors.Clear(Link_to_Vector(fac.dg));
! 174: Purify(R,tol);
! 175: return R;
! 176: end;
! 177: end if;
! 178: Standard_Natural_Vectors.Clear(Link_to_Vector(ltp.dg));
! 179: Standard_Natural_Vectors.Clear(Link_to_Vector(tq.dg));
! 180: end;
! 181: end if;
! 182: declare
! 183: temp : Poly;
! 184: begin
! 185: Copy(p,temp);
! 186: return temp;
! 187: end;
! 188: end Rpoly;
! 189:
! 190: end Reduction_of_Polynomials;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to reduction_of_polynomials.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