{
    Subroutine to compute the logarithm (bases 2,e, and 10)
    of its floating point input.

	    Y=LN(X) or
	    Y=LOG(X) or
	    Y=LOG2(X)

	Calling Registers
    F0 = Input Value
    l_reg=0;

	Result Registers
    F0 = Logarithm of input

	Altered Registers
    F0, F1, F6, F7, F8, F10, F11, F12
    i_reg

	Computation Time
    49 Cycles

    Version 0.04    6/24/91	Gordon A. Sterling
}

#include "asm_glob.h"

.SEGMENT/PM	Assembly_Library_Code_Space;

.PRECISION=MACHINE_PRECISION;

.GLOBAL	    log, log10, logb2;

logb2:		CALL logs_core (DB);		{Enter same routine in two cycles}
		R11=LOGB F0, F1=F0;		{Extract the exponent}
		F12=ABS F1;			{Get absolute value}
		RTS (DB);
		F11=mem(1,i_reg);		{1/Log(2)}
		F0=F11*F0;			{F12 = log2(X)}

log10:		CALL logs_core (DB);		{Enter same routine in two cycles}
		R11=LOGB F0, F1=F0;		{Extract the exponent}
		F12=ABS F1;			{Get absolute value}
		RTS (DB);
		F11=mem(i_reg,1);		{1/log(10)}
		F0=F11*F0;			{F12 = log10(X)}

log:		R11=LOGB F0, F1=F0;
		F12=ABS F1;

logs_core:	i_reg=logs_data;		{Point to data array}
		R11=R11+1;			{Increment exponent}
		R7=-R11, F10=mem(i_reg,1);	{Negate exponent}
		F12=SCALB F12 BY R7;		{F12= .5<=f<1}
		COMP(F12,F10), F10=mem(i_reg,1);{Compare f > C0}
		IF GT JUMP adjust_z (DB);
		F7=F12-F10;			{znum = f-.5}
		F8=F7*F10;			{znum * .5}
		JUMP compute_r (DB);		
		F12=F8+F10;			{zden = znum * .5 + .5}	
		R11=R11-1;			{N = N - 1}

adjust_z:	F7=F7-F10;			{znum = f - .5 - .5}	
		F8=F12*F10;			{f * .5}
		F12=F8+F10;			{zden = f * .5 + .5}

compute_r:	F0=RECIPS F12;			{Get 4 bit seed R0=1/D}
		F12=F0*F12, F10=mem(i_reg,1);	{D(prime) = D*R0}
		F7=F0*F7, F0=F10-F12;		{F0=R1=2-D(prime), F7=N*R0}
		F12=F0*F12;			{F12=D(prime)=D(prime)*R1}
		F7=F0*F7, F0=F10-F12;		{F7=N*R0*R1, F0=R2=2-D(prime)}
		F12=F0*F12;			{F12=D(prime)=D(prime)*R2}
		F7=F0*F7, F0=F10-F12;		{F7=N*R0*R1*R2, F0=R3=2-D(prime)}
		F6=F0*F7;			{F7=N*R0*R1*R2*R3}
		F0=F6*F6, F8=mem(i_reg,1);	{w = z^2}
		F12=F8+F0, F8=mem(i_reg,1);	{B(W) = w + b0}
		F7=F8*F0, F8=mem(i_reg,1);	{ w*a1}
		F7=F7+F8, F8=F0;		{A(W) = w * a1 + a0}

		F0=RECIPS F12;			{Get 4 bit seed R0=1/D}
		F12=F0*F12;			{D(prime) = D*R0}
		F7=F0*F7, F0=F10-F12;		{F0=R1=2-D(prime), F7=N*R0}
		F12=F0*F12;			{F12=D(prime)=D(prime)*R1}
		F7=F0*F7, F0=F10-F12;		{F7=N*R0*R1, F0=R2=2-D(prime)}
		F12=F0*F12;			{F12=D(prime)=D(prime)*R2}
		F7=F0*F7, F0=F10-F12;		{F7=N*R0*R1*R2, F0=R3=2-D(prime)}
		F7=F0*F7;			{F7=N*R0*R1*R2*R3}
		F7=F7*F8;			{Compute r(z^2)=w*A(w)/B(w)}

compute_R:	F7=F6*F7;			{ z*r(z^2)}
		F12=F6+F7;			{R(z) = z + z * r(z^2)}
		F0=FLOAT R11, F7=mem(i_reg,1);	{F0=XN, F7=C2}
		F10=F0*F7, F7=mem(i_reg,1);	{F10=XN*C2, F7=C1}
		RTS (DB);
		F7=F0*F7, F0=F10+F12;		{F0=XN*C2+R(z), F7=XN*C1}
		F0=F0+F7;			{F0 = ln(X)}

.ENDSEG;

.SEGMENT/SPACE	Assembly_Library_Data_Space;

.VAR logs_data[10] =	 0.70710678118654752440,	{C0 = sqrt(.5)}
			 0.5,				{Constant used}
			 2.0,				{Constant used}
			-5.578873750242,		{b0}
			 0.1360095468621E-1,		{a1}
			-0.4649062303464,		{a0}
			-2.121944400546905827679E-4,	{C2}
			 0.693359375,			{C1}
  			 0.43429448190325182765,	{1/log(10)}
			 1.442695041;                   {1/Log(2)}

.ENDSEG;