math.functions: implement "frexp" and support log of really big numbers. Fixes #160.

diff --git a/basis/math/functions/functions-docs.factor b/basis/math/functions/functions-docs.factor
index 414228b54e9ac0cde8671440dae47d0bc59e4107..90d6906292cd19ed4b90b01cc5dc96379bc36fba 100644 (file)
--- a/basis/math/functions/functions-docs.factor
+++ b/basis/math/functions/functions-docs.factor
@@ -110,6 +110,10 @@ HELP: exp
 { $values { "x" number } { "y" number } }
 { $description "Exponential function, " { $snippet "y=e^x" } "." } ;
 
+HELP: frexp
+{ $values { "x" number } { "y" float } { "exp" integer } }
+{ $description "Break the number " { $snippet "x" } " into a normalized fraction " { $snippet "y" } " and an integral power of 2 " { $snippet "exp" } "." $nl "The function returns a number " { $snippet "y" } " in the interval [1/2, 1) or 0, and a number " { $snippet "exp" } " such that " { $snippet "x = y*(2**exp)" } "." } ;
+
 HELP: log
 { $values { "x" number } { "y" number } }
 { $description "Natural logarithm function. Outputs negative infinity if " { $snippet "x" } " is 0." } ;

diff --git a/basis/math/functions/functions-tests.factor b/basis/math/functions/functions-tests.factor
index 507258f0d1bd730bf91f5332b08d208c5a4a243f..a439c4c41a9507904f903b2b084e3b00081126a9 100644 (file)
--- a/basis/math/functions/functions-tests.factor
+++ b/basis/math/functions/functions-tests.factor
@@ -1,5 +1,6 @@
-USING: kernel math math.constants math.functions math.order
-math.private math.libm tools.test ;
+USING: kernel math math.constants math.functions math.libm
+math.order math.ranges math.private sequences tools.test ;
+
 IN: math.functions.tests
 
 [ t ] [ 4 4 .00000001 ~ ] unit-test
@@ -37,15 +38,33 @@ IN: math.functions.tests
 [ 0 ] [ 0 3.0 ^ ] unit-test
 [ 0 ] [ 0 3 ^ ] unit-test
 
+: factorial ( n -- n! ) [ 1 ] [ [1,b] 1 [ * ] reduce ] if-zero ;
+
+[ 0.0 0 ] [ 0 frexp ] unit-test
+[ 0.5 1 ] [ 1 frexp ] unit-test
+[ -0.5 1 ] [ -1 frexp ] unit-test
+[ 0.5 2 ] [ 2 frexp ] unit-test
+[ -0.5 2 ] [ -2 frexp ] unit-test
+[ 0.64 -6 ] [ 0.01 frexp ] unit-test
+[ -0.64 -6 ] [ -0.01 frexp ] unit-test
+[ 0.75 0 ] [ 0.75 frexp ] unit-test
+[ -0.75 0 ] [ -0.75 frexp ] unit-test
+[ 1/0. 0 ] [ 1/0. frexp ] unit-test
+[ -1/0. 0 ] [ -1/0. frexp ] unit-test
+[ 0.6588418960767314 8530 t ] [ 1000 factorial [ frexp ] [ bignum? ] bi ] unit-test
+[ -0.6588418960767314 8530 t ] [ 1000 factorial neg [ frexp ] [ bignum? ] bi ] unit-test
+
 [ 0.0 ] [ 1 log ] unit-test
 [ 0.0 ] [ 1.0 log ] unit-test
 [ 1.0 ] [ e log ] unit-test
+[ 5912.128178488163 t ] [ 1000 factorial [ log ] [ bignum? ] bi ] unit-test
 
 [ 0.0 ] [ 1.0 log10 ] unit-test
 [ 1.0 ] [ 10.0 log10 ] unit-test
 [ 2.0 ] [ 100.0 log10 ] unit-test
 [ 3.0 ] [ 1000.0 log10 ] unit-test
 [ 4.0 ] [ 10000.0 log10 ] unit-test
+[ 2567.604644222133 t ] [ 1000 factorial [ log10 ] [ bignum? ] bi ] unit-test
 
 [ t ] [ 1 exp e 1.e-10 ~ ] unit-test
 [ f ] [ 1 exp 0/0. 1.e-10 ~ ] unit-test

diff --git a/basis/math/functions/functions.factor b/basis/math/functions/functions.factor
index 1b120831c0355d3ff2138c68215b629b77ded4df..865997a34d07a85375b3aba044585b912aa6198d 100644 (file)
--- a/basis/math/functions/functions.factor
+++ b/basis/math/functions/functions.factor
@@ -1,7 +1,8 @@
 ! Copyright (C) 2004, 2010 Slava Pestov.
 ! See http://factorcode.org/license.txt for BSD license.
 USING: math kernel math.constants math.private math.bits
-math.libm combinators fry math.order sequences ;
+math.libm combinators fry math.order sequences
+combinators.short-circuit math.bitwise ;
 IN: math.functions
 
 : >fraction ( a/b -- a b )
@@ -155,6 +156,23 @@ M: real absq sq ; inline
 : >=1? ( x -- ? )
     dup complex? [ drop f ] [ 1 >= ] if ; inline
 
+GENERIC: frexp ( x -- y exp )
+
+M: float frexp
+    dup { [ fp-special? ] [ zero? ] } 1|| [ 0 ] [
+        double>bits
+        [ HEX: 800f,ffff,ffff,ffff bitand 0.5 double>bits bitor bits>double ]
+        [ -52 shift 11 on-bits bitand 1022 - ] bi
+    ] if ; inline
+
+M: integer frexp
+    [ 0.0 0 ] [
+        dup 0 > [ 1 ] [ abs -1 ] if swap dup log2 [
+            52 swap - shift 52 on-bits bitand
+            0.5 double>bits bitor bits>double
+        ] [ 1 + ] bi [ * ] dip
+    ] if-zero ; inline
+
 GENERIC: log ( x -- y )
 
 M: float log dup 0.0 >= [ flog ] [ 0.0 rect> log ] if ; inline
@@ -163,6 +181,24 @@ M: real log >float log ; inline
 
 M: complex log >polar [ flog ] dip rect> ; inline
 
+<PRIVATE
+
+CONSTANT: most-positive-finite-float $[ 1/0. prev-float >integer ]
+CONSTANT: most-negative-finite-float $[ -1/0. next-float >integer ]
+
+MACRO: bignum-loghelper ( quot: ( x -- y ) -- quot )
+    dup 2 over call( x -- y ) '[
+        dup
+        most-positive-finite-float
+        most-negative-finite-float
+        between?
+        [ >float @ ] [ frexp [ @ ] [ _ * ] bi* + ] if
+    ] ;
+
+PRIVATE>
+
+M: bignum log [ log ] bignum-loghelper ;
+
 GENERIC: log1+ ( x -- y )
 
 M: object log1+ 1 + log ; inline
@@ -177,6 +213,8 @@ M: real log10 >float flog10 ; inline
 
 M: complex log10 log 10 log / ; inline
 
+M: bignum log10 [ log10 ] bignum-loghelper ;
+
 GENERIC: cos ( x -- y ) foldable
 
 M: complex cos