Renamed pow2 to fast-pow.

README

@@ -61,6 +61,9 @@ fast-fib:
 fast-now:
     faster, cached-per-millisecond version of "now"
 
+fast-pow:
+    faster version of pow2, pow, log, exp for doubles
+
 fizzbuzz:
     implementations of FizzBuzz problem
 
@@ -118,9 +121,6 @@ plagiarism:
 port-scan:
     simple port scanner
 
-pow:
-    faster version of 2^n
-
 power-of-2:
     various methods of implementing "power-of-2?"
 

fast-pow/fast-pow-tests.factor

@@ -0,0 +1,37 @@
+
+USING: fast-pow kernel math math.functions memory random
+sequences tools.test tools.time ;
+
+IN: fast-pow.tests
+
+[ 0.5 ] [ -1 pow2 ] unit-test
+[ 1.0 ] [ 0 pow2 ] unit-test
+[ 2.0 ] [ 1 pow2 ] unit-test
+[ 4.0 ] [ 2 pow2 ] unit-test
+[ 1024.0 ] [ 10 pow2 ] unit-test
+[ 0.0009765625 ] [ -10 pow2 ] unit-test
+[ 0.81225239593676 ] [ -0.3 pow2 ] unit-test
+[ 3.24900958374704 ] [ 1.7 pow2 ] unit-test
+
+: pow2-seq ( n -- seq )
+    [ -20 20 uniform-random-float ] replicate ;
+
+: pow2-test ( seq -- new old )
+    [ [ pow2 drop ] [ each ] benchmark ]
+    [ 2 swap [ ^ drop ] with [ each ] benchmark ] bi ;
+
+! check its at least 25% faster
+[ t ] [ 10000 pow2-seq gc pow2-test / 0.75 < ] unit-test
+
+: relative-error ( approx value -- relative-error )
+    [ - abs ] keep / ;
+
+[ t ] [
+    10000 pow2-seq
+    [ [ pow2 ] [ 2 swap ^ ] bi relative-error ] map
+    supremum 1e-9 <
+] unit-test
+
+! "orig" print [ 1000000 [ 2 16.3 ^ drop ] times ] time
+! "pow2" print [ 1000000 [ 16.3 pow2 drop ] times ] time
+

pow2/pow2.factor → fast-pow/fast-pow.factor

@@ -4,7 +4,23 @@
 USING: kernel literals math math.functions sequences
 sequences.private ;
 
-IN: pow2
+IN: fast-pow
+
+! fast-pow:
+! http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
+
+: fast-log ( x -- y )
+    double>bits -32 shift 1072632447 - 1512775 / ;
+
+: fast-exp ( x -- e^x )
+    1512775 * 1072693248 60801 - + 32 shift bits>double ;
+
+: fast-pow ( a b -- a^b )
+    [ double>bits -32 shift 1072632447 - ]
+    [ * 1072632447 + >integer 32 shift bits>double ] bi* ;
+
+! fast-pow2:
+! http://falasol.net/2-pow-x-optimization-for-double-type
 
 : float>parts ( x -- float int )
     dup >integer [ - ] keep ; inline
@@ -52,25 +68,8 @@ CONSTANT: FRAC3 $[ 2 PRECISION1 iota [ PRECISION3 / ^ ] with map ]
 : pow2 ( n -- 2^n )
     >float 2^int 2^frac * >float ;
 
-! "orig" print [ 1000000 [ 2 16.3 ^ drop ] times ] time
-! "pow2" print [ 1000000 [ 16.3 pow2 drop ] times ] time
 
-USE: random
-USE: tools.time
+! Other sources:
+! http://www.hxa.name/articles/content/fast-pow-adjustable_hxa7241_2007.html
+! http://jrfonseca.blogspot.com/2008/09/fast-sse2-pow-tables-or-polynomials.html
 
-: pow2-seq ( n -- seq )
-    [ -20 20 uniform-random-float ] replicate ;
-
-: pow2-test ( seq -- new old )
-    [ [ pow2 drop ] [ each ] benchmark ]
-    [ 2 swap [ ^ drop ] with [ each ] benchmark ] bi ;
-
-: fast-pow ( a b -- a^b )
-    [ double>bits -32 shift 1072632447 - ]
-    [ * 1072632447 + >integer 32 shift bits>double ] bi* ;
-
-: fast-exp ( x -- e^x )
-    1512775 * 1072693248 60801 - + 32 shift bits>double ;
-
-: fast-ln ( x -- y )
-    double>bits -32 shift 1072632447 - 1512775 / ;