PR#4638: added expm1 and log1p.

Cleaned up some doc comments in stdlib/pervasives.mli


git-svn-id: http://caml.inria.fr/svn/ocaml/trunk@9216 f963ae5c-01c2-4b8c-9fe0-0dff7051ff02

6 files changed, 88 insertions, 21 deletions

diff --git a/byterun/floats.c b/byterun/floats.c
index c708cbe0d1..21ba411a77 100644
--- a/byterun/floats.c
+++ b/byterun/floats.c
@@ -248,11 +248,8 @@ CAMLprim value caml_log10_float(value f)
 
 CAMLprim value caml_modf_float(value f)
 {
-#if __SC__
-  _float_eval frem;       /* Problem with Apple's <math.h> */
-#else
   double frem;
-#endif
+
   CAMLparam1 (f);
   CAMLlocal3 (res, quo, rem);
 
@@ -329,6 +326,46 @@ CAMLprim value caml_ceil_float(value f)
   return caml_copy_double(ceil(Double_val(f)));
 }
 
+/* These emulations of expm1() and log1p() are due to William Kahan.
+   See http://www.plunk.org/~hatch/rightway.php */
+
+CAMLexport double caml_expm1(double x)
+{
+#ifdef HAS_EXPM1_LOG1P
+  return expm1(x);
+#else
+  double u = exp(x);
+  if (u == 1.)
+    return x;
+  if (u - 1. == -1.)
+    return -1.;
+  return (u - 1.) * x / log(u);
+#endif
+}
+
+CAMLexport double caml_log1p(double x)
+{
+#ifdef HAS_EXPM1_LOG1P
+  return log1p(x);
+#else
+  double u = 1. + x;
+  if (u == 1.)
+    return x;
+  else
+    return log(u) * x / (u - 1.);
+#endif
+}
+
+CAMLprim value caml_expm1_float(value f)
+{
+  return caml_copy_double(caml_expm1(Double_val(f)));
+}
+
+CAMLprim value caml_log1p_float(value f)
+{
+  return caml_copy_double(caml_log1p(Double_val(f)));
+}
+
 CAMLprim value caml_eq_float(value f, value g)
 {
   return Val_bool(Double_val(f) == Double_val(g));
diff --git a/config/s-templ.h b/config/s-templ.h
index bbbffb4fa3..bc4d9643f9 100644
--- a/config/s-templ.h
+++ b/config/s-templ.h
@@ -52,6 +52,11 @@
 /* Define SUPPORT_DYNAMIC_LINKING if dynamic loading of C stub code
    via dlopen() is available. */
 
+#define HAS_EXPM1_LOG1P
+
+/* Define HAS_EXPM1_LOG1P if the math functions expm1() and log1p()
+   are available.  (Standard C99 but not C89.) */
+
 /* 2. For the Unix library. */
 
 #define HAS_SOCKETS
diff --git a/configure b/configure
index 7d2fe7aeef..1c35ee2e53 100755
--- a/configure
+++ b/configure
@@ -808,7 +808,14 @@ else
   fi
 fi
 
-# For the sys module
+# For the Pervasives module
+
+if sh ./trycompile expm1.c $mathlib; then
+  echo "expm1() and log1p() found."
+  echo "#define HAS_EXPM1_LOG1P" >> s.h
+fi
+
+# For the Sys module
 
 if sh ./hasgot getrusage; then
   echo "getrusage() found."
diff --git a/otherlibs/threads/pervasives.ml b/otherlibs/threads/pervasives.ml
index 4e3c3fe76f..09c5381414 100644
--- a/otherlibs/threads/pervasives.ml
+++ b/otherlibs/threads/pervasives.ml
@@ -87,6 +87,7 @@ external ( *. ) : float -> float -> float = "%mulfloat"
 external (/.) : float -> float -> float = "%divfloat"
 external ( ** ) : float -> float -> float = "caml_power_float" "pow" "float"
 external exp : float -> float = "caml_exp_float" "exp" "float"
+external expm1 : float -> float = "caml_expm1_float" "caml_expm1" "float"
 external acos : float -> float = "caml_acos_float" "acos" "float"
 external asin : float -> float = "caml_asin_float" "asin" "float"
 external atan : float -> float = "caml_atan_float" "atan" "float"
@@ -95,6 +96,7 @@ external cos : float -> float = "caml_cos_float" "cos" "float"
 external cosh : float -> float = "caml_cosh_float" "cosh" "float"
 external log : float -> float = "caml_log_float" "log" "float"
 external log10 : float -> float = "caml_log10_float" "log10" "float"
+external log1p : float -> float = "caml_log1p_float" "caml_log1p" "float"
 external sin : float -> float = "caml_sin_float" "sin" "float"
 external sinh : float -> float = "caml_sinh_float" "sinh" "float"
 external sqrt : float -> float = "caml_sqrt_float" "sqrt" "float"
diff --git a/stdlib/pervasives.ml b/stdlib/pervasives.ml
index 9e34cf2fd4..94fe485178 100644
--- a/stdlib/pervasives.ml
+++ b/stdlib/pervasives.ml
@@ -83,6 +83,7 @@ external ( *. ) : float -> float -> float = "%mulfloat"
 external (/.) : float -> float -> float = "%divfloat"
 external ( ** ) : float -> float -> float = "caml_power_float" "pow" "float"
 external exp : float -> float = "caml_exp_float" "exp" "float"
+external expm1 : float -> float = "caml_expm1_float" "caml_expm1" "float"
 external acos : float -> float = "caml_acos_float" "acos" "float"
 external asin : float -> float = "caml_asin_float" "asin" "float"
 external atan : float -> float = "caml_atan_float" "atan" "float"
@@ -91,6 +92,7 @@ external cos : float -> float = "caml_cos_float" "cos" "float"
 external cosh : float -> float = "caml_cosh_float" "cosh" "float"
 external log : float -> float = "caml_log_float" "log" "float"
 external log10 : float -> float = "caml_log10_float" "log10" "float"
+external log1p : float -> float = "caml_log1p_float" "caml_log1p" "float"
 external sin : float -> float = "caml_sin_float" "sin" "float"
 external sinh : float -> float = "caml_sinh_float" "sinh" "float"
 external sqrt : float -> float = "caml_sqrt_float" "sqrt" "float"
diff --git a/stdlib/pervasives.mli b/stdlib/pervasives.mli
index 7ead634fbb..72567af001 100644
--- a/stdlib/pervasives.mli
+++ b/stdlib/pervasives.mli
@@ -264,45 +264,59 @@ external log : float -> float = "caml_log_float" "log" "float"
 external log10 : float -> float = "caml_log10_float" "log10" "float"
 (** Base 10 logarithm. *)
 
+external expm1 : float -> float = "caml_expm1_float" "caml_expm1" "float"
+(** [expm1 x] computes [exp x -. 1.0], giving numerically-accurate results
+    even if [x] is close to [0.0]. *)
+
+external log1p : float -> float = "caml_log1p_float" "caml_log1p" "float"
+(** [log1p x] computes [log(1.0 +. x)] (natural logarithm),
+    giving numerically-accurate results even if [x] is close to [0.0]. *)
+
 external cos : float -> float = "caml_cos_float" "cos" "float"
-(** See {!Pervasives.atan2}. *)
+(** Cosine.  Argument is in radians. *)
 
 external sin : float -> float = "caml_sin_float" "sin" "float"
-(** See {!Pervasives.atan2}. *)
+(** Sine.  Argument is in radians. *)
 
 external tan : float -> float = "caml_tan_float" "tan" "float"
-(** See {!Pervasives.atan2}. *)
+(** Tangent.  Argument is in radians. *)
 
 external acos : float -> float = "caml_acos_float" "acos" "float"
-(** See {!Pervasives.atan2}. *)
+(** Arc cosine.  The argument must fall within the range [[-1.0, 1.0]]. 
+    Result is in radians and is between [0.0] and [pi]. *)
 
 external asin : float -> float = "caml_asin_float" "asin" "float"
-(** See {!Pervasives.atan2}. *)
+(** Arc sine.  The argument must fall within the range [[-1.0, 1.0]]. 
+    Result is in radians and is between [-pi/2] and [pi/2]. *)
 
 external atan : float -> float = "caml_atan_float" "atan" "float"
-(** See {!Pervasives.atan2}. *)
+(** Arc tangent.  
+    Result is in radians and is between [-pi/2] and [pi/2]. *)
 
 external atan2 : float -> float -> float = "caml_atan2_float" "atan2" "float"
-(** The usual trigonometric functions. *)
+(** [atan x y] returns the arc tangent of [y /. x].  The signs of [x]
+    and [y] are used to determine the quadrant of the result.
+    Result is in radians and is between [-pi] and [pi]. *)
 
 external cosh : float -> float = "caml_cosh_float" "cosh" "float"
-(** See {!Pervasives.tanh}. *)
+(** Hyperbolic cosine. *)
 
 external sinh : float -> float = "caml_sinh_float" "sinh" "float"
-(** See {!Pervasives.tanh}. *)
+(** Hyperbolic sine. *)
 
 external tanh : float -> float = "caml_tanh_float" "tanh" "float"
-(** The usual hyperbolic trigonometric functions. *)
+(** Hyperbolic tangent. *)
 
 external ceil : float -> float = "caml_ceil_float" "ceil" "float"
-(** See {!Pervasives.floor}. *)
+(** Round above to an integer value.
+    [ceil f] returns the least integer value greater than or equal to [f].
+    The result is returned as a float. *)
 
 external floor : float -> float = "caml_floor_float" "floor" "float"
-(** Round the given float to an integer value.
-   [floor f] returns the greatest integer value less than or
-   equal to [f].
-   [ceil f] returns the least integer value greater than or
-   equal to [f]. *)
+(** Round below to an integer value.
+    [floor f] returns the greatest integer value less than or
+    equal to [f].
+    The result is returned as a float. *)
 
 external abs_float : float -> float = "%absfloat"
 (** Return the absolute value of the argument. *)