/*
	This C language subroutine computes the exponential
	value of its 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

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

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;
 
.GLOBAL	    _expf;

_expf:		i_reg=exponential_data;		/*Skip one cycle after this*/
		R2=PASS R4, put(R10);		/*Hold input*/
		F12=ABS F2, F1=mem(i_reg,1);	/*Fetch maximum input*/
		COMP(F2,F1), F1=mem(i_reg,1);	/*Error if greater than max*/
		IF GT JUMP (PC, output_too_large); /*Return XMAX with error*/
		COMP(F2,F1), F1=mem(i_reg,1);	/*Test for input to small*/
		IF LT JUMP (PC, output_too_small); /*Return 0 with error*/
		COMP(F12,F1), F1=mem(i_reg,1);/*Check for output 1*/
		IF LT JUMP (PC, output_one);	/*Simply return 1*/
		F12=F2*F1, F8=F2;		/*Compute N = X/ln(C)*/
		R4=FIX F12;			/*Round to nearest*/
		F4=FLOAT R4, F0=mem(i_reg,1);	/*Back to floating point*/

compute_g:	F12=F0*F4, F0=mem(i_reg,1);	/*Compute XN*C1*/
		F0=F0*F4, F12=F8-F12;		/*Compute |X|-XN*C1, and XN*C2*/
		F8=F12-F0, F0=mem(i_reg,1);	/*Compute g=(|X|-XN*C1)-XN*C2*/

compute_R:	F1=F8*F8;		        /*Compute z=g*g*/
		F2=F1*F0, F0=mem(i_reg,1);	/*Compute p1*z*/
		F2=F2+F0, F0=mem(i_reg,1);	/*Compute p1*z + p0*/
		F2=F8*F2;		        /*Compute g*P(z) = (p1*z+p0)*g*/
		F12=F0*F1, F0=mem(i_reg,1);	/*Compute q2*z*/
		F12=F0+F12, F8=mem(i_reg,1);	/*Compute q2*z + q1*/
		F12=F1*F12;			/*Compute (q2*z+q1)*z)*/
		F12=F8+F12;		        /*Compute Q(z)=(q2*z+q1)*z+q0*/
		F12=F12-F2, F10=mem(i_reg,1);	/*Compute Q(z) - g*P(z)*/

		F0=RECIPS F12, F1=F4;		/*Get 4 bit seed R0=1/D*/
		F12=F0*F12, F4=F2;		/*D' = D*R0*/
		F4=F0*F4, F0=F10-F12;		/*F0=R1=2-D', F7=N*R0*/
		F12=F0*F12;			/*F12=D'=D'*R1*/
		F4=F0*F4, F0=F10-F12;		/*F7=N*R0*R1, F0=R2=2-D'*/
		F12=F0*F12;			/*F12=D'=D'*R2*/
		F2=F0*F4, F0=F10-F12;		/*F7=N*R0*R1*R2, F0=R3=2-D'*/
		F2=F0*F2;			/*F7=N*R0*R1*R2*R3*/
		F2=F2+F8;			/*R(g) = .5 + (g*P(z))/(Q(z)-g*P(z))*/
		R4=FIX F1;			/*Get N in fixed point again*/
		R4=R4+1, FETCH_RETURN
		F0=SCALB F2 BY R4, get(R10,1);	/*R(g) * 2^(N+1)*/

restore_state:	RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

output_too_large:JUMP (PC, restore_state) (DB);
		F0=10000;			/*Set some error here*/
		FETCH_RETURN

output_too_small:JUMP (PC, restore_state) (DB);
		R0 = 0x00800000;
		FETCH_RETURN

output_one:	JUMP (PC, restore_state) (DB);
		F0 = 1.0;
		FETCH_RETURN

.ENDSEG;

.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR exponential_data[12] = 127.0,			/*BIGX */
			   -88.0,			/*SMALLX*/
			      0.000000001,		/*eps*/
			      1.4426950408889634074,	/*1/ln(2.0)*/
			      0.693359375,		/*C1*/
			     -2.1219444005469058277E-4, /*C2*/
			      0.59504254977591E-2,	/*P1*/
			      0.24999999999992,		/*P0*/
			      0.29729363682238E-3,	/*Q2*/
			      0.53567517645222E-1,	/*Q1*/
			      0.5,			/*Q0 and others*/
			      2.0;			/*Used in divide*/

.ENDSEG;