[genius] Tue Jun 11 17:19:55 2013 Jiri (George) Lebl <jirka 5z com>

[George Lebl <jirka src gnome org>]

commit 7bdd5f6e48ece567cb47bcb49c0092c3e000dad8
Author: Jiri (George) Lebl <jirka 5z com>
Date:   Tue Jun 11 17:20:05 2013 -0500

    Tue Jun 11 17:19:55 2013  Jiri (George) Lebl <jirka 5z com>
    
        * lib/number_theory/factoring.gel: Improve Factors so that for larger
          numbers it uses "Factorize" rather than a slow naive
          implementation.
    
        * lib/number_theory/primes.gel: minor mersenne update

 ChangeLog                       |    8 +++++
 lib/number_theory/factoring.gel |   66 +++++++++++++++++++++++++-------------
 lib/number_theory/primes.gel    |    6 ++--
 3 files changed, 54 insertions(+), 26 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 8291a3c..780385e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+Tue Jun 11 17:19:55 2013  Jiri (George) Lebl <jirka 5z com>
+
+       * lib/number_theory/factoring.gel: Improve Factors so that for larger
+         numbers it uses "Factorize" rather than a slow naive
+         implementation.
+
+       * lib/number_theory/primes.gel: minor mersenne update
+
 Tue Jun 11 16:12:48 2013  Jiri (George) Lebl <jirka 5z com>
 
        * src/gnome-genius.c: add completion to "help on function"
diff --git a/lib/number_theory/factoring.gel b/lib/number_theory/factoring.gel
index ad07549..6a5f777 100644
--- a/lib/number_theory/factoring.gel
+++ b/lib/number_theory/factoring.gel
@@ -67,39 +67,59 @@ function CombineFactorizations(a,b) = (
 )
 
 # Returns all factors of n, i.e., all numbers between 1 and n that divide
-# n.
-# FIXME: this is a naive implementation. Better would be to compute
-# the prime factorization for n and then combine it in all possible ways.
+# n.  For small n does a naive implementation by trying all factors
+# up to square root, for large numbers does Factorize first and then
+# combines the factors in all possible ways.
 SetHelp("Factors", "number_theory", "Return all factors of a number")
-function Factors(n) =
- (
+function Factors(n) = (
   if not IsInteger(n) then
        (error("Factors: argument not an integer");bailout);
   if n == 0 then return null;
 
-  front_list=[0];
-  back_list=[0];
+  list = null;
 
   if n < 0 then (
     n = -n;
-    front_list=[0,-1]
+    list = [-1]
   );
 
-  for loop = 1 to floor(sqrt(n-1)) do
-   (
-    if Divides(loop,n) then
-     (
-      front_list=[front_list,loop];
-      back_list=[n/loop,back_list]
-     )
-   );
-
-  if IsPerfectSquare(n) then
-    r = [front_list,sqrt(n),back_list]
-  else
-    r = [front_list,back_list];
-  r@(2:elements(r)-1)
- )
+  # for small n do naive implementation which is faster
+  if n < 100000 then (
+    back_list=null;
+
+    for loop = 1 to floor(sqrt(n-1)) do (
+      if Divides(loop,n) then (
+        list=[list,loop];
+        back_list=[n/loop,back_list]
+      )
+    );
+
+    if IsPerfectSquare(n) then
+      r = [list,sqrt(n),back_list]
+    else
+      r = [list,back_list];
+    r
+  ) else (
+    fnc = fn = Factorize (n);
+    for j=2 to columns(fnc) do fnc@(2,j) = 0;
+  
+    do (
+      list = [list, (prod j=2 to columns(fnc) do fnc@(1,j)^fnc@(2,j))];
+    
+      gotnext = false;
+      for j = 2 to columns(fnc) do (
+        if fnc@(2,j) < fn@(2,j) then (
+          increment fnc@(2,j);
+          gotnext = true;
+          break
+        ) else ( #fnc@(2,j) == fn@(2,j)
+          fnc@(2,j) := 0
+        )
+      )
+    ) while (gotnext); 
+    SortVector(list)
+  ) 
+)
 
 SetHelp("FermatFactorization", "number_theory", "Attempt Fermat factorization of n into (t-s)*(t+s), returns 
t and s as a vector if possible, null otherwise")
 function FermatFactorization(n,tries) = (
diff --git a/lib/number_theory/primes.gel b/lib/number_theory/primes.gel
index 1068ad4..de598a9 100644
--- a/lib/number_theory/primes.gel
+++ b/lib/number_theory/primes.gel
@@ -121,10 +121,10 @@ function IsMersennePrimeExponent(p) = (
                return true;
 
        # http://www.mersenne.org / GIMPS doublechecked everything up
-       # to 26,083,693 (on Feb. 6, 2013)
-       if p <= 26083693 then
+       # to 26,114,633 (on June 11, 2013)
+       if p <= 26114633 then
                return false;
 
-       error("IsMersennePrimeExponent: Number too large (known values up to: " + 26083693 + ")");
+       error("IsMersennePrimeExponent: Number too large (known values up to: " + 26114633 + ")");
        bailout
 );