math.functions: implement ldexp.

basis/math/functions/functions-docs.factor

@@ -116,6 +116,11 @@ 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 "e^" } "." $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: ldexp
+{ $values { "x" number } { "exp" number } { "y" number } }
+{ $description "Multiply " { $snippet "x" } " by " { $snippet "2^exp" } "." }
+{ $notes { $link ldexp } " is the inverse of " { $link frexp } "." } ;
+
 HELP: log
 { $values { "x" number } { "y" number } }
 { $description "Natural logarithm function. Outputs negative infinity if " { $snippet "x" } " is 0." } ;

basis/math/functions/functions-tests.factor

@@ -55,6 +55,16 @@ IN: math.functions.tests
 [  0.75 10,002 t ] [  3 10,000 2^ * [ frexp ] [ bignum? ] bi ] unit-test
 [ -0.75 10,002 t ] [ -3 10,000 2^ * [ frexp ] [ bignum? ] bi ] unit-test
 
+{ 0.0 } [ 0.0 1 ldexp ] unit-test
+{ -0.0 } [ -0.0 1 ldexp ] unit-test
+{ 1/0. } [ 1/0. 1 ldexp ] unit-test
+{ -1/0. } [ -1/0. 1 ldexp ] unit-test
+{ t } [ NAN: 90210 dup 1 ldexp [ fp-nan-payload ] same? ] unit-test
+{ 49152.0 } [ 12.0 12 ldexp ] unit-test
+{ 0x1.8p-9 } [ 12.0 -12 ldexp ] unit-test
+{ 49152 } [ 12 12 ldexp ] unit-test
+{ 0 } [ 12 -12 ldexp ] unit-test
+
 [ 0.0 ] [ 1 log ] unit-test
 [ 0.0 ] [ 1.0 log ] unit-test
 [ 1.0 ] [ e log ] unit-test

basis/math/functions/functions.factor

@@ -204,6 +204,28 @@ M: integer frexp
         ] [ 1 + ] bi [ * ] dip
     ] if-zero ; inline
 
+DEFER: copysign
+
+GENERIC# ldexp 1 ( x exp -- y )
+
+M: float ldexp
+    over fp-special? [ over zero? ] unless* [ drop ] [
+        [ double>bits dup -52 shift 0x7ff bitand 1023 - ] dip +
+        {
+            { [ dup -1074 < ] [ drop 0 copysign ] }
+            { [ dup 1023 > ] [ drop 0 < -1/0. 1/0. ? ] }
+            [
+                dup -1022 < [ 52 + -52 2^ ] [ 1 ] if
+                [ -0x7ff0,0000,0000,0001 bitand ]
+                [ 1023 + 52 shift bitor bits>double ]
+                [ * ] tri*
+            ]
+        } cond
+    ] if ;
+
+M: integer ldexp
+    2dup [ zero? ] either? [ 2drop 0 ] [ shift ] if ;
+
 GENERIC: log ( x -- y )
 
 M: float log dup 0.0 >= [ flog ] [ 0.0 rect> log ] if ; inline