Merge pull request #532 from StephenWampler/fastprime

Support large integer primes, fix typos in matrix_util.icn

Author: Jafaral

uni/lib/fastprime.icn

@@ -0,0 +1,92 @@
+#<p>
+# Efficient support for large primes.
+#</p>
+
+link factors   # Need to get <i>prime</i> generator
+
+#<p>
+# Generate all primes >= <i>n</i>.
+# <[param n start looking for primes >= <i>n</> (defaults to 2)]>
+#</p>
+procedure genprimes(n)
+   /n := 2
+   if n = 2 then suspend 2
+   if n%2 = 0 then n +:= 1
+   suspend is_prime(seq(n,2))
+end
+
+#<p>
+# Is <i>n</i> probably prime at least to some large probability?
+# Fully accurate for <i>n</i> < 341550071728321.  Accuracy above that
+# value depends on value of <i>precision</i>.  The default value for
+# <i>precision</i> is generally sufficient.
+# Google the <b>Miller-Rabin primality test</b> for details on the
+# algorithm.
+# <[param n integer value to test for primality.]> 
+# <[param precision degree of confidence to use (defaults to 16).]>
+#</p>
+procedure is_prime(n, precision)
+   static known_primes, k_primes_set
+   initial {
+      /precision := 16
+      every put(known_primes := [], prime()\17)
+      k_primes_set := set(known_primes)
+      }
+    /precision := 16
+      
+    if n <= 1 then fail
+    # Fast check for small n
+    if member(k_primes_set, n) then return n
+    if (n % known_primes[1 to precision]) = 0 then fail
+
+    # d and s are needed for primetest()
+    (d := n-1, s := 0)
+    while d%2 = 0 do (d /:= 2, s +:= 1)
+    # Returns exact according to http://primes.utm.edu/prove/prove2_3.html
+    if n < 1373653 then return  primetest(n, known_primes[1:2], d, s)
+    if n < 25326001 then return primetest(n, known_primes[1:3], d, s)
+    if n < 118670087467 then
+        if n ~= 3215031751 then return primetest(n, known_primes[1:4], d, s)
+        else fail
+    if n < 2152302898747 then return primetest(n, known_primes[1:5], d, s)
+    if n < 3474749660383 then return primetest(n, known_primes[1:6], d, s)
+    if n < 341550071728321 then return primetest(n, known_primes[1:7], d, s)
+    # otherwise
+    return primetest(n, known_primes[1:precision], d, s)
+end
+
+#<p>
+# Efficient implementation of <i>(b^e)%m</i>.
+# <[param b base]>
+# <[param e exponent]>
+# <[param m modulus]>
+#</p>
+procedure pow(b,e,m)
+   r := 1
+   while e > 0 do
+      if e%2 = 0 then (e /:= 2, b := (b*b)%m)
+      else (e -:= 1, r := (r*b)%m, e /:= 2, b := (b*b)%m)
+   return r
+end
+
+#<p>
+# Is <i>n</i> probably prime, based on the known primes in <i>A</i>?
+#  (Used to simplify the code in <b>is_prime</b> as it just wraps
+#   and inverts <b>try_composite</b>.)
+# <b>Internal use only unless you know how to compute <i>d</i> and <i>s</i>.</b>
+#</p>
+procedure primetest(n,A,d,s)
+   if try_composite(n,!A,d,s) then fail
+   return n
+end
+
+#<p>
+# Succeed if <i>n</i> might be a composite, based solely on the value
+# of <i>a</i>.
+# <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
+    return n
+end

uni/lib/matrix_util.icn

@@ -28,10 +28,10 @@ end
 #<p>
 #  <[param m number of rows]>
 #  <[param n number of columns]>
-#  <[param x initial value of every element]>
+#  <[param x initial value of every element (defaults to 0)]>
 #  <[returns a <tt>m</tt> x <tt>n<tt> matrix with <tt>x</tt> everywhere]>
 #</p>
-procedure m_constant(m,n,x:0.0)
+procedure m_constant(m,n,x:0)
    local A
    /n := m
 
