DataMuseum.dk

top - metrics - download

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

└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
    └─⟦b2ec5d50f⟧ 
        └─⟦09b4e9619⟧ »thcømat« 
            └─⟦this⟧ 

TextFile

; 12 09 73 Anders Lindgård
; 07 05 80 Ole J. Heilmann
\f


message polfitw

polfitw=algol message.no
external
real procedure polfitw( X, R, F, n, C, S, r);
value n, r; integer n, r; array X, R, F, C, S;

if n<=0 or r<0 then
begin
  for r:=r step -1 until 0 do
  begin
    C(r):=0; S(r):=maxreal
  end;
  polfitw:=-1
end
else

begin
  integer j, k, l;
  real f, fj, fx, f1, rf, a, b, c, Q, x, y, xn, xk;
  array  F1(1:n), A, B, D(0:r);
  comment initialisering;
  if r>=n then 
  begin
    for j:=r step -1 until n do 
    begin
      C(j):=0; S(j):=maxreal 
    end;
    r:=n-1 
  end;
  f:=fx:= rf:= 0;
  for j:= 1 step 1 until n do 
  begin
    fj:=F(j); F1(j):= 0; rf:=R(j)*fj+rf;
    fj:=fj*fj; f:=fj+f; fx:=fx+X(j)*fj 
  end;
  f1:=1; Q:=f*smallreal**2;
  comment beregn polynomiet;
  for k:= 0 step 1 until r do 
  begin
    a:=A(k):=fx/f; b:=B(k):=f/f1;
    c:=C(k):=rf/f; S(k):=1/f;
    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;
    if f<=Q then 
    begin
      for j:=k+1 step 1 until r do 
      begin
        C(j):=0; S(j):=maxreal 
      end;
      r:=k 
    end 
  end;
  comment omform til sædvanligt polynomium og beregn usikkerheder;
  if n-r>1 then 
  begin
    Q:=0;
    for j:=1 step 1 until n do Q:=R(j)**2+Q;
    Q:=polfitw:=Q/(n-r-1) 
  end
  else 
  begin
    Q:=1; polfitw:=0 
  end;
  for j:=0 step 1 until r do D(j):=0;
  for k:=0 step 1 until r do 
  begin
    c:=D(k); D(k):=a:=1; b:=0; fj:=C(k); rf:=S(k);
    for j:=k+1 step 1 until r do 
    begin
      f:=-A(j-1)*a-B(j-1)*b+c;
      c:=D(j); b:=a; D(j):=a:=f;
      fj:=C(j)*a+fj; rf:=S(j)*a*a+rf 
    end;
    C(k):=fj; S(k):=sqrt(rf*Q) 
  end
end polfitw
; end
▶EOF◀