/*
	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>
	double log(double x);		<Prototype>
	double log10(double x);		<Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_log, _log10;

_log10:		JUMP (PC, logs_core) (DB);
		pm_mdf = 1;
		R2=MR2F;

_log:		JUMP (PC, logs_core) (DB);
		pm_mdf = 0;
		R2=MR2F;

logs_core:	put(R2);
		R2=MR1F;
		put(R2);
		R2=MR0F;
		put(R2);
		put(R3);
		put(R5);
		put(R6);
		put(R7);
		put(R9);
		put(R10);
		put(R11);
		put(R13);
		put(R14);
		put(R15);

		reads(R4,1);
		reads(R8,2);

		i_reg=logs_data;		/* Point to data array		*/
		R4=PASS R4, l_reg=dm_0;		/* Test for <= 0		*/
		IF EQ JUMP (PC, zero_input);	/* Return error and max		*/
		IF LT JUMP (PC, set_edom);	/* Return error and 0		*/

estimate_N:	R2=FEXT R4 BY 20:11;		/* Get exponent-1 from input	*/
		R12=D_EXP_BIAS;			/* Get double exponent bias	*/
		R2=R2-R12, R11=R8;		/* Unbiased exponent-1, copy LSW*/
		R2=R2+1;			/* Compute exponent		*/
		dm_mdf=R2;			/* Hold N for later		*/

estimate_f:	R13=0x800fffff;			/* Mantissa mask		*/
		R13=R4 AND R13;			/* Mantissa with no exponent	*/
		R12=R12-1;			/* Used D_EXP_BIAS - 1 as exp	*/
		R13=R13 OR ASHIFT R12 BY 20;	/* Place exponent		*/

compare_to_C0:	R9=0x3fe00000;			/* Load MSW of 0.5 for later	*/
		R10=PASS R13, R12=mem(i_reg,1);	/* Fetch MSW of C0		*/
		COMP(R13,R12), R12=mem(i_reg,1);/* Test for f > C0, fetch LSW	*/
		IF GT JUMP big_f;		/* f is larger than C0		*/

small_f:	R0=dm_mdf;			/* Need to compute N-1		*/
		R0=R0-1;
		dm_mdf=R0;
		CALL (PC, ___lib_dsub) (DB);	/* Subtract in double precision	*/
		R14=PASS R10, R15=R11;		/* Copy f			*/
		R12=PASS R9, R13=dm_0;		/* Load 0.5			*/

		R10=PASS R0, R11=R1;		/* Hold znum for now		*/
		CALL (PC, ___lib_dmult) (DB);	/* Add input double precision	*/
		R12=PASS R9, R13=dm_0;		/* Load 0.5			*/
		R14=PASS R0, R15=R1;		/* Compute znum * 0.5		*/
		
		CALL (PC, ___lib_dadd) (DB);	/* Add input double precision	*/
		R14=PASS R0, R15=R1;		/* Compute znum * 0.5 + 0.5	*/
		R12=PASS R9, R13=dm_0;		/* Load 0.5			*/
		
		JUMP (PC, evaluate_R);

big_f:		CALL (PC, ___lib_dsub) (DB);	/* Subtract in double precision	*/
		R14=PASS R10, R15=R11;		/* Copy f			*/
		R12=PASS R9, R13=dm_0;		/* Copy 0.5			*/

		CALL (PC, ___lib_dsub) (DB);	/* Compute (f-0.5)-0.5		*/
		R14=PASS R0, R15=R1;		/* Copy f-0.5			*/
		R12=PASS R9, R13=dm_0;		/* Copy 0.5			*/

		R12=PASS R9, R13=dm_0;		/* Copy 0.5			*/
		CALL (PC, ___lib_dmult) (DB);	/* Multiply in double precision	*/
		R10=PASS R0, R14=R10;		/* Hold znum, get f		*/
		R11=PASS R1, R15=R11;		/* Same for LSW			*/

		CALL (PC, ___lib_dadd) (DB);	/* Compute f * .5 + .5		*/
		R14=PASS R0, R15=R1;		/* Copy f * .5			*/
		R12=PASS R9, R13=dm_0;		/* Copy 0.5			*/
		
