log1+ word (maps to C99 log1p function)

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." } ;

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
 

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 ;