// cl_cosh_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_cosh_sinh_t cl_cosh_sinh (const cl_R& x)
{
// Methode:
// x rational -> bei x=0 (1,0) als Ergebnis, sonst x in Float umwandeln.
// x Float -> Genauigkeit erhöhen,
//   e := Exponent aus (decode-float x), d := (float-digits x)
//   falls x=0.0 oder e<=(1-d)/2 liefere (1.0,x)
//     (denn bei e<=(1-d)/2 ist
//      1 <= sinh(x)/x < cosh(x) = 1+x^2/2+... < 1+2^(-d),
//      also ist cosh(x), auf d Bits gerundet, gleich 1.0
//      und sinh(x), auf d Bits gerundet, gleich x).
//   falls e<0:
//     y:=(sinh(x)/x)^2 errechnen,
//     cosh(x) = sqrt(1+x^2*y) und sinh(x) = x*sqrt(y) errechnen.
//   falls e>=0: y:=exp(x) errechnen,
//     (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden.
//   Genauigkeit wieder verringern.

	var cl_F xx;
	if (rationalp(x)) {
		DeclareType(cl_RA,x);
		if (zerop(x)) // x=0 -> (1,0) als Ergebnis
			return cl_cosh_sinh_t(1,0);
		xx = cl_float(x); // sonst in Float umwandeln
	} else {
		DeclareType(cl_F,x);
		xx = x;
	}
	// x Float
	var sintL e = float_exponent(xx);
	if (e < 0) { // Exponent e abtesten
		// e<0
		if (zerop(xx) || (e <= (1-(sintL)float_digits(xx))>>1))
			// e <= (1-d)/2 <==> e <= -ceiling((d-1)/2)
			return cl_cosh_sinh_t(cl_float(1,xx),xx);
		// 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_cosh_sinh_t(
					cl_float(hyp.cosh,xx),
					cl_float(hyp.sinh,xx)
				       );
			} else
			#endif
			if (TheLfloat(xx)->len >= 585) {
				// verwende exp(x), schneller als coshsinh_ratseries
				var cl_LF xxx = extend(xx,TheLfloat(xx)->len+ceiling((uintL)(-e),intDsize));
				var cl_F y = exp(xxx);
				var cl_F y_inv = recip(y);
				return cl_cosh_sinh_t(
					cl_float(scale_float(y + y_inv, -1), xx),
					cl_float(scale_float(y - y_inv, -1), xx)
				       );
			} else {
				var cl_LF xxx = The(cl_LF)(cl_F_extendsqrt(xx));
				var cl_LF y = sinhx_naive(xxx);
				var cl_LF z = sqrt(y);
				if (minusp(xxx))
					z = -z;
				return cl_cosh_sinh_t(
					cl_float(sqrt(1+y),xx), // sqrt(1+y)
					cl_float(z,xx)
				       );
			}
		} else {
			var cl_F xxx = cl_F_extendsqrt(xx);
			var cl_F y = sinhxbyx_naive(xxx);
			return cl_cosh_sinh_t(
				cl_float(sqrt(1+square(xxx)*y),xx), // sqrt(1+x^2*y)
				cl_float(xxx*sqrt(y),xx)
			       );
		}
	} else {
		// e>=0 -> verwende exp(x)
		var cl_F y = exp(xx);
		var cl_F y_inv = recip(y);
		return cl_cosh_sinh_t(
			scale_float(y+y_inv,-1),
			scale_float(y-y_inv,-1)
		       );
	}
}