evaluate_R:	put(R1);			/* Store zden LSW		*/
		put(R0);			/* Store zden MSW		*/
		put(R11);			/* Store LSW of znum		*/
		put(R10);			/* Store MSW of znum		*/
		CALLER_HOLD(R2)
		CALLER_SWAP
		RTLCALL (PC, ___ddiv) (DB);	/* Compute znum/zden		*/
		SAVE_OLD_FRAME(R2)
		SAVE_RET_ADDR

		alter(4);			/* Clean parameters		*/
		R10=PASS R0, R11=R1;		/* Hold z for later		*/
	
		CALL (PC, ___lib_dmult) (DB);	/* Compute w = z * z		*/
		R14=PASS R0, R15=R1;
		R12=PASS R0, R13=R1;

		R8=PASS R0, R9=R1;		/* Hold w for later		*/
		CALL (PC, ___lib_dadd) (DB);	/* Compute w + b2		*/
		R14=PASS R8, R12=mem(i_reg,1);	/* Read b2			*/
		R15=PASS R9, R13=mem(i_reg,1);

		CALL (PC, ___lib_dmult) (DB);	/* Compute (w+b2)*w		*/
		R12=PASS R0, R13=R1;		
		R14=PASS R8, R15=R9;		/* Reload w			*/

		CALL (PC, ___lib_dadd) (DB);	/* Compute (w+b2)*w+b1		*/
		R14=PASS R0, R12=mem(i_reg,1);	/* Read b1			*/
		R15=PASS R1, R13=mem(i_reg,1);

		CALL (PC, ___lib_dmult) (DB);	/* Compute ((w+b2)*w+b1)*w	*/
		R14=PASS R8, R12=R0;		
		R15=PASS R9, R13=R1;

		CALL (PC, ___lib_dadd) (DB);	/* B(w) = ((w+b2)*w+b1)*w+b0	*/
		R14=PASS R0, R12=mem(i_reg,1);
		R15=PASS R1, R13=mem(i_reg,1);

		R6=PASS R0, R7=R1;		/* Hold B(w)			*/
		CALL (PC, ___lib_dmult) (DB);	/* Compute w*a2			*/
		R14=PASS R8, R12=mem(i_reg,1);	/* Read MSW a2			*/
		R15=PASS R9, R13=mem(i_reg,1);

		CALL (PC, ___lib_dadd) (DB);	/* Compute w*a2+a1		*/
		R14=PASS R0, R12=mem(i_reg,1);
		R15=PASS R1, R13=mem(i_reg,1);

		CALL (PC, ___lib_dmult) (DB);	/* Compute (w*a2+a1)*w		*/
		R14=PASS R8, R12=R0;
		R15=PASS R9, R13=R1;

		CALL (PC, ___lib_dadd) (DB);	/* A(w) = (w*a2+a1)*w+a0	*/
		R14=PASS R0, R12=mem(i_reg,1);
		R15=PASS R1, R13=mem(i_reg,1);

		R14=PASS R8, R15=R9;		/* Move w from scratch registers*/

		put(R7);			/* Store LSW of B(w)		*/
		put(R6);			/* Store MSW of B(w)		*/
		put(R1);			/* Store A(w) LSW		*/
		put(R0);			/* Store A(w) MSW		*/
		CALLER_HOLD(R2)
		CALLER_SWAP
		RTLCALL (PC, ___ddiv) (DB);	/* Compute A(w)/B(w)		*/
		SAVE_OLD_FRAME(R2)
		SAVE_RET_ADDR

		CALL (PC, ___lib_dmult) (DB);	/* Compute r(z*z)=w*A(w)/B(w)	*/
		R12=PASS R0, R13=R1;		/* Copy A(w)/B(w), (w ==R14:R15)*/
		alter(4);			/* Clean stack			*/

		CALL (PC, ___lib_dmult) (DB);	/* Compute z*r(z*z)		*/
		R12=PASS R0, R13=R1;
		R14=PASS R10, R15=R11;

		CALL (PC, ___lib_dadd) (DB);	/* Compute R(z) = z + z*r(z*z)	*/
		R12=PASS R0, R13=R1;
		R14=PASS R10, R15=R11;

		CALL (PC, ___lib_itod) (DB);
		R10=PASS R0, R11=R1;		/* Hold R(z)			*/
		R0=dm_mdf;			/* FLOAT(N) = XN		*/

		R8=PASS R0, R9=R1;		/* Hold XN			*/
		CALL (PC, ___lib_dmult) (DB);	/* Compute XN * C1		*/
		R14=PASS R0, R12=mem(i_reg,1);
		R15=PASS R1, R13=mem(i_reg,1);
		
		CALL (PC, ___lib_dadd) (DB);	/* Compute XN * C1 + R(z)	*/
		R14=PASS R0, R12=R10;
		R15=PASS R1, R13=R11;

		R10=PASS R0, R11=R1;
		CALL (PC, ___lib_dmult) (DB);	/* Compute XN * C2		*/
		R14=PASS R8, R12=mem(i_reg,1);
		R15=PASS R9, R13=mem(i_reg,1);

		CALL (PC, ___lib_dadd) (DB);	/* Result = XN*C1+R(z) + XN*C2	*/
		R14=PASS R0, R12=R10;
		R15=PASS R1, R13=R11;

		R5=pm_mdf;			/* Get input flag		*/
		R5=PASS R5;			/* Test for flag set		*/
		IF EQ JUMP (PC, restore_state);	/* LN requested, so return	*/

		CALL (PC, ___lib_dmult) (DB);	/* Compute C3 * result		*/
		R14=PASS R0, R12=mem(i_reg,1);
		R15=PASS R1, R13=mem(i_reg,1);

restore_state:	FETCH_RETURN
		get(R15,1);
		get(R14,2);
		get(R13,3);
		get(R11,4);
		get(R10,5);
		get(R9,6);
		get(R7,7);
		get(R6,8);
		get(R5,9);
		get(R2,11);
		MR0F=R2, get(R2,12);
		MR1F=R2, get(R2,13);
		MR2F=R2, get(R3,10);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

zero_input:	ram_ireg=_errno;
		JUMP (PC, restore_state) (DB);
		rammem(ram_ireg,ram_0)=EDOM;
		R1=R1-R1, R0=dm_0;

set_edom:	ram_ireg=_errno;
		JUMP (PC, restore_state) (DB);
		rammem(ram_ireg,ram_0)=EDOM;
		R1=R1-R1, R0=dm_0;

.ENDSEG;


.SEGMENT/SPACE	Data_Space_Name;

.VAR logs_data[20] =	0x3fe6a09e,			/* C0 = sqrt(0.5)	*/
			0x667f3bcc,
			0xc041d580,			/* b2			*/
			0x4b67ce0f,
			0x40738083,			/* b1			*/
			0xfa15267e,
			0xc0880bfe,			/* b0			*/
			0x9c0d9077,
			0xbfe94415,			/* a2			*/
			0xb356db28,
			0x4030624a,			/* a1			*/
			0x2016afec,
			0xc05007ff,			/* a0			*/
			0x12b3b59b,
			0x3fe63000,			/* C1			*/
			0x00000000,
			0xbf2bd010,			/* C2			*/
			0x5c610ca8,
			0x3fdbcb7b,			/* C3			*/
			0x1526e50e;


.ENDSEG;