{ 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;