%% Copyright Massachusetts Institute of Technology 1982, 1989

%% arctan(x) function.
%% Uses arctan approximation #4900 from
%% "Computer Approximations", John F. Hart, et. al.
%% This procedure does the appropriate range reductions and then
%% call atan4900 to do the approximation.
 
%% For the value of the result will be in the range [-(PI/2), PI/2]

ATan = proc (X: real) returns (real);
    %% To use atan4900 the argument must be reduced to the range
    %% [-tan(PI/32), tan(PI/32)]
    %% See Hart, et. al., page 125 for description how this is done.

    rs = sequence[real]
    %% Partition points.
    PP1 =  0.098491403357  %% tan(PI/32)
    PP2 =  0.303346683607  %% tan(3PI/32)
    PP3 =  0.534511135951  %% tan(5PI/32)
    PP4 =  0.820678790829  %% tan(7PI/32)
    PP5 =  1.218503525588  %% tan(9PI/32)
    PP6 =  1.870868411789  %% tan(11PI/32)
    PP7 =  3.296558208938  %% tan(13PI/32)
    PP8 = 10.15317387609   %% tan(15PI/32)

    %% tan**-1 of central points.
    CP1 = 1./0.198912367380  %% 1/tan(2PI/32)  = 1/tan(PI/16)
    CP2 = 1./0.414213562373  %% 1/tan(4PI/32)  = 1/tan(PI/8)
    CP3 = 1./0.668178637919  %% 1/tan(6PI/32)  = 1/tan(3PI/16)
    CP4 = 1./1.000000000000  %% 1/tan(8PI/32)  = 1/tan(PI/4)
    CP5 = 1./1.496605762665  %% 1/tan(10PI/32) = 1/tan(5PI/16)
    CP6 = 1./2.414213562373  %% 1/tan(12PI/32) = 1/tan(3PI/8)
    CP7 = 1./5.027339492126  %% 1/tan(14PI/32) = 1/tan(7PI/16)
    CP8 = 0.0                %% 1/tan(16PI/32) = 1/tan(PI/2)
    CPs = rs$[ CP1, CP2, CP3, CP4, CP5, CP6, CP7, CP8 ];
    %% (tan**-2 + 1) of central points.
    CPsSqP1 = rs$[
		  (CP1*CP1 + 1.),
		  (CP2*CP2 + 1.),
		  (CP3*CP3 + 1.),
		  (CP4*CP4 + 1.),
		  (CP5*CP5 + 1.),
		  (CP6*CP6 + 1.),
		  (CP7*CP7 + 1.),
		  (CP8*CP8 + 1.)
		  ];

    %% arctans of central points.
    PI = 3.141592653589793
    ATs = rs$[
	      PI/16.,         %%  PI/16
	      PI/8.,          %%  PI/8
	      (3.*PI)/16.,    %% 3PI/16
	      PI/4.,          %%  PI/4
	      (5.*PI)/16.,    %% 5PI/16
	      (3.*PI)/8.,     %% 3PI/8
	      (7.*PI)/16.,    %% 7PI/16
	      PI/2.           %%  PI/2
	      ];

    Is_Neg: bool := X < 0.;
    X := real$Abs(X);
    
    %% Search for proper range.   Note that the range points have odd
    %% indexes, while the mid points have even.
    I: int;
    %% Search tree.
    if (X < PP4)
       then if (X < PP2)
	       then if (X < PP1)
		       then I := 0;
		       else I := 1;
		       end;
	       else if (X < PP3)
		       then I := 2;
		       else I := 3;
		       end;
	       end;
       else if (X < PP6)
	       then if (X < PP5)
		       then I := 4;
		       else I := 5;
		       end;
	       else if (X < PP7)
		       then I := 6;
		       else if (X < PP8)
			       then I := 7;
			       else I := 8;
			       end;
		       end;
	       end;
       end;
    
    %% t = CPs[I] - (CPs[I]*CPs[I] + 1) / (CPs[I] + X)
    %% t in range [-PI/2, PI/2]
    T: real;
    ATX: real;		%% arctan(CP[I])
    if (I = 0)
       then T   := X;
	    ATX := 0.;
       else T   := CPs[I] - (CPsSqP1[I] / (CPs[I] + X));
	    ATX := ATs[I];
       end;

    ATX := ATX + ATan4900(T);
    if (Is_Neg)
       then ATX := -ATX;
       end;
    return (ATX);
    end ATan;