use radix literals

[factor.git] / basis / math / functions / functions.factor

diff --git a/basis/math/functions/functions.factor b/basis/math/functions/functions.factor
index d91b4b6b92a0c5904d418037be7a6999e573322c..f1e68c9f86efd6ddc00aac19d8bd8f1898725900 100644 (file)
--- a/basis/math/functions/functions.factor
+++ b/basis/math/functions/functions.factor
@@ -1,7 +1,7 @@
-! Copyright (C) 2004, 2008 Slava Pestov.
+! 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 math.order sequences ;
+math.libm combinators fry math.order sequences ;
 IN: math.functions
 
 : >fraction ( a/b -- a b )
@@ -13,12 +13,13 @@ IN: math.functions
 GENERIC: sqrt ( x -- y ) foldable
 
 M: real sqrt
-    >float dup 0.0 < [ neg fsqrt 0.0 swap rect> ] [ fsqrt ] if ; inline
+    >float dup 0.0 <
+    [ neg fsqrt [ 0.0 ] dip rect> ] [ fsqrt ] if ; inline
 
 : factor-2s ( n -- r s )
     #! factor an integer into 2^r * s
     dup 0 = [ 1 ] [
-        0 swap [ dup even? ] [ [ 1 + ] [ 2/ ] bi* ] while
+        [ 0 ] dip [ dup even? ] [ [ 1 + ] [ 2/ ] bi* ] while
     ] if ; inline
 
 <PRIVATE
@@ -26,13 +27,13 @@ M: real sqrt
 GENERIC# ^n 1 ( z w -- z^w ) foldable
 
 : (^n) ( z w -- z^w )
-    make-bits 1 [ [ dupd * ] when [ sq ] dip ] reduce nip ; inline
+    make-bits 1 [ [ over * ] when [ sq ] dip ] reduce nip ; inline
 
 M: integer ^n
     [ factor-2s ] dip [ (^n) ] keep rot * shift ;
 
 M: ratio ^n
-    [ >fraction ] dip [ ^n ] curry bi@ / ;
+    [ >fraction ] dip '[ _ ^n ] bi@ / ;
 
 M: float ^n (^n) ;
 
@@ -62,7 +63,7 @@ M: float exp fexp ; inline
 
 M: real exp >float exp ; inline
 
-M: complex exp >rect swap fexp swap polar> ; inline
+M: complex exp >rect [ exp ] dip polar> ; inline
 
 <PRIVATE
 
@@ -81,26 +82,25 @@ M: complex exp >rect swap fexp swap polar> ; inline
 : real^? ( x y -- ? )
     2dup [ real? ] both? [ drop 0 >= ] [ 2drop f ] if ; inline
 
-: 0^ ( x -- z )
-    [ 0/0. ] [ 0 < 1/0. 0 ? ] if-zero ; inline
+: 0^ ( zero x -- z )
+    swap [ 0/0. ] swap '[ 0 < 1/0. _ ? ] if-zero ; inline
 
-: (^mod) ( n x y -- z )
-    make-bits 1 [
-        [ dupd * pick mod ] when [ sq over mod ] dip
-    ] reduce 2nip ; inline
+: (^mod) ( x y n -- z )
+    [ make-bits 1 ] dip dup
+    '[ [ over * _ mod ] when [ sq _ mod ] dip ] reduce nip ; inline
 
 : (gcd) ( b a x y -- a d )
     over zero? [
         2nip
     ] [
         swap [ /mod [ over * swapd - ] dip ] keep (gcd)
-    ] if ;
+    ] if ; inline recursive
 
 PRIVATE>
 
 : ^ ( x y -- z )
     {
-        { [ over 0 = ] [ nip 0^ ] }
+        { [ over zero? ] [ 0^ ] }
         { [ dup integer? ] [ integer^ ] }
         { [ 2dup real^? ] [ [ >float ] bi@ fpow ] }
         [ ^complex ]
@@ -109,13 +109,13 @@ PRIVATE>
 : nth-root ( n x -- y ) swap recip ^ ; inline
 
 : gcd ( x y -- a d )
-    [ 0 1 ] 2dip (gcd) dup 0 < [ neg ] when ; foldable
+    [ 0 1 ] 2dip (gcd) dup 0 < [ neg ] when ; inline
 
 : lcm ( a b -- c )
     [ * ] 2keep gcd nip /i ; foldable
 
 : divisor? ( m n -- ? )
-    mod 0 = ;
+    mod 0 = ; inline
 
 ERROR: non-trivial-divisor n ;
 
@@ -125,11 +125,8 @@ ERROR: non-trivial-divisor n ;
     [ non-trivial-divisor ] if ; foldable
 
 : ^mod ( x y n -- z )
-    over 0 < [
-        [ [ neg ] dip ^mod ] keep mod-inv
-    ] [
-        -rot (^mod)
-    ] if ; foldable
+    over 0 <
+    [ [ [ neg ] dip ^mod ] keep mod-inv ] [ (^mod) ] if ; foldable
 
 GENERIC: absq ( x -- y ) foldable
 
@@ -158,6 +155,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? [ dup zero? ] unless* [ 0 ] [
+        double>bits
+        [ 0x800f,ffff,ffff,ffff bitand 0.5 double>bits bitor bits>double ]
+        [ -52 shift 0x7ff bitand 1022 - ] bi
+    ] if ; inline
+
+M: integer frexp
+    [ 0.0 0 ] [
+        dup 0 > [ 1 ] [ abs -1 ] if swap dup log2 [
+            52 swap - shift 0x000f,ffff,ffff,ffff 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
@@ -166,6 +180,29 @@ M: real log >float log ; inline
 
 M: complex log >polar [ flog ] dip rect> ; inline
 
+<PRIVATE
+
+: most-negative-finite-float ( -- x )
+    -0x1.ffff,ffff,ffff,fp1023 >integer ; inline
+: most-positive-finite-float ( -- x )
+    0x1.ffff,ffff,ffff,fp1023 >integer ; inline
+CONSTANT: log-2   0x1.62e42fefa39efp-1
+CONSTANT: log10-2 0x1.34413509f79ffp-2
+
+: (representable-as-float?) ( x -- ? )
+    most-negative-finite-float
+    most-positive-finite-float between? ; inline
+
+: (bignum-log) ( n log-quot: ( x -- y ) log-2 -- log )
+    [ dup ] dip '[
+        dup (representable-as-float?)
+        [ >float @ ] [ frexp [ @ ] [ _ * ] bi* + ] if
+    ] call ; inline
+
+PRIVATE>
+
+M: bignum log [ log ] log-2 (bignum-log) ;
+
 GENERIC: log1+ ( x -- y )
 
 M: object log1+ 1 + log ; inline
@@ -180,6 +217,8 @@ M: real log10 >float flog10 ; inline
 
 M: complex log10 log 10 log / ; inline
 
+M: bignum log10 [ log10 ] log10-2 (bignum-log) ;
+
 GENERIC: cos ( x -- y ) foldable
 
 M: complex cos