// exp().

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

// Specification.
#include "cl_real.h"


// Implementation.

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

#undef MAYBE_INLINE
#define MAYBE_INLINE inline
#include "cl_LF_minusp.cc"
#include "cl_LF_exponent.cc"

// Division durch ln(2).
inline cl_F_div_t cl_floor_ln2 (const cl_F& x)
{
	// Bei 0<=x<1/2 kann man sofort q:=0 setzen.
	if (!minusp(x) && (float_exponent(x) < 0))
		return cl_F_div_t(0,x);
	else
		return floor2(x,cl_ln2(x));
}
inline cl_LF_div_t cl_floor_ln2 (const cl_LF& x)
{
	// Bei 0<=x<1/2 kann man sofort q:=0 setzen.
	if (!minusp(x) && (float_exponent(x) < 0))
		return cl_LF_div_t(0,x);
	else
		return floor2(x,The(cl_LF)(cl_ln2(x)));
}

cl_R exp (const cl_R& x)
{
// Methode:
// x rational -> bei x=0 1 als Ergebnis, sonst x in Float umwandeln.
// x Float ->
//   d := (float-digits x),
//   Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erhöhen,
//   (q,r) := (floor x ln(2))
//   Ergebnis ist exp(q*ln(2)+r) = (scale-float exp(r) q).

	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;
	}
	// Rechengenauigkeit erhöhen und durch ln(2) dividieren:
	if (longfloatp(xx) && (TheLfloat(xx)->len >= 160)) {
		DeclareType(cl_LF,xx);
		var cl_LF_div_t q_r = cl_floor_ln2(extend(xx,TheLfloat(xx)->len+1));
		var cl_I& q = q_r.quotient;
		var cl_LF& r = q_r.remainder;
		return cl_float(scale_float(expx_ratseries(r),q),xx);
	} else {
		var cl_F_div_t q_r = cl_floor_ln2(cl_F_extendsqrtx(xx));
		var cl_I& q = q_r.quotient;
		var cl_F& r = q_r.remainder;
		return cl_float(scale_float(expx_naive(r),q),xx);
	}
}
