math.functions, speed up truncate for floats

basis/math/functions/functions-tests.factor

@@ -1,5 +1,5 @@
 USING: kernel literals math math.constants math.functions math.libm
-math.order math.ranges math.private sequences tools.test ;
+math.order math.ranges math.private sequences tools.test math.floats.env ;
 
 IN: math.functions.tests
 
@@ -168,6 +168,12 @@ CONSTANT: log10-factorial-1000 0x1.40f3593ed6f8ep11
 { -4.0 } [ -4.0 floor ] unit-test
 { -4.0 } [ -4.0 ceiling ] unit-test
 
+! first floats without fractional part
+{ 0x1.0p52 } [ 0x1.0p52 truncate ] unit-test
+{ 0x1.0000000000001p52 } [ 0x1.0000000000001p52 truncate ] unit-test
+{ -0x1.0p52 } [ -0x1.0p52 truncate ] unit-test
+{ -0x1.0000000000001p52 } [ -0x1.0000000000001p52 truncate ] unit-test
+
 { -5 } [ -9/2 round ] unit-test
 { -4 } [ -22/5 round ] unit-test
 { 5 } [ 9/2 round ] unit-test
@@ -188,6 +194,12 @@ CONSTANT: log10-factorial-1000 0x1.40f3593ed6f8ep11
 { t } [ -0.3 round double>bits -0.0 double>bits = ] unit-test
 { t } [ 0.3 round double>bits 0.0 double>bits = ] unit-test
 
+! A signaling NaN should raise an exception
+{ { +fp-invalid-operation+ } } [ [ NAN: 4000000000000 truncate drop ] collect-fp-exceptions ] unit-test
+{ { +fp-invalid-operation+ } } [ [ NAN: 4000000000000 round drop ] collect-fp-exceptions ] unit-test
+{ { +fp-invalid-operation+ } } [ [ NAN: 4000000000000 ceiling drop ] collect-fp-exceptions ] unit-test
+{ { +fp-invalid-operation+ } } [ [ NAN: 4000000000000 floor drop ] collect-fp-exceptions ] unit-test
+
 { 6 59967 } [ 3837888 factor-2s ] unit-test
 { 6 -59967 } [ -3837888 factor-2s ] unit-test
 

basis/math/functions/functions.factor

@@ -348,10 +348,30 @@ M: real atan >float atan ; inline
 
 : acot ( x -- y ) recip atan ; inline
 
-: truncate ( x -- y ) dup dup 1 mod - over float? [
-    over [ -1.0 > ] [ 0.0 < ] bi and
-    [ swap copysign ] [ nip ] if
- ] [ nip ] if ; inline
+GENERIC: truncate ( x -- y )
+
+M: real truncate dup 1 mod - ;
+
+M: float truncate
+    dup double>bits
+    dup -52 shift 0x7ff bitand 0x3ff -
+    ! check for floats without fractional part (>= 2^52)
+    dup 52 < [
+        [ drop ] 2dip
+        dup 0 < [
+            ! the float is between -1.0 and 1.0,
+            ! the result is +/-0.0
+            drop -63 shift zero? 0.0 -0.0 ?
+        ] [
+            ! Put zeroes in the correct part of the mantissa
+            0x000fffffffffffff swap neg shift bitnot bitand
+            bits>double
+        ] if
+    ] [
+        ! check for nans and infinities and do an operation on them
+        ! to trigger fp exceptions if necessary
+        nip 0x400 = [ dup + ] when
+    ] if ; inline
 
 GENERIC: round ( x -- y )