Annotation of OpenXM_contrib/PHC/Ada/Continuation/ts_vlprs.adb, Revision 1.1.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_Floating_Numbers_io; use Standard_Floating_Numbers_io;
4: with Standard_Random_Numbers; use Standard_Random_Numbers;
5: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
6: with Standard_Integer_Vectors;
7: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
8: with Standard_Floating_Vectors;
9: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
10: with vLpRs_Algorithm; use vLpRs_Algorithm;
11:
12: procedure ts_vlprs is
13:
14: -- DESCRIPTION :
15: -- Test on the extrapolation algorithm on power series.
16:
17: -- GENERATION AND EVALUATION OF A POWER SERIES :
18:
19: function Power_Series_Exponents
20: ( n : natural ) return Standard_Integer_Vectors.Vector is
21:
22: -- DESCRIPTION :
23: -- Asks for a vector of n exponents.
24:
25: res : Standard_Integer_Vectors.Vector(1..n) := (1..n => 0);
26:
27: begin
28: put("Give "); put(n,1); put_line(" integer numbers : ");
29: for i in res'range loop
30: get(res(i));
31: end loop;
32: return res;
33: end Power_Series_Exponents;
34:
35: function Power_Series_Coefficients
36: ( n : natural ) return Standard_Floating_Vectors.Vector is
37:
38: -- DESCRIPTION :
39: -- Returns a vector of range 0..n of randomly generated floating-point
40: -- numbers, which represent the coefficients of a power series.
41:
42: res : Standard_Floating_Vectors.Vector(1..n);
43:
44: begin
45: for i in res'range loop
46: res(i) := Random;
47: end loop;
48: return res;
49: end Power_Series_Coefficients;
50:
51: function Eval_Power_Series
52: ( s : double_float; a : Standard_Floating_Vectors.Vector;
53: m : Standard_Integer_Vectors.Vector ) return double_float is
54:
55: -- DESCRIPTION :
56: -- Returns the sum of a_i*s**m(i).
57:
58: res : double_float := 0.0;
59:
60: begin
61: for i in m'range loop
62: res := res + a(i)*s**m(i);
63: end loop;
64: return res;
65: end Eval_Power_Series;
66:
67: procedure Eval_Log_Power_Series
68: ( s : in double_float; a : in Standard_Floating_Vectors.Vector;
69: m : in Standard_Integer_Vectors.Vector;
70: logs,logx : out double_float ) is
71:
72: -- DESCRIPTION :
73: -- On return are the logarithms of s and the result of the power series.
74:
75: begin
76: logs := LOG10(s);
77: logx := LOG10(Eval_Power_Series(s,a,m));
78: end Eval_Log_Power_Series;
79:
80: procedure Eval_Log_Power_Series
81: ( s,a : in Standard_Floating_Vectors.Vector;
82: m : in Standard_Integer_Vectors.Vector;
83: logs,logx : out Standard_Floating_Vectors.Vector ) is
84:
85: -- DESCRIPTION :
86: -- On return are the logarithms of s and the result of the power series.
87:
88: begin
89: for i in s'range loop
90: logs(i) := LOG10(s(i));
91: logx(i) := LOG10(abs(Eval_Power_Series(s(i),a,m)));
92: end loop;
93: end Eval_Log_Power_Series;
94:
95: -- AUXILIARY :
96:
97: procedure Write ( l,v : in Standard_Floating_Vectors.Vector;
98: m : in integer ) is
99:
100: w,err : double_float;
101:
102: begin
103: for i in v'range loop
104: put(v(i),3,3,3); put(l(i),3,3,3);
105: w := v(i)/l(i); put(" "); put(w);
106: err := abs(w - double_float(m));
107: put(err,3,3,3); new_line;
108: end loop;
109: end Write;
110:
111: -- CALLING the EXTRAPOLATOR :
112:
113: procedure Extrapolate ( n,m,r : in natural;
114: a : in Standard_Floating_Vectors.Vector;
115: exp : in Standard_Integer_Vectors.Vector ) is
116:
117: s,logs,logx : Standard_Floating_Vectors.Vector(0..m-1);
118: srp,dsp : Standard_Floating_Vectors.Vector(1..r-1) := (1..r-1 => 0.0);
119: p : Standard_Floating_Vectors.Vector(0..r-1) := (0..r-1 => 0.0);
120: l,v : Standard_Floating_Vectors.Vector(0..r) := (0..r => 0.0);
121: rt1,rt2 : Matrix(1..r-1,1..r-1);
122:
123: begin
124: put("Give an decreasing sequence of "); put(m,1);
125: put_line(" floating-point numbers :");
126: for i in s'range loop
127: get(s(i));
128: end loop;
129: Eval_Log_Power_Series(s,a,exp,logs,logx);
130: put_line("The logs of the s-values : ");
131: for i in logs'range loop
132: put(logs(i)); new_line;
133: end loop;
134: put_line("The logs of the x-values : ");
135: for i in logx'range loop
136: put(logx(i)); new_line;
137: end loop;
138: rt1(1,1) := 0.0; rt2(1,1) := 0.0;
139: vlprs_full(r,s,logs,logx,srp,dsp,p,l,v,rt1,rt2);
140: put_line("The extrapolated values : ");
141: Write(l,v,exp(exp'first));
142: put_line("The errors of the extrapolated values : ");
143: vlprs_pipe(Standard_Output,r,s,logs,logx,srp,dsp,p,l,v,rt1,rt2);
144: put_line("The extrapolated values : ");
145: Write(l,v,exp(exp'first));
146: end Extrapolate;
147:
148: procedure Test_vlprs is
149:
150: n,m,r : natural;
151:
152: begin
153: put("Give the order of the extrapolator : "); get(r);
154: loop
155: put("Give the number of points : "); get(m);
156: exit when m > r;
157: put("should be larger than "); put(r,1); put_line(". Please retry.");
158: end loop;
159: put("Give the length of the power series : "); get(n);
160: declare
161: exp : Standard_Integer_Vectors.Vector(1..n) := Power_Series_Exponents(n);
162: pcf : Standard_Floating_Vectors.Vector(1..n)
163: := Power_Series_Coefficients(n);
164: begin
165: put("the exponents : "); put(exp); new_line;
166: put_line("the coefficients : ");
167: for i in pcf'range loop
168: put(pcf(i)); new_line;
169: end loop;
170: Extrapolate(n,m,r,pcf,exp);
171: end;
172: end Test_vlprs;
173:
174: begin
175: Test_vlprs;
176: end ts_vlprs;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>