[pythagresult,a,b]=pythag(a,b);

function [pythagresult,a,b]=pythag(a,b);
pythagresult=[];
persistent p q r s t ; 

;
%***BEGIN PROLOGUE  PYTHAG
%***SUBSIDIARY
%***PURPOSE  Compute the complex square root of a complex number without
%            destructive overflow or underflow.
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (PYTHAG-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     Finds sqrt(A**2+B**2) without overflow or destructive underflow
%
%***SEE ALSO  EISDOC
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   811101  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  PYTHAG
%
if isempty(p), p=0; end;
if isempty(q), q=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
%***FIRST EXECUTABLE STATEMENT  PYTHAG
p = max(abs(a),abs(b));
q = min(abs(a),abs(b));
if( q==0.0e0 )
pythagresult = p;
else;
while( true );
r =(q./p).^2;
t = 4.0e0 + r;
if( t==4.0e0 )
break;
end;
s = r./t;
p = p + 2.0e0.*p.*s;
q = q.*s;
end;
pythagresult = p;
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end %function pythag
%DECK QAGE