DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »polftxt«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦720b7e52e⟧ »calprog« 
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
    └─⟦4334b4c0b⟧ 
        └─⟦720b7e52e⟧ »calprog« 
            └─⟦this⟧ 

TextFile

; 12 09 73 Anders Lindgård
\f


message                  polfit

clear polfit
polfit=set 5
permanent polfit.5
polfit=algol

external


real procedure polfit(i, wi, xi, yi, C, S, r);
value r; integer i, r; real xi, yi, wi; array C,S;

if i<r or i=0 or r<0 then
begin
for r:=r step -1 until 0 do C(r):=S(r):=0;
polfit:=0;
end else
begin integer j, k , l , h, n;
real f, fj, fx, f1, rf, a, b, c, Q, x, y, xn, xk,norm;
array F, F1, X, R(1:i), A, B(0:r), M(1:(r+1)*(r+2)//2);

comment initialisering;
n:=i;
f:= fx:= rf:= 0;
for i:= 1 step 1 until n do
begin fj:= F(i):= sqrt(wi); R(i):= yi * fj;
      X(i):= xi; F1(i):= 0; rf:= rf + R(i) * fj;
      f:= f + fj * fj; fx:= fx + X(i) * fj * fj
end i;
norm:=n/f;
a:= A(0):= fx/f;
b:= 0;
c:= C(0):= rf/f;

comment beregn polynomiet;
for k:= 1 step 1 until r do
begin f1:= f; f:= fx:= rf:= 0;
      for j:= 1 step 1 until n do
      begin fj:= (X(j) - a) * F(j) - b * F1(j);
            R(j):= R(j) - F(j) * c;
            F1(j):= F(j); F(j):= fj;
            rf:= rf + R(j) * fj;
            f:= f + fj * fj; fx:= fx + X(j) * fj * fj
      end j;
      a:= A(k):= fx/f;
      b:= B(k):= f/f1;
      c:= C(k):= rf/f
end beregn;

comment omform til sædvanligt polynomium;
for k:= 1 step 1 until r do
begin C(r-1):= C(r-1) - A(r-k) * C(r);
      for j:= r - 1 step -1 until k do
        C(j-1):= C(j-1) - A(j-k) * C(j) - B(j+1-k) * C(j+1)
end;

comment beregn usikkerheder;
Q:=0;
if n-r<>1 then begin
l:=(r+2)*(r+1)//2;
M(1):=n;
for k:=2 step 1 until l do M(k):=0;
for i:=1 step 1 until n do
  begin
  y:=C(r);
  x:=X(i);
  for h:=r-1 step - 1 until 0 do y:=y*x+C(h);
  f:=y-yi;
  Q:=Q+f*f*wi;
  xk:=1;
  l:=1;
  for k:=2 step 1 until r+1 do
  begin
  xn:=xk:=xk*x;
  for h:=1 step 1 until k do
    begin
    l:=l+1;
    M(l):=M(l)+xn;
    xn:=xn*x;
  end;
    end;
  end;
if -,symin(r+1,M) then system(9,0,<:determinant:>);
Q:=Q/(n-r-1)*norm
end;
for i:=1 step 1 until r+1 do
  begin
  S(i-1):=M(i*(i+1)//2)*Q;
  S(i-1):=if S(i-1)>0 then sqrt(S(i-1)) else
          -sqrt(-S(i-1));
  end;
polfit:=Q;
i:=n;

end polfit;

end
▶EOF◀