#FEATURE
authorClaus Gittinger <cg@exept.de>
Wed, 10 Feb 2016 01:18:19 +0100
changeset 19161 a7f0570d0c3a
parent 19160 e31658d8c07c
child 19162 dcefb6a62d18
#FEATURE class: Integer added: #factorialEvenOdd #factorialHalf #factorialIter comment/format in: #factorialR changed: #factorial #factorialEO2 #factorialF (send #raisedToInteger: instead of #raisedTo:) #factorialPM
Integer.st
--- a/Integer.st	Tue Feb 09 21:49:06 2016 +0100
+++ b/Integer.st	Wed Feb 10 01:18:19 2016 +0100
@@ -751,6 +751,8 @@
     "Modified: / 15.11.1999 / 20:35:20 / cg"
 ! !
 
+
+
 !Integer class methodsFor:'class initialization'!
 
 initialize
@@ -796,6 +798,7 @@
     "Modified: 18.7.1996 / 12:26:38 / cg"
 ! !
 
+
 !Integer class methodsFor:'prime numbers'!
 
 flushPrimeCache
@@ -1163,6 +1166,7 @@
     ^ self == Integer
 ! !
 
+
 !Integer methodsFor:'*Roe'!
 
 acceptRoeVisitor: aVisitor
@@ -1433,6 +1437,7 @@
     "
 ! !
 
+
 !Integer methodsFor:'bcd conversion'!
 
 decodeFromBCD
@@ -3031,12 +3036,11 @@
 
 factorial
     "return fac(self) (i.e. 1*2*3...*self) using an iterative algorithm.
-     This is slightly faster than the recursive algorithm, and does not
-     suffer from stack overflow problems (with big receivers)"
-
-    |p i|
-
-    (self < 2) ifTrue:[
+     This chooses a good algorithm, based on the receiver.
+     Some heuristics here, which has to do with the speed of largeInteger
+     arrithmetic."
+
+    (self <= 20) ifTrue:[
         self < 0 ifTrue:[
             "/
             "/ requested factorial of a negative number
@@ -3048,8 +3052,188 @@
                 arguments:#()
                 errorString:'factorial of negative number'
         ].
-        ^ 1
+        ^ #(1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800
+          479001600 6227020800 87178291200 1307674368000 20922789888000 
+          355687428096000 6402373705728000 121645100408832000
+          2432902008176640000) at:self+1
+    ].
+    
+    self < 80000 ifTrue:[
+        ^ self factorialHalf
+    ].    
+    ^ self factorialEvenOdd
+
+    "
+     10 factorial
+     100 factorial
+     1000 factorial
+     10000 factorial
+     100000 factorial
+     200000 factorial
+     300000 factorial
+     1000000 factorial
+
+     Time millisecondsToRun:[10000 factorial]40
+     Time millisecondsToRun:[100000 factorial]3220
+     Time millisecondsToRun:[1000000 factorial]357120
+    "
+!
+
+factorialEvenOdd
+    "an recursive odd-even algorithm, which processes smaller largeInts in the loop."
+
+    |pO i s2 t stop|
+
+    (self <= 20) ifTrue:[
+        ^ #(1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800
+          479001600 6227020800 87178291200 1307674368000 20922789888000 
+          355687428096000 6402373705728000 121645100408832000
+          2432902008176640000) at:self+1
+    ].
+
+    "/
+    "/ 3 * 4 * 5 * 6 *7 * 8 .... * n
+    "/ odd numbers:
+    "/   3 5 7 9 ... n
+    "/ half even:
+    "/   2 4 6 8 ... n
+    "/   1 2 3 4 ... n//2
+    "/ is (n/2)!! << n-1
+                 
+    pO := 1.
+    i := 3.
+ 
+    "/ odds only in pairs as
+    "/      i * (i + 2)
+    "/ to get to the next pair,
+    "/      (i+4)(i+6)
+    "/ we add: 8i + 24
+    "/ (i+4)(i+6)-(i*(i+2))
+    "/ i^2 + 10i + 24 - i^2 - 2i
+    "/ 8i + 24
+    stop := self-2.
+    t := i*(i+2).
+    [i <= stop] whileTrue:[
+        "/ odd*next odd
+        pO := pO * t.
+        t := t + (i*8) + 24.
+        i := i + 4.
+    ].
+    
+    [i <= self] whileTrue:[
+        "/ odd
+        pO := pO * i.
+        i := i + 2.
     ].