@@ -44,11 +44,11 @@ end
 #  Creates a m x n matrix with ones on the diagonal and zeros everywhere else
 #  <[param m number of rows]>
 #  <[param n number of columns]>
-#  <[param zero (<i>optional</i>) can be used in place of <tt>0.0</tt>]>
-#  <[param one (<i>optional</i>) can be used in place of <tt>1.0</tt>]>
+#  <[param zero (<i>optional</i>) can be used in place of <tt>0</tt>]>
+#  <[param one (<i>optional</i>) can be used in place of <tt>1</tt>]>
 #  <[returns <tt>m</tt>X<tt>n</tt> identify matrix]>
 #</p>
-procedure m_identity(m,n, zero:0.0,one:1.0)
+procedure m_identity(m,n, zero:0,one:1)
    local A, i
    /n := m
 
@@ -106,17 +106,18 @@ end
 #  <[param M1 first matrix]>
 #  <[param M2 second matrix]>
 #  <[param op binary operation to invoke]>
+#  <[param zero (<i>optional</i>) can be used in place of <tt>0</tt>]>
 #  <[returns <tt>M1 op M2</tt>, with <tt>op</tt> defaulting to
 #    <tt>proc("+",2)</tt>]>
 #</p>
-procedure m_binop(M1, M1, op)
+procedure m_binop(M1, M1, op, zero:0)
    local R, i, j
    /op := ::proc("+", 2)
 
    if (*M1 ~= *M2) | (*M1[1] ~= *M2[1]) then
       ::runerr(500, "nonequal matrix dimensions to m_binop")
 
-   R = m_constant(*M1, *M1[1])
+   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
@@ -148,14 +149,15 @@ end
 #  Generic unary operation across a matrix.
 #  <[param M matrix]>
 #  <[param op unary operation to invoke]>
+#  <[param zero (<i>optional</i>) can be used in place of <tt>0</tt>]>
 #  <[returns <tt>op M</tt>, with <tt>op</tt> defaulting to
 #            <tt>proc("-",1)</tt>]>
 #<p>
-procedure m_unaryop(M, op)
+procedure m_unaryop(M, op, zero:0)
    local R, i, j
    /op := ::proc("-", 1)
 
-   R = m_constant(*M, *M[1])
+   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
@@ -185,7 +187,7 @@ procedure m_multiply(M1, M2, addident:0, addfcn, mulfcn)
    if (*M1[1] ~= *M2) then
       ::runerr(500, "incompatible matrix dimensions to m_multiply")
 
-   R := m_constant(*M1, *M2[1])
+   R := m_constant(*M1, *M2[1], addident)
    every (i := 1 to *R, j := 1 to *R[1]) do {
       R[i,j] := addident
       every k := 1 to *M2 do
@@ -204,7 +206,7 @@ end
 #  by overriding various functions and constants used internally.
 #</p>
 procedure m_lupDecomposition(M, subfcn, mulfcn, divfcn,
-                             invertMetric, addident:0.0, mulident:1.0)
+                             invertMetric, addident:0, mulident:1)
    local n, i, b, k, k2, p, j, L, U
    /subfcn := ::proc("-", 2)
    /mulfcn := ::proc("*", 2)
@@ -252,7 +254,7 @@ end
 #  Additional optional arguments allow you to customize the operation
 #  by overriding various functions and constants used internally.
 #</p>
-procedure m_lupSolve(L,U,p,b,addfun,subfun,mulfun,divfun,addident:0.0)
+procedure m_lupSolve(L,U,p,b,addfun,subfun,mulfun,divfun,addident:0)
    local n, i, j, y, x
    /addfun := ::proc("+", 2)
    /subfun := ::proc("-", 2)
@@ -346,7 +348,7 @@ procedure m_write(args[])   # Matrices and output files
       case ::type(arg) of {
          "file": outfile := arg
          "list": {    # Assume it's a matrix!
-            every i := !M do {
+            every i := !arg do {
                every ::writes(outfile, !i," ")
                ::write(outfile)
                }