DataMuseum.dk

top - metrics - download

    Length: 3072 (0xc00)
    Types: TextFile
    Names: »tspercre«

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

TextFile

;rene moss 12-1-1973/6-8-1980

spercre=algol
\f


external procedure spercre(AR,n,AI);
value n; integer n; real array AR,AI;
begin
real u,v,w,k,m,f,fm,fc,xm,ym,xr,yr,xc,yc,dx,dy;
integer i,j,p,ovf;
  RED:
  if n>0 then
  begin
    if AR(n)=0 then
    begin
      AI(n):=0; n:=n-1;
      goto RED
    end;
    ovf:=overflows; overflows:=0;
    xc:=0.0;
    p:=n-1;
    f:=fc:=AR(n)*AR(n);
    dx:=abs(AR(n)/AR(0))**(1/n);
    m:=dx*dx;
    for k:=dx+dx,-k do
    begin
      fm:=f;
      u:=v:=0.0;
      for j:=0 step 1 until p do
      begin
        w:=AR(j)+k*u-m*v;
        v:=u; u:=w
      end;
      f:=(AR(n)+u*k/2-m*v)**2
    end;
    if f<fm then
    begin
      dx:=-dx; fm:=f
    end
    else k:=-k;
    for m:=4*m while fm<fc do
    begin
      xc:=dx; dx:=k; k:=dx+dx;
      fc:=fm; u:=v:=0;
      for j:=0 step 1 until p do
      begin
        w:=AR(j)+k*u-m*v;
        v:=u; u:=w
      end;
      fm:=(AR(n)+u*dx-m*v)**2
    end;
    yc:=abs(xc); dx:=abs(dx);
    for dx:=0.5*dx while yc+dx>yc do
    begin
      for xr:=xc+dx,xc-dx do
      begin
        u:=v:=0;
        k:=xr+xr;
        m:=xr*xr;
        for j:=0 step 1 until p do
        begin
          w:=AR(j)+k*u-m*v;
          v:=u; u:=w
        end;
        f:=(AR(n)+u*xr-m*v)**2;
        if f<fc then
        begin
          xc:=xr; yc:=abs(xc);
          fc:=f;
          goto ITERE
        end
      end;
      ITERE:
    end;
    overflows:=ovf;
    dy:=yc:=0;
    dx:=(fc/AR(0)**2)**(1/2/n);
    goto TEST;
    ITER:
    if overflows>0 then
    begin
      n:=0;
      goto RED
    end;
    fm:=fc+fc;
    for i:=1,2,3,4 do
    begin
      u:=-dy; dy:=dx; dx:=u;
      xr:=xc+dx; yr:=yc+dy;
      u:=v:=0;
      k:=xr+xr;
      m:=xr*xr+yr*yr;
      for j:=0 step 1 until p do
      begin
        w:=AR(j)+k*u-m*v; 
        v:=u; u:=w
      end;
      f:=(AR(n)+u*xr-m*v)**2+u*u*yr*yr;
      if f<fm then
      begin
        xm:=xr; ym:=yr; fm:=f
      end
    end;
    if fm<=fc then
    begin
      dx:=dx*1.5; dy:=dy*1.5;
      xc:=xm; yc:=ym; fc:=fm
    end
    else
    begin
      u:=0.4*dx-0.3*dy; dy:=0.4*dy+0.3*dx; dx:=u
    end;
    TEST:
    u:=abs(xc)+abs(yc);
    if u+abs(dx)+abs(dy)>u and fc<>0 then goto ITER;
    u:=v:=0;
    k:=xc+xc;
    m:=xc*xc;
    for j:=0 step 1 until p do
    begin
      w:=AR(j)+k*u-m*v;
      v:=u;  u:=w
    end;
    if (AR(n)+u*xc-m*v)**2<=fc then
    begin
      u:=0;
      for j:=0 step 1 until p do u:=AR(j):=u*xc+AR(j);
      AR(n):=xc; AI(n):=0
    end
    else
    begin
      u:=v:=0;
      k:=xc+xc;
      m:=xc*xc+yc*yc;
      p:=n-2;
      for j:=0 step 1 until p do
      begin
        w:=AR(j):=AR(j)+k*u-m*v;
        v:=u; u:=w
      end;
      AR(n-1):=AR(n):=xc;
      AI(n-1):=-yc; AI(n):=yc
    end;
    n:=p;
    goto RED
  end
end spercre
; end
▶EOF◀