1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
{$n+,g+,a+}
unit math_rx;
interface
const deg_mx = 20; { degre max du fit }
data_mx = 16383; { nombre max de data }
type
a0 = array[0..deg_mx,0..deg_mx] of extended;
b0 = array[0..deg_mx] of extended;
xy = array[1..data_mx] of single;
function polynome_fits(
px,py : pointer; { pointe les tableaux de x,y de type single dont
le nombre d'element ne doit pas
exceder data_mx }
pf : pointer; { pointe 1 tableau 0..dg of extended pour retour coef }
dg,ndat : byte; { degre ( doit etre <= dg_mx) et nombre de data (doit etre <= ndata_mx )}
var ki2 : single; { pour retour du residu }
_weight : pointer) : { pointeur sur tableau de poids (NIL si non souhait‚)}
boolean;
implementation
const epsilon1 = 1e-25; { arbitrairement petit }
var
a_mat : a0;
arot,
piv,
b_mat : b0;
brot : extended;
function puiss0(s : extended; v : byte) : extended;
var se : extended;
begin
if v=0 then
puiss0:=1
else
if abs(s) < 1e-40 then
puiss0:=0
else
begin
se:=exp(v*ln(abs(s)));
if (se < 0) then if ((v and 1) = 1) then se:=-se;
puiss0:=se;
end;
end;
function polynome_fits( px,py : pointer; pf : pointer ; dg,ndat : byte;
var ki2 : single ; _weight : pointer) : boolean;
procedure make_a_b;
var i,j,k : integer;
begin
if _weight=nil then
begin
for i:=0 to dg do
begin
b_mat[i]:=0;
for k:=1 to ndat do
b_mat[i]:=b_mat[i] + xy(py^)[k] *
puiss0( xy(px^)[k],i);
for j:=0 to dg do
begin
a_mat[i,j]:=0;
for k:=1 to ndat do
a_mat[i,j]:=a_mat[i,j] +
puiss0( xy(px^)[k],i+j)
end
end
end
else
begin
for i:=0 to dg do
begin
b_mat[i]:=0;
for k:=1 to ndat do
begin
b_mat[i]:=b_mat[i] +
xy(py^)[k] * puiss0( xy(px^)[k],i) *
xy(_weight^)[k];
end;
for j:=0 to dg do
begin
a_mat[i,j]:=0;
for k:=1 to ndat do
begin
a_mat[i,j]:=a_mat[i,j] +
puiss0(xy(px^)[k],i+j) *
xy(_weight^)[k];
end;
end
end
end
end;
var toto : boolean;
ktest : integer;
procedure permut (i,ns : integer);
var j,k : integer;
begin
inc(ktest);
brot:=b_mat[i];
for j:=i to dg do arot[j]:=a_mat[i,j];
for j:=i to dg-1 do
begin
for k:=i to dg do a_mat[j,k]:=a_mat[j+1,k];
b_mat[j]:=b_mat[j+1]
end;
for k:=i to dg do
a_mat[dg,k]:=arot[k];
b_mat[dg]:=brot;
toto:=false
end;
function triangle : boolean;
var i,j,k : integer;
begin
triangle:=false;
for i:=0 to dg-1 do
begin
ktest:=0;
repeat
toto:=true;
if ktest > dg then exit;
if a_mat[i,i]=0 then permut(i,dg)
until toto;
for j:=i+1 to dg do piv[j]:=a_mat[i,j]/a_mat[i,i];
for j:=i+1 to dg do
begin
for k:=i to dg do
a_mat[j,k]:=a_mat[j,k]-piv[j]*a_mat[i,k];
b_mat[j]:=b_mat[j]-piv[j]*b_mat[i]
end
end;
triangle:=true
end;
procedure solu;
var u : extended;
i,j : integer;
begin
for i:=dg downto 0 do
begin
u:=b_mat[i];
if i < dg then
for j:=i+1 to dg do u:=u-a_mat[i,j] * b0(pf^)[j];
u:=u/a_mat[i,i];
b0(pf^)[i]:=u;
end
end;
var s,n : extended;
i,j : integer;
begin
polynome_fits:=false;
if dg > deg_mx then exit;
if ndat <= dg then exit;
make_a_b;
if not triangle then exit;
solu;
ki2:=0;
n:=0;
for i:=1 to ndat do
begin
s:=0;
for j:=0 to dg do
s := s + puiss0(xy(px^)[i],j) * b0(pf^)[j];
n:=n + sqr(xy(py^)[i]);
ki2:=ki2 + sqr(xy(py^)[i]-s)
end;
if n > 1e-100 then
ki2:=sqrt(ki2/n)*100
else
ki2:=1e35;
polynome_fits:=true
end;
end. |
polynômes
Envoyé par ford_escort
(ou xmgrace voire meme excel on ne sait jamais)
Répondre avec citation
Partager