diff --git a/basis/math/functions/functions-docs.factor b/basis/math/functions/functions-docs.factor index 414228b54e..90d6906292 100644 --- 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 507258f0d1..a439c4c41a 100644 --- 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 1b120831c0..865997a34d 100644 --- 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 +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