{
    Subroutine to compute the exponential
    of its floating point input.

	    Y=EXP(X) 

	Calling Registers
    F0 = Input Value
    l_reg = 0;

	Result Registers
   F0 = exponential of input

	Altered Registers
    F0, F1, F4, F7, F8, F10, F12
    i_reg

	Computation Time
    38 Cycles

    Version 0.03    8/6/90	Gordon A. Sterling
}
   
#include "asm_glob.h"

.SEGMENT/PM	Assembly_Library_Code_Space;

.PRECISION=MACHINE_PRECISION;

.GLOBAL	    exponential;

output_too_large:RTS (DB);
		F0=10000;			{Set some error here}
		F0=10000;

output_too_small:RTS (DB);
		F0 = .000001;
		F0 = .000001;

output_one:	RTS (DB);
		F0 = 1.0;
		F0 = 1.0;

exponential:	i_reg=exponential_data;		{Skip one cycle after this}
		F1=PASS F0;			{Copy into F1}
		F12=ABS F1, F10=mem(i_reg,1);	{Fetch maximum input}
		COMP(F1,F10), F10=mem(i_reg,1);	{Error if greater than max}
		IF GT JUMP output_too_large;	{Return XMAX with error}
		COMP(F1,F10), F10=mem(i_reg,1);	{Test for input to small}
		IF LT JUMP output_too_small;	{Return 0 with error}
		COMP(F12,F10), F10=mem(i_reg,1);{Check for output 1}
		IF LT JUMP output_one;		{Simply return 1}
		F12=F1*F10, F8=F1;		{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:	F10=F8*F8;			{Compute z=g*g}
		F7=F10*F0, F0=mem(i_reg,1);	{Compute p1*z}
		F7=F7+F0, F0=mem(i_reg,1);	{Compute p1*z + p0}
		F7=F8*F7;			{Compute g*P(z) = (p1*z+p0)*g}
		F12=F0*F10, F0=mem(i_reg,1);	{Compute q2*z}
		F12=F0+F12, F8=mem(i_reg,1);	{Compute q2*z + q1}
		F12=F10*F12;			{Compute (q2*z+q1)*z)}
		F12=F8+F12;			{Compute Q(z)=(q2*z+q1)*z+q0}
		F12=F12-F7, F10=mem(i_reg,1);	{Compute Q(z) - g*P(z)}

		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;			{R(g) = .5 + (g*P(z))/(Q(z)-g*P(z))}
		R4=FIX F4;			{Get N in fixed point again}
		RTS (DB);
		R4=R4+1;
		F0=SCALB F7 BY R4;		{R(g) * 2^(N+1)}

.ENDSEG;

.SEGMENT/SPACE	Assembly_Library_Data_Space;

.PRECISION=MEMORY_PRECISION;

.VAR exponential_data[12] = 127.0,			{BIGX }
			   -127.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;