// sinh().

// General includes.
#include "cl_sysdep.h"

// Specification.
#include "cl_real.h"


// Implementation.

#include "cl_R_tran.h"
#include "cl_F.h"
#include "cl_lfloat.h"
#include "cl_LF.h"

cl_R sinh (const cl_R& x)
{
// Methode:
// x rational -> bei x=0 0 als Ergebnis, sonst x in Float umwandeln.
// x Float -> Genauigkeit erhöhen,
//   e := Exponent aus (decode-float x)
//   falls e<0: (sinh(x)/x)^2 errechnen, Wurzel ziehen, mit x multiplizieren.
//   falls e>=0: y:=exp(x) errechnen, (scale-float (- y (/ y)) -1) bilden.

	var cl_F xx;
	if (rationalp(x)) {
		DeclareType(cl_RA,x);
		if (zerop(x)) // x=0 -> 0 als Ergebnis
			return 0;
		xx = cl_float(x); // sonst in Float umwandeln
	} else {
		DeclareType(cl_F,x);
		xx = x;
	}
	// x Float
	if (float_exponent(xx) < 0) { // Exponent e abtesten
		// e<0
		// Rechengenauigkeit erhöhen
		if (longfloatp(xx)) {
			DeclareType(cl_LF,xx);
			#if 0
			if (TheLfloat(xx)->len >= infty) {
				var cl_LF xxx = extend(xx,TheLfloat(xx)->len+1);
				var cl_LF_cosh_sinh_t hyp = coshsinh_ratseries(xxx);
				return cl_float(hyp.sinh,xx);
			} else
			#endif
			if ((TheLfloat(xx)->len >= 500)
			    && (float_exponent(xx) > (-(sintL)float_digits(xx))>>1)) {
				// verwende exp(x), schneller als coshsinh_ratseries
				// (aber nur bei 0 > e > -d/2, denn wir müssen, um
				// Auslöschung zu verhindern, |e| Bits dazunehmen)
				var cl_LF xxx = extend(xx,TheLfloat(xx)->len+ceiling((uintL)(-float_exponent(xx)),intDsize));
				var cl_F y = exp(xxx);
				var cl_F z = scale_float(y - recip(y), -1); // (/ (- y (/ y)) 2)
				return cl_float(z,xx);
			} else {
				var cl_LF xxx = The(cl_LF)(cl_F_extendsqrt(xx));
				// Wurzel aus sinh(x)^2 bilden
				var cl_LF z = sqrt(sinhx_naive(xxx));
				if (minusp(xxx))
					z = -z;
				return cl_float(z,xx);
			}
		} else {
			var cl_F xxx = cl_F_extendsqrt(xx);
			// Wurzel aus (sinh(x)/x)^2 mit x multiplizieren und wieder runden
			return cl_float(sqrt(sinhxbyx_naive(xxx))*xxx,xx);
		}
	} else {
		// e>=0 -> verwende exp(x)
		var cl_F y = exp(xx);
		return scale_float(y - recip(y), -1); // (/ (- y (/ y)) 2)
	}
}

// Timings of the two algorithms, on an i486 33 MHz, running Linux,
// applied to x = sqrt(2)-1 = 0.414...
//   N      naive  ratseries
//   10     0.008   0.037
//   25     0.034   0.115
//   50     0.13    0.33
//  100     0.50    1.07
//  250     3.3     5.2
//  500    14.2    18.8
// 1000    59      61
// 2500   297     247
// ==> ratseries faster for N >= 1300.
