/*
	This C language subroutine computes the hyperbolic
	sine or hyperbolic cosine 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 6/10/90
	Updated on 5/94 by AS

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

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

.SEGMENT/CODE	Code_Space_Name;
.FILE RTL_FILENAME;

.EXTERN		_expf;

.GLOBAL		_sinhf, _coshf;

_sinhf:		put(R3);
		put(R7);
		F7=-1.0;		/*Set SINH flag	*/
		F1=-F7, put(R11);	/* F1=1.0	*/

		COMP(F4,F1);		/*Test for Y > 1.0 */
		IF GT JUMP (PC, cosh_path); /*Follow cosh if Y > 1.0*/

		F1=eps;			/*Test for small input*/
		COMP(F4,F1), F0=F4;	/*If small, return input*/
		IF LT JUMP (PC, restore_state) (DB);
		F3=F4*F4, FETCH_RETURN	/*Compute f=x*x */
		F1=p1;
		F7=F1*F3;		/*Compute p1*x */
		F1=p0;
		F7=F7+F1;		/*Compute P(f)=p1*x+p0 */

		F1=q0;
		CALL (PC, ___float_divide) (DB);
		F12=F3+F1;		/*Compute Q(f)=f+q0 */
		F11=2.0;

		F7=F7*F3;		/*Compute R(f)=f*P(f)/Q(f)*/
		JUMP (PC, restore_state) (DB);
		F0=F7*F4;		/*Compute X*R(f)*/
		F0=F0+F4;		/*Compute X+X*R(f)*/


_coshf:		put(R3);
		put(R7);
		F3=ABS F4, put(R11);

		F7=1.0;			/*Set COSH flag*/
		
cosh_path:	F1=ybar;   
		COMP(F3,F1);		/*Test for Y > ybar*/
		IF LE JUMP (PC, exp_path);
		F1=ln_v;
		F4=F3-F1;		/*W=Y-ln(v)*/
		F1=wmax;
		COMP(F0,F1);
		IF GT JUMP (PC, error_return);
	    	CALLER_HOLD(R2)
		CALLER_SWAP
		RTLCALL (PC, _expf) (DB);
		SAVE_OLD_FRAME(R2)
		SAVE_RET_ADDR

		F1=v_over_2_minus_1;
		F4=F1*F0;		/*Compute (v/2-1.0)*Z*/
		F0=F4+F0;		/*Compute Result=Z+(v/2-1.0)*Z*/
tst_sign:	F3=PASS F3, FETCH_RETURN/*Test for sign of result*/
		IF LT F0=-F0;		/*Negate result if necessary*/

restore_state:	get(R11,1);
		get(R7,2);
		get(R3,3);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

exp_path:	R4=PASS R3, CALLER_HOLD(R2)
		CALLER_SWAP
		RTLCALL (PC, _expf) (DB);
		SAVE_OLD_FRAME(R2)
		SAVE_RET_ADDR

		CALL (PC, ___float_divide) (DB);
		F12=PASS F0, F4=F0;
		F11=2.0;

		JUMP (PC, tst_sign) (DB);
		F7=F7+F4, R12=dm_M1;	/*Compute Z +/- 1/Z */
		F0=SCALB F7 BY R12;	/*Compute (Z +/- 1.0/Z)/2 */

error_return:	JUMP (PC, restore_state) (DB);
		F0=0;
		FETCH_RETURN

.ENDSEG;