// cosh().

// 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 cosh (const cl_R& x)
{
// Methode:
// x rational -> bei x=0 1 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
//     (denn bei e<=(1-d)/2 ist 1 <= cosh(x) = 1+x^2/2+... < 1+2^(-d),
//      also ist cosh(x), auf d Bits gerundet, gleich 1.0).
//   falls e<0:
//     y := x/2 = (scale-float x -1), (sinh(y)/y)^2 errechnen,
//     cosh(x) = 1+x*y*(sinh(y)/y)^2 errechnen.
//   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 -> 1 als Ergebnis
			return 1;
		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))
			return cl_float(1,xx);
		var uintL d = float_digits(xx);
		if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
			return cl_float(1,xx); // ja -> 1.0 als Ergebnis
		// 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.cosh,xx);
			} else
			#endif
			if (TheLfloat(xx)->len >= 600) {
				// verwende exp(x), schneller als coshsinh_ratseries
				var cl_LF xxx = extend(xx,TheLfloat(xx)->len+1);
				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));
				var cl_LF y = scale_float(xxx,-1);
				// 1 + 2*sinh(y)^2, und wieder runden
				return cl_float(1 + scale_float(sinhx_naive(y),1), xx);
			}
		} else {
			var cl_F xxx = cl_F_extendsqrt(xx);
			var cl_F y = scale_float(xxx,-1);
			// 1 + 2*y^2*(sinh(y)/y)^2, und wieder runden
			return cl_float(1 + scale_float(square(y) * sinhxbyx_naive(y),1), 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 three algorithms, on an i486 33 MHz, running Linux,
// applied to x = sqrt(2)-1 = 0.414...
//   N      naive  ratseries exp&recip
//   10     0.008   0.037     0.012
//   25     0.032   0.117     0.047
//   50     0.11    0.33      0.017
//  100     0.40    1.06      0.63
//  250     2.65    5.2       3.3
//  500    11.1    18.7      11.5
// 1000    46      61        35
// 2500   238     250       143
// ==> exp&recip fastest for N >= 600.
