Numerical evaluation of exponential and atan2 functions

check/features.frm

@@ -1225,6 +1225,27 @@ assert result("ATAN") =~ expr("
        + 7.85398163397448309615661e-01*atan(1)
 ")
 *--#] evaluate_atan : 
+*--#[ evaluate_atan2 :
+#Startfloat 21d
+CFunction atan2;
+Local ATAN2 = atan2(0,0)*atan2_(0,0)
+			+atan2(0,24)*atan2_(0,24)
+			+atan2(24/13,0)*atan2_(24/13,0)
+			+atan2(3,-1.45)*atan2_(3,-1.45)
+			+atan2(0.54321,-1.2345)*atan2_(0.54321,-1.2345)
+			+atan2(5.4321,-45/11)*atan2_(5.4321,-45/11);
+Evaluate;
+Print +s;
+.end
+#pend_if wordsize == 2
+assert succeeded?
+assert result("ATAN2") =~ expr("
+         + 1.57079632679489661923e+00*atan2(24/13,0)
+       + 2.72706541948852419832e+00*atan2(5.4321e-01, - 1.2345e+00)
+       + 2.21627784862167698618e+00*atan2(5.4321e+00, - 45/11)
+       + 2.02102192305561165724e+00*atan2(3, - 1.45e+00)
+")
+*--#] evaluate_atan2 : 
 *--#[ evaluate_sqrt :
 #Startfloat 21d
 CFunction sqrt;
@@ -1259,6 +1280,23 @@ assert result("LN") =~ expr("
        + 9.99999327347282003e-01*ln(2.71828e+00)
 ")
 *--#] evaluate_ln : 
+*--#[ evaluate_eexp :
+#Startfloat 84b
+CFunction exp;
+Local EXP = exp(0)*eexp_(0)+exp(1)*eexp_(1)+exp(-10)*eexp_(-10)+exp(2/11)*eexp_(2/11)+exp(36124.5)*eexp_(36124.5);
+Evaluate;
+Print +s;
+.end
+#pend_if wordsize == 2
+assert succeeded?
+assert result("EXP") =~ expr("
+         + 1.199396102035385909645662e+00*exp(2/11)
+       + 4.688258114961839863324201e+15688*exp(3.61245e+04)
+       + 4.539992976248485153559152e-05*exp(-10)
+       + 1.0e+00*exp(0)
+       + 2.718281828459045235360287e+00*exp(1)
+")
+*--#] evaluate_eexp : 
 *--#[ evaluate_li2 :
 #Startfloat 18d
 CFunction li2;

doc/manual/float.tex

@@ -124,9 +124,9 @@ \chapter{Floating point}
 the floating point numbers are printed in raw internal format. 
 \end{description}
 In addition to the above commands there are the following functions that 
-can be evaluated sqrt\_, ln\_, li2\_, gamma\_, agm\_, sin\_, cos\_, tan\_,
-asin\_, acos\_, atan\_, sinh\_, cosh\_, tanh\_, asinh\_, acosh\_, atanh\_.
-For the functions atan2\_ and lin\_ there is currently no code.
+can be evaluated sqrt\_, ln\_, eexp\_, li2\_, gamma\_, agm\_, sin\_, cos\_, tan\_,
+asin\_, acos\_, atan\_, atan2\_, sinh\_, cosh\_, tanh\_, asinh\_, acosh\_, atanh\_.
+For the function lin\_ there is currently no code.
 The agm\_ function is the arithmetic geometric mean of its two input 
 values.
 

doc/manual/functions.tex

@@ -1315,6 +1315,9 @@ \section{Extra reserved names}
 \leftvitem{3cm}{ln\_}\index{ln\_}\index{function!ln\_}
 \rightvitem{13cm}{The natural logarithm.}
 
+\leftvitem{3cm}{eexp\_}\index{eexp\_}\index{function!eexp\_}
+\rightvitem{13cm}{The exponential function.}
+
 \leftvitem{3cm}{sin\_}\index{sin\_}\index{function!sin\_}
 \rightvitem{13cm}{The sine function.}
 

sources/evaluate.c

@@ -373,7 +373,7 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 		if ( pars[2] == *t ) {	/* have to do this one if possible */
 TestArgument:
 /*
-			There must be a single argument, except for the AGM function
+			There must be a single argument, except for the AGM or atan2 functions
 */
 			tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);
 			if( *t == SYMBOL) {
@@ -397,7 +397,7 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 				first = 0;
 				goto nextfun;
 			}
-			if ( tt != tnext && *t != AGMFUNCTION ) goto nextfun;
+			if ( tt != tnext && *t != AGMFUNCTION && *t != ATAN2FUNCTION) goto nextfun;
 			if ( *t == SINFUNCTION ) {
 				pimul = GetPiArgument(BHEAD t+FUNHEAD);
 				if ( pimul >= 0 && pimul < 24 ) {
@@ -511,7 +511,7 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 				}
 			}
 
-			if ( *t == AGMFUNCTION ) {
+			if ( *t == AGMFUNCTION || *t == ATAN2FUNCTION ) {
 				if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;
 				if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;
 			}
@@ -531,6 +531,10 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 					else mpfr_log(auxr3,auxr1,RND);
 					mpfr_mul(auxr2,auxr2,auxr3,RND);
 				break;
+				case EXPFUNCTION:
+					mpfr_exp(auxr3,auxr1,RND);
+					mpfr_mul(auxr2,auxr2,auxr3,RND);
+				break;
 				case LI2FUNCTION: /* should be between -1 and +1 */
 					if ( nsgn == 0 ) goto getout;
 					mpfr_abs(auxr3,auxr1,RND);
@@ -615,6 +619,14 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 					mpfr_atan(auxr3,auxr1,RND);
 					mpfr_mul(auxr2,auxr2,auxr3,RND);
 				break;
+				case ATAN2FUNCTION:
+					nsgn = mpfr_sgn(auxr1);
+					nsgn2 = mpfr_sgn(auxr3);
+					// We follow the conventions of mpfr here:
+					if ( nsgn == 0 && nsgn2 >= 0) goto getout;
+					mpfr_atan2(auxr3,auxr1,auxr3,RND);
+					mpfr_mul(auxr2,auxr2,auxr3,RND);
+				break;
 				case SINFUNCTION:
 					mpfr_sin(auxr3,auxr1,RND);
 					mpfr_mul(auxr2,auxr2,auxr3,RND);
@@ -645,9 +657,11 @@ int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
 				case LNFUNCTION:
 				case LI2FUNCTION:
 				case LINFUNCTION:
+				case EXPFUNCTION:
 				case ASINFUNCTION:
 				case ACOSFUNCTION:
 				case ATANFUNCTION:
+				case ATAN2FUNCTION:
 				case SINHFUNCTION:
 				case COSHFUNCTION:
 				case TANHFUNCTION:

sources/float.c

@@ -2768,6 +2768,7 @@ int CoEvaluate(UBYTE *s)
 */
 			case SQRTFUNCTION:
 			case LNFUNCTION:
+			case EXPFUNCTION:
 			case SINFUNCTION:
 			case COSFUNCTION:
 			case TANFUNCTION:

sources/ftypes.h

@@ -477,7 +477,10 @@ typedef int (*TFUN1)();
 #define MZVHALF 121
 #define AGMFUNCTION 122
 #define GAMMAFUN 123
-#define MAXBUILTINFUNCTION 123
+#define EXPFUNCTION 124
+#define HPLFUNCTION 125
+#define MPLFUNCTION 126
+#define MAXBUILTINFUNCTION 126
 #else
 #define MAXBUILTINFUNCTION 115
 #endif

sources/inivar.h

@@ -252,6 +252,9 @@ static struct fixedfun {
 	,{"mzvhalf_"    ,0 ,0             ,0               ,0} /* MZVHALF */
 	,{"agm_"        ,0 ,0             ,0               ,0} /* AGMFUNCTION */
 	,{"gamma_"      ,0 ,0             ,0               ,0} /* GAMMAFUN */
+	,{"eexp_"       ,0 ,0             ,0               ,0} /* EXPFUNCTION */
+	,{"hpl_"        ,0 ,0             ,0               ,0} /* HPLFUNCTION */
+	,{"mpl_"        ,0 ,0             ,0               ,0} /* MPLFUNCTION */
 #endif
 };