From c31791d470104d189453eb58a279afc14f1a5cc4 Mon Sep 17 00:00:00 2001
From: Joe Groff <arcata@gmail.com>
Date: Fri, 11 Sep 2009 20:28:15 -0500
Subject: [PATCH] log1+ word (maps to C99 log1p function)

---
 basis/math/functions/functions-docs.factor | 5 +++++
 basis/math/functions/functions.factor      | 8 +++++++-
 basis/math/libm/libm.factor                | 3 +++
 3 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/basis/math/functions/functions-docs.factor b/basis/math/functions/functions-docs.factor
index 114b92ecde..134cbd398c 100644
--- a/basis/math/functions/functions-docs.factor
+++ b/basis/math/functions/functions-docs.factor
@@ -51,6 +51,7 @@ ARTICLE: "power-functions" "Powers and logarithms"
 { $subsection exp }
 { $subsection cis }
 { $subsection log }
+{ $subsection log1+ }
 { $subsection log10 }
 "Raising a number to a power:"
 { $subsection ^ }
@@ -125,6 +126,10 @@ HELP: log
 { $values { "x" number } { "y" number } }
 { $description "Natural logarithm function. Outputs negative infinity if " { $snippet "x" } " is 0." } ;
 
+HELP: log1+
+{ $values { "x" number } { "y" number } }
+{ $description "Takes the natural logarithm of " { $snippet "1 + x" } ". Outputs negative infinity if " { $snippet "1 + x" } " is zero. This word may be more accurate than " { $snippet "1 + log" } " for very small values of " { $snippet "x" } "." } ;
+
 HELP: log10
 { $values { "x" number } { "y" number } }
 { $description "Logarithm function base 10. Outputs negative infinity if " { $snippet "x" } " is 0." } ;
diff --git a/basis/math/functions/functions.factor b/basis/math/functions/functions.factor
index 92f16764c0..0cf9467795 100644
--- a/basis/math/functions/functions.factor
+++ b/basis/math/functions/functions.factor
@@ -163,7 +163,13 @@ M: float log dup 0.0 >= [ flog ] [ 0.0 rect> log ] if ; inline
 
 M: real log >float log ; inline
 
-M: complex log >polar swap flog swap rect> ; inline
+M: complex log >polar [ flog ] dip rect> ; inline
+
+GENERIC: log1+ ( x -- y )
+
+M: object log1+ 1 + log ; inline
+
+M: float log1+ dup -1.0 >= [ flog1+ ] [ 1.0 + 0.0 rect> log ] if ; inline
 
 : 10^ ( x -- y ) 10 swap ^ ; inline
 
diff --git a/basis/math/libm/libm.factor b/basis/math/libm/libm.factor
index e2bd2ef6eb..0ebc96a05c 100644
--- a/basis/math/libm/libm.factor
+++ b/basis/math/libm/libm.factor
@@ -46,6 +46,9 @@ IN: math.libm
     "double" "libm" "sqrt" { "double" } alien-invoke ;
     
 ! Windows doesn't have these...
+: flog1p ( x -- y )
+    "double" "libm" "log1p" { "double" } alien-invoke ;
+
 : facosh ( x -- y )
     "double" "libm" "acosh" { "double" } alien-invoke ;