+
+    "/ the factorial of the evens...
+    s2 := (self//2).
+    ^ (s2 factorialEvenOdd * pO) << s2.
+
+    "
+     (6 to:2000) conform:[:i | i factorialIter = i factorialEvenOdd]
+     
+     Time millisecondsToRun:[20000 factorialIter]
+     Time millisecondsToRun:[50000 factorialIter]
+     Time millisecondsToRun:[70000 factorialIter]
+     Time millisecondsToRun:[100000 factorialIter]
+     Time millisecondsToRun:[200000 factorialIter] 
+
+     Time millisecondsToRun:[20000 factorialEvenOdd]
+     Time millisecondsToRun:[50000 factorialEvenOdd]
+     Time millisecondsToRun:[70000 factorialEvenOdd]
+     Time millisecondsToRun:[100000 factorialEvenOdd]
+     Time millisecondsToRun:[200000 factorialEvenOdd]
+    "
+!
+
+factorialHalf
+    "an algorithm, which processes does it with half the number of multiplications.
+     this is faster than factorialPM to roughly 60000."
+
+    |p i d|
+
+    i := self.
+    self odd ifTrue:[
+        i := i - 1.
+    ].
+    
+    p := i.
+    d := i - 2.
+    [d >= 2] whileTrue:[
+        i := i + d.
+        p := p * i.
+        d := d - 2.
+    ].
+    self odd ifTrue:[
+        p := p * self
+    ].    
+    ^ p
+
+    "
+     10 factorial 3628800
+     10 factorialHalf 3628800
+     
+     11 factorial 39916800
+     11 factorialHalf 39916800
+
+     12 factorial 479001600
+     12 factorialHalf 479001600 
+
+     10000 factorial = 10000 factorialHalf
+
+     (6 to:2000) conform:[:i | i factorialIter = i factorialHalf]
+     
+     Time microsecondsToRun:[30 factorialIter]
+     Time microsecondsToRun:[30 factorialHalf]
+     Time microsecondsToRun:[50 factorialIter]
+     Time microsecondsToRun:[50 factorialHalf]
+     Time microsecondsToRun:[75 factorialIter]
+     Time microsecondsToRun:[75 factorialHalf]
+     Time microsecondsToRun:[100 factorialIter]
+     Time microsecondsToRun:[100 factorialHalf]
+     Time microsecondsToRun:[500 factorialIter]
+     Time microsecondsToRun:[500 factorialHalf]
+     Time microsecondsToRun:[1000 factorialIter]
+     Time microsecondsToRun:[1000 factorialHalf]
+     Time microsecondsToRun:[2000 factorialIter]
+     Time microsecondsToRun:[2000 factorialHalf]
+
+     Time microsecondsToRun:[500 factorial]118 120 120
+     Time microsecondsToRun:[1000 factorial]339 355 406
+     Time microsecondsToRun:[5000 factorial]15703 13669 7715
+     Time millisecondsToRun:[10000 factorial]40 30 50
+     Time millisecondsToRun:[20000 factorial]140 150 150
+     Time millisecondsToRun:[40000 factorial]600 570 560 570
+     Time millisecondsToRun:[60000 factorial]1220 1240 1340
+     Time millisecondsToRun:[80000 factorial]2600 2580 2540
+     Time millisecondsToRun:[100000 factorial]4680 4810 5280
+     Time millisecondsToRun:[120000 factorial]8100 8010 7920
+     Time millisecondsToRun:[150000 factorial]13830 14040 13360
+     Time millisecondsToRun:[200000 factorial]23880 23740 
+
+     Time microsecondsToRun:[500 factorialHalf]150 142 192
+     Time microsecondsToRun:[1000 factorialHalf]383 527 684
+     Time microsecondsToRun:[5000 factorialHalf]6654 9221 4629
+     Time millisecondsToRun:[10000 factorialHalf]20 30 20
+     Time millisecondsToRun:[20000 factorialHalf]110 110 110
+     Time millisecondsToRun:[40000 factorialHalf]490 490 490
+     Time millisecondsToRun:[60000 factorialHalf]1100 1090 1070
+     Time millisecondsToRun:[80000 factorialHalf]1920 1920 1880
+     Time millisecondsToRun:[100000 factorialHalf]3030 3010 3000
+     Time millisecondsToRun:[120000 factorialHalf]4830 4770 4760
+     Time millisecondsToRun:[150000 factorialHalf]14510 13940 13900
+     Time millisecondsToRun:[200000 factorialHalf]28730 28160 
+    "
+!
+
+factorialIter
+    "return fac(self) (i.e. 1*2*3...*self) using an iterative algorithm.
+     This is slightly faster than the recursive algorithm, and does not
+     suffer from stack overflow problems (with big receivers)"
+
+    |p i|
+
     p := 2.
     i := 3.
     [i <= self] whileTrue:[
@@ -3076,7 +3260,7 @@
     "return fac(self) (i.e. 1*2*3...*self) using a recursive algorithm.
 
      This is included to demonstration purposes - if you really need
-     factorial numbers, use the iterative #factorial, which is slightly
+     factorial numbers, use the iterative #factorial, which is
      faster and does not suffer from stack overflow problems (with big receivers)."
 
     (self >= 2) ifTrue:[
@@ -4574,6 +4758,7 @@
     "Created: / 09-01-2012 / 17:18:06 / cg"
 ! !
 
+
 !Integer methodsFor:'special modulo arithmetic'!
 
 add_32:anInteger