Annotation of OpenXM_contrib/PHC/Ada/Continuation/vlprs_tables.adb, Revision 1.1.1.1
1.1 maekawa 1: package body vLpRs_Tables is
2:
3: -- I. The v-table :
4:
5: procedure v_pipe ( v : in out Vector; p : in Vector;
6: vr : in double_float ) is
7:
8: tmp0,tmp1 : double_float;
9:
10: begin
11: tmp0 := v(0);
12: v(0) := vr;
13: for i in 1..v'last loop
14: tmp1 := v(i);
15: v(i) := tmp0 - p(i-1)*v(i-1);
16: tmp0 := tmp1;
17: end loop;
18: end v_pipe;
19:
20: procedure v_pipe ( v,p : in Vector; vr : in double_float;
21: vrp : in out Vector ) is
22: begin
23: vrp(0) := vr;
24: for i in 1..v'last loop
25: vrp(i) := v(i-1) - p(i-1)*vrp(i-1);
26: end loop;
27: end v_pipe;
28:
29: -- II. The L-table :
30:
31: procedure L_pipe ( l : in out Vector; p : in Vector;
32: lr : in double_float ) is
33:
34: tmp0,tmp1 : double_float;
35:
36: begin
37: tmp0 := l(0);
38: l(0) := lr;
39: for i in 1..l'last loop
40: tmp1 := l(i);
41: l(i) := tmp0 - p(i-1)*l(i-1);
42: tmp0 := tmp1;
43: end loop;
44: end L_pipe;
45:
46: procedure L_pipe ( l,p : in Vector; lr : in double_float;
47: lrp : in out Vector ) is
48: begin
49: lrp(0) := lr;
50: for i in 1..l'last loop
51: lrp(i) := l(i-1) - p(i-1)*lrp(i-1);
52: end loop;
53: end L_pipe;
54:
55: -- The combined full version for v- and L-table :
56:
57: procedure vL_full ( s,l,v : in Vector; srp,dsp,p,lrp,vrp : out Vector;
58: rt1,rt2 : in out Matrix ) is
59:
60: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
61: tmpp : Vector(0..s'last-1);
62: prev_lrp,new_lrp,prev_vrp,new_vrp : Vector(s'range);
63: k_last : integer;
64:
65: begin
66: srm1(srm1'first) := s(0);
67: for i in srm1'first+1..srm1'last loop
68: srm1(i) := s(0)*srm1(i-1);
69: end loop;
70: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
71: for j in rt1'range(2) loop
72: rt1(1,j) := tmpdsp(j);
73: end loop;
74: tmpp(0) := 1.0;
75: prev_lrp(0) := l(0); new_lrp(0) := l(1);
76: prev_vrp(0) := v(0); new_vrp(0) := v(1);
77: new_vrp(1) := prev_vrp(0) - new_vrp(0);
78: new_lrp(1) := prev_lrp(0) - new_lrp(0);
79: prev_lrp(0..1) := new_lrp(0..1);
80: prev_vrp(0..1) := new_vrp(0..1);
81: for i in 2..s'last loop
82: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
83: for j in rt2'range(2) loop
84: rt2(1,j) := tmpdsp(j);
85: end loop;
86: if i < s'last
87: then k_last := i; -- compute one additional row
88: else k_last := i-1;
89: end if;
90: for k in 1..k_last loop
91: tmpp(k) := rt1(k,k)/rt2(k,k);
92: for j in k+1..rt2'last(2) loop
93: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
94: end loop;
95: end loop;
96: rt1 := rt2;
97: new_vrp(0) := v(i);
98: new_lrp(0) := l(i);
99: for k in 1..i loop
100: new_vrp(k) := prev_vrp(k-1) - tmpp(k-1)*new_vrp(k-1);
101: new_lrp(k) := prev_lrp(k-1) - tmpp(k-1)*new_lrp(k-1);
102: end loop;
103: prev_vrp(0..i) := new_vrp(0..i);
104: prev_lrp(0..i) := new_lrp(0..i);
105: end loop;
106: srp := tmpsrp;
107: dsp := tmpdsp;
108: p := tmpp;
109: lrp := new_lrp;
110: vrp := new_vrp;
111: end vL_full;
112:
113: -- III. The p-table :
114:
115: procedure p_full ( s : in Vector; srp,dsp,p : out Vector;
116: rt1,rt2 : in out Matrix ) is
117: begin
118: RR_full(s,srp,dsp,p,rt1,rt2);
119: end p_full;
120:
121: procedure p_pipe ( rt1,rt2 : in Matrix; p : out Vector ) is
122: begin
123: p(0) := 1.0;
124: for i in p'first+1..p'last loop
125: p(i) := rt1(i,i)/rt2(i,i);
126: end loop;
127: end p_pipe;
128:
129: -- IV. The R-table :
130:
131: procedure R_full ( s : in Vector; srp,dsp,p : out Vector;
132: rt1,rt2 : in out Matrix ) is
133:
134: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
135: tmpp : Vector(0..s'last-1);
136:
137: begin
138: srm1(srm1'first) := s(0);
139: for i in srm1'first+1..srm1'last loop
140: srm1(i) := s(0)*srm1(i-1);
141: end loop;
142: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
143: for j in tmpdsp'range loop
144: rt1(1,j) := tmpdsp(j);
145: end loop;
146: tmpp(0) := 1.0;
147: for i in 2..s'last loop
148: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
149: for j in tmpdsp'range loop
150: rt2(1,j) := tmpdsp(j);
151: end loop;
152: for k in 1..i-1 loop
153: tmpp(k) := rt1(k,k)/rt2(k,k);
154: for j in k+1..rt2'last(2) loop
155: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
156: end loop;
157: end loop;
158: rt1 := rt2;
159: end loop;
160: srp := tmpsrp;
161: dsp := tmpdsp;
162: p := tmpp;
163: end R_full;
164:
165: procedure RR_full ( s : in Vector; srp,dsp,p : out Vector;
166: rt1,rt2 : in out Matrix ) is
167:
168: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
169: tmpp : Vector(0..s'last-1);
170: k_last,j_first : integer;
171:
172: begin
173: srm1(srm1'first) := s(0);
174: for i in srm1'first+1..srm1'last loop
175: srm1(i) := s(0)*srm1(i-1);
176: end loop;
177: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
178: for j in tmpdsp'range loop
179: rt1(1,j) := tmpdsp(j);
180: end loop;
181: tmpp(0) := 1.0;
182: for i in 2..s'last loop
183: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
184: for j in tmpdsp'range loop
185: rt2(1,j) := tmpdsp(j);
186: end loop;
187: if i < s'last
188: then k_last := i; -- compute one additional row
189: else k_last := i-1;
190: end if;
191: for k in 1..k_last loop
192: tmpp(k) := rt1(k,k)/rt2(k,k);
193: if k < rt2'last(2)
194: then j_first := k;
195: else j_first := k+1;
196: end if;
197: for j in j_first..rt2'last(2) loop -- start earlier with columns
198: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
199: end loop;
200: end loop;
201: rt1 := rt2;
202: end loop;
203: srp := tmpsrp;
204: dsp := tmpdsp;
205: p := tmpp;
206: end RR_full;
207:
208: procedure R_pipe ( rt1 : in Matrix; s,p : in Vector; rt2 : in out Matrix ) is
209: begin
210: rt2(1,1) := s(1);
211: for j in 2..s'last loop
212: rt2(1,j) := s(j);
213: for i in 2..j loop
214: rt2(i,j) := rt1(i-1,j) - p(i-1)*rt2(i-1,j);
215: end loop;
216: end loop;
217: end R_pipe;
218:
219: procedure RR_pipe ( rt1 : in Matrix; s,p : in Vector; rt2 : in out Matrix ) is
220:
221: i_last : integer;
222:
223: begin
224: rt2(1,1) := s(1);
225: for j in 2..s'last loop
226: rt2(1,j) := s(j);
227: if j < rt2'last(1)
228: then i_last := j+1; -- compute one additional row
229: else i_last := j;
230: end if;
231: for i in 2..i_last loop
232: rt2(i,j) := rt1(i-1,j) - p(i-1)*rt2(i-1,j);
233: end loop;
234: end loop;
235: end RR_pipe;
236:
237: -- V. The s-table :
238:
239: procedure s_full ( s : in Vector; srp,dsp : out Vector ) is
240:
241: srm1,tmpsrp : Vector(1..s'last-1);
242:
243: begin
244: srm1(srm1'first) := s(0);
245: for i in srm1'first+1..srm1'last loop
246: srm1(i) := s(0)*srm1(i-1);
247: end loop;
248: for i in 1..s'last loop
249: s_pipe(srm1,s(i),tmpsrp,dsp);
250: srm1 := tmpsrp;
251: end loop;
252: srp := tmpsrp;
253: end s_full;
254:
255: procedure s_pipe ( srp : in out Vector; sr : in double_float;
256: dsp : out Vector ) is
257:
258: tmp : double_float := srp(1);
259: srpow : double_float := sr;
260:
261: begin
262: srp(1) := srpow;
263: dsp(1) := srpow - tmp;
264: for i in 2..srp'last loop
265: srpow := srpow*sr;
266: tmp := srp(i);
267: srp(i) := srpow;
268: dsp(i) := srpow - tmp;
269: end loop;
270: end s_pipe;
271:
272: procedure s_pipe ( sr1 : in Vector; sr : in double_float;
273: srp,dsp : out Vector ) is
274:
275: srpow : double_float := sr;
276:
277: begin
278: srp(1) := srpow;
279: dsp(1) := srpow - sr1(1);
280: for i in sr1'first+1..sr1'last loop
281: srpow := srpow*sr;
282: srp(i) := srpow;
283: dsp(i) := srpow - sr1(i);
284: end loop;
285: end s_pipe;
286:
287: end vLpRs_Tables;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to vlprs_tables.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation