Annotation of OpenXM_contrib/PHC/Ada/Homotopy/reduction_of_polynomials.adb, Revision 1.1.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