/*
	This C language subroutine computes the natural
	or base-10 logarithm of the floating point input.
	
	The Run Time Library for the C Language.

	Gordon A. Sterling	(617) 461 - 3076
	DSP Development Tools Engineering

	Created on 4/28/90
	Updated on 5/94 by AS

	#include <math.h>		<Header>
	float logf(float x);		<Prototype>
	float log10f(float x);		<Prototype>
*/

#include "lib_glob.h"
#include "mth_glob.h"

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.PRECISION=MACHINE_PRECISION;

.GLOBAL	    	_logf, _log10f;

_log10f:	JUMP (PC, logs_core) (DB);
		F0=0.43429448190325182765;
		R1=PASS R4, put(R6);

_logf:		JUMP (PC, logs_core) (DB);
		F0=1.0;
		R1=PASS R4, put(R6);

logs_core:	F12=ABS F1, put(R15);
		F15=PASS F0, put(R11);
		R11=LOGB F1, put(R7);
		F1=PASS F1, put(R10);		/*Test input value*/

		IF EQ JUMP (PC, zero_input);	/*Return error and max*/
		IF LT JUMP (PC, set_edom);	/*Return error*/
		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 (PC, adjust_z) (DB);
		F7=F12-F10;			/*znum = f-.5*/
		F8=F7*F10;			/*znum * .5*/
		JUMP (PC, 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' = D*R0*/
		F7=F0*F7, F0=F10-F12;		/*F0=R1=2-D', F7=N*R0*/
		F12=F0*F12;			/*F12=D'=D'*R1*/
		F7=F0*F7, F0=F10-F12;		/*F7=N*R0*R1, F0=R2=2-D'*/
		F12=F0*F12;			/*F12=D'=D'*R2*/
		F7=F0*F7, F0=F10-F12;		/*F7=N*R0*R1*R2, F0=R3=2-D'*/
		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' = D*R0*/
		F7=F0*F7, F0=F10-F12;		/*F0=R1=2-D', F7=N*R0*/
		F12=F0*F12;			/*F12=D'=D'*R1*/
		F7=F0*F7, F0=F10-F12;		/*F7=N*R0*R1, F0=R2=2-D'*/
		F12=F0*F12;			/*F12=D'=D'*R2*/
		F7=F0*F7, F0=F10-F12;		/*F7=N*R0*R1*R2, F0=R3=2-D'*/
		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*/
		F7=F0*F7, F0=F10+F12;		/*F0=XN*C2+R(z), F7=XN*C1*/
		F0=F0+F7;			/*F0 = ln(X)*/

		F0=F0*F15;			/*Compute log base N*/

restore_state:	FETCH_RETURN
		get(R10,1);
		get(R7,2);
		get(R11,3);
		get(R15,4);
		get(R6,5);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

zero_input:     ram_ireg=_errno;
                JUMP (PC, restore_state) (DB);
                rammem(ram_ireg,ram_0)=EDOM;
                F0=0.0;

set_edom:       ram_ireg=_errno;
                JUMP (PC, restore_state) (DB);
                rammem(ram_ireg,ram_0)=EDOM;
                F0=0.0;

.ENDSEG;


.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR logs_data[8] =	 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*/
.ENDSEG;