{ Subroutine to compute the ArcTangent values of a floating point input.

atan-       H=ATAN(Y)
atan2-      H=ATAN(Y/X)         
                            where H is in radians

        Calling Registers
    F0 = Input Value Y=[6E-20, 6E20] 
    F1 = Input Value X=[6E-20, 6E20] (atan2 only)
    l_reg=0

        Result Registers
    F0 = ArcTangent of input
            =[-pi/2,pi/2]   for atan
            =[-pi,pi]       for atan2

        Altered Registers
    F0, F1, F2, F4, F7, F8, F11, F12, F15
    i_reg
    ms_reg

        Computation Time
    atan    61 Cycles maximum
    atan2   82 cycles maximum

    Version 0.01   3/20/91      Ronnin Yee
    (based on Gordon Sterlings c lib routine)
}

{The divide Macro used in the arctan routine}
{-----------------------------------------------------------------------------
DIVIDE - Divide Macro
        
Register Usage:
    q       = f0-f15        Quotient
    n       = f4-f7         Numerator
    d       = f12-f15       Denominator
    two     = f8-f11        must have 2.0 pre-stored
    tmp     = f0-f3

Indirectly affected registers:
    ASTAT,STKY

Looping:    none

Special Cases:
    q may be any register, all others must be distinct.
    
Cycles: 8

Created:    3/19/91
-----------------------------------------------------------------------------}

#define DIVIDE(q,n,d,two,tmp)                                           \
        n=RECIPS d, tmp=n;          {Get 8 bit seed R0=1/D}             \
        d=n*d;                      {D(prime) = D*R0}                   \
        tmp=tmp*n, n=two-d;         {N=2-D(prime), TMP=N*R0}            \
        d=n*d;                      {D=D(prime)=D(prime)*R1}            \
        tmp=tmp*n, n=two-d;         {TMP=N*R0*R1, N=R2=2-D(prime)}      \
        d=n*d;                      {D=D(prime)=D(prime)*R2}            \
        tmp=tmp*n, n=two-d;         {TMP=N*R0*R1*R2, N=R3=2-D(prime)}   \
        q=tmp*n


#include "asm_glob.h"

.SEGMENT/PM     Assembly_Library_Code_Space;

.PRECISION=MACHINE_PRECISION;

.GLOBAL     atan, atan2;

atan2:          i_reg=atan_data;
                F11= 2.0;
                F2= 0.0;
                F1=PASS F1;
                IF EQ JUMP denom_zero;      {if Denom. = 0, goto special case}
                IF LT F2=mem(11,i_reg);     {if Denom. < 0, F2=pi (use at end)}

overflow_tst:   R4=LOGB F0, F7=F0;          {Detect div overflow for atan2}
                R1=LOGB F1, F15=F1;
                R1=R4-R1;                   {Roughly exp. of quotient}
                R4=124;                     {Max exponent - 3}
                COMP(R1,R4);
                IF GE JUMP overflow;        {Over upper range? Goto overflow}
                R4=-R4;
                COMP(R1,R4);
                IF LE JUMP underflow;       {Over lower range? Goto underflow}

do_division:    DIVIDE(F0,F7,F15,F11,F1);
                JUMP re_entry (DB);

atan:           R10= 0;                     {Flags multiple of pi to add at end}
                F15=ABS F0;             

                i_reg=atan_data;            {This init is redundant for atan2f}
                F11= 2.0;                       {Needed for divide}
                F2= 0.0;                        {Result is not in Quad 2 or 3 }

re_entry:       F7 = 1.0;
                COMP(F15,F7), F4=mem(0,i_reg);   {F4=2-sqrt(3)}
                IF LE JUMP tst_f;           {If input<=1, do arctan(input)}
                                                {else do arctan(1/input)+const}
                DIVIDE(F15,F7,F15,F11,F1);  {do x=1/x}

                R10 = 2;                    {signal to add const at end}

tst_f:          COMP(F15,F4);               {Note F4 prev. loaded from memory}
                IF LT JUMP tst_for_eps; 
                R10=R10+1,  F4=mem(1,i_reg);    {F4=sqrt(3)}    
                F12=F4*F15;             
                F7=F12-F7;                  {F7=F12-1.0}            
                F15=F4+F15;             
                DIVIDE(F15,F7,F15,F11,F1);  { = (sqrt(3)*x-1)/(x+sqrt(3))}

tst_for_eps:    F7=ABS F15, F4=mem(2,i_reg);    {F4=eps (i.e. small)}   
                COMP(F7,F4);
                IF LE JUMP tst_N;           {if x<=eps, then h=x}
                F1=F15*F15, F4=mem(3,i_reg);{else . . .}        
                F7=F1*F4,   F4=mem(4,i_reg);            
                F7=F7+F4,   F4=mem(5,i_reg);            
                F7=F7*F1;               
                F12=F1+F4,  F4=mem(6,i_reg);            
                F12=F12*F1;
                F12=F12+F4;
                DIVIDE(F7,F7,F12,F11,F1);   {e=((p1*x^2 +p0)x^2)/(x^2 +q1)x^2+q0}
                F7=F7*F15;              
                F15=F7+F15;                 {h=e*x+x}

tst_N:          R1=R10-1, R7=mem(i_reg,7);  {if R10 > 1, h=-h; dummy read}
                ms_reg=R10;
                IF GT F15=-F15; 
                F4 = mem(ms_reg,i_reg);     {index correct angle addend to h }
                F15=F15+F4;                 {h=h+ a*pi      }

tst_sign_y:     F2=PASS F2;                 {if (atan2f denom <0) h=pi-h else h=h  }
                IF NE F15=F2-F15;           
tst_sign_x:     RTS (DB), F0=PASS F0;       {if (numer<0) h=-h else h=h}
                IF LT F15=-F15;
                F0=PASS F15;                {return with result in F0!}
                                

underflow:      JUMP tst_sign_y (DB);
                F15=0;
                NOP;


overflow:       JUMP tst_sign_x (DB);
                F15=mem(9,i_reg);           {Load pi/2}
                NOP;

                
denom_zero:     F0=PASS F0;                 {Whooh doggie, if Num !=0, then overflow}
                IF NE JUMP overflow;
    
error:          RTS;                        {Yikes! Its a 0/0!}


.ENDSEG;


.SEGMENT/SPACE  Assembly_Library_Data_Space;

.PRECISION=MEMORY_PRECISION;

.VAR atan_data[11] =    0.26794919243112270647,     {2-sqrt(3)}
                        1.73205080756887729353,     {sqrt(3)}
                        0.000244140625,             {eps}
                        -0.720026848898E+0,         {p1}
                        -0.144008344874E+1,         {p0}
                        0.475222584599E+1,          {q1}
                        0.432025038919E+1,          {q0}
                        0.00000000000000000000,     {0*pi}
                        0.52359877559829887308,     {pi/6}
                        1.57079632679489661923,     {pi/2}
                        1.04719755119659774615,     {pi/3}
                        3.14159265358979323846;     {pi}
.ENDSEG;