Merge pull request #535 from StephenWampler/WorkBranch

Avoid name collision with ipl, add fast factoring of large integers

Author: Jafaral

uni/lib/fastprime.icn

@@ -24,6 +24,7 @@ end
 # algorithm.
 # <[param n integer value to test for primality.]> 
 # <[param precision degree of confidence to use (defaults to 16).]>
+# <[return n if n is prime, fails otherwise.]>
 #</p>
 procedure is_prime(n, precision)
    static known_primes, k_primes_set
@@ -60,8 +61,9 @@ end
 # <[param b base]>
 # <[param e exponent]>
 # <[param m modulus]>
+# <[return (b^e)%m]>
 #</p>
-procedure pow(b,e,m)
+procedure powm(b,e,m)
    r := 1
    while e > 0 do
       if e%2 = 0 then (e /:= 2, b := (b*b)%m)
@@ -86,7 +88,53 @@ end
 # <b>Internal use only unless you know how to compute <i>d</i> and <i>s</i>.</b>
 #</p>
 procedure try_composite(n, a, d, s)
-    if pow(a,d,n) = 1 then fail
-    every i := 0 to (s-1) do if pow(a,(2^i)*d,n) = n-1 then fail
+    if powm(a,d,n) = 1 then fail
+    every i := 0 to (s-1) do if powm(a,(2^i)*d,n) = n-1 then fail
     return n
 end
+
+#<p>
+# Produce a list of all the factors of <i>n</i> using Pollard's rho algorithm.
+# Note that some factors may repeat, e.g. if <i>n = 5^2</i>.
+# <[param n number to factor]>
+# <[return list of all the factors of <i>n</i>]>
+#</p>
+procedure all_factors(n)
+   flist := []
+   f := getfactor(n)
+   if f = n then return (put(flist,n), flist)
+   flist |||:= all_factors(f) ||| all_factors(n/f)
+   suspend sort(flist)
+end
+
+#<p>
+# Generate the prime factors of <i>n</i> using Pollard's rho algorithm.
+# <[param n number to factor]>
+# <[generate the prime factors of n]>
+#</p>
+procedure gen_factors(n)
+   suspend !sort(set(all_factors(n)))
+end
+
+#<p>
+# Produce a factor of <i>n</i> using Pollard's rho algorithm.
+# <[param n number to factor]>
+# <[return one factor of <i>n</i>, or <i>n</i> if unable to find a factor.]>
+#</p>
+procedure getfactor(n)
+static sprimes
+initial sprimes := [: prime()\100 :]
+   if n == 1 then return n
+   if n%(i := !sprimes) == 0 then return i  # Speed up for small factors
+   y := x := rand_num()%(n-2) + 2
+   c := rand_num()%(n-1) + 1
+   d := 1
+   while d = 1 do {
+      x := (powm(x,2,n)+c+n)%n
+      every (y|y) := (powm(y,2,n)+c+n)%n
+      k := +(x-y)%n
+      every (z := k) *:= |k\99
+      d := gcd(z%n,n)
+      }
+   return d
+end

uni/lib/matrix_util.icn

@@ -117,7 +117,7 @@ procedure m_binop(M1, M1, op, zero:0)
    if (*M1 ~= *M2) | (*M1[1] ~= *M2[1]) then
       ::runerr(500, "nonequal matrix dimensions to m_binop")
 
-   R = m_constant(*M1, *M1[1],zero)
+   R := m_constant(*M1, *M1[1],zero)
    every (i := 1 to *M1, j := 1 to *M1[1]) do
       R[i,j] := invokeFcn(op,M1[i,j],M2[i,j])
    return R
@@ -157,7 +157,7 @@ procedure m_unaryop(M, op, zero:0)
    local R, i, j
    /op := ::proc("-", 1)
 
-   R = m_constant(*M, *M[1], zero)
+   R := m_constant(*M, *M[1], zero)
    every (i := 1 to *M, j := 1 to *M[1]) do
       R[i,j] := invokeFcn(op, M[i,j])
    return R
@@ -358,3 +358,15 @@ procedure m_write(args[])   # Matrices and output files
       }
    ::write(outfile)
 end
+
+#<p>
+# A simplification of m_write() that (a) only writes a single matrix to stdout
+# and (b) allows column width specifiction for pretty-printing the matrix.
+# <[param m matrix to output]>
+# <[param w width of each column (default: 5)]>
+#</p>
+procedure writem(m,w)
+   local i
+   /w := 5
+   every i := !m do every ::writes(::right(!i,3)|"\n")
+end