From 0611e558739273824681943d03fd8519f9288df4 Mon Sep 17 00:00:00 2001
From: Peter Schaefer <schaeferpm@gmail.com>
Date: Wed, 27 Jun 2012 21:54:55 +0200
Subject: [PATCH] Edit:primes,euklid*,gcd,lcm,phi,sqmult*
 Added:isprime,fermat,mersenne,fermtest,solstratest,milrabtest,pepin*,lucl*,polrhofak,fermatfak

---
 UE/lib.py | 135 +++++++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 119 insertions(+), 16 deletions(-)

diff --git a/UE/lib.py b/UE/lib.py
index 068f0f8..88369c0 100644
--- a/UE/lib.py
+++ b/UE/lib.py
@@ -36,12 +36,20 @@ def primes(n):
     p = range(1, n + 1, 2)
     q = len(p);
     p[0] = 2;
-    for k in range(3, int(n ** .5 + 1), 2):
+    k = 3
+    while k < int(n ** .5 + 1):
         if p[(k - 1) / 2]:
-            for i in range(((k * k - 1) / 2), q, k):
+            i = ((k * k - 1) / 2)
+            while i < q:
                 p[i] = 0
+                i = i+k
+        k = k+2
     return filter(lambda x:x > 0, p)
 
+# Gibt zurueck ob x eine Primzahl ist
+def isprime(x):
+    return (x in primes(x))
+
 # Gibt alle Primfaktoren von n zurueck
 def factor(n):
     if n < 4:
@@ -107,12 +115,16 @@ def divisors(num):
 
 #erweiterter Euklid
 def euklid(a, b):
-    if a < b:
-        tmp = a
-        a = b
-        b = tmp
+    if not a*b:
+        return [0,0,0];
     x = [1, 0]
     y = [0, 1]
+    if a<0:
+        a = -a
+        x[0] = -1
+    if b<0:
+        b = -b
+        y[1] = -1
     q = a / b
     while a % b:
         [a , b] = [b, a - q * b]
@@ -123,12 +135,16 @@ def euklid(a, b):
 
 #erweiterter Euklid mit TeX Ausgabe
 def euklid_tex(a, b):
-    if a < b:
-        tmp = a
-        a = b
-        b = tmp
+    if not a*b:
+        return [0,0,0];
     x = [1, 0]
     y = [0, 1]
+    if a<0:
+        a = -a
+        x[0] = -1
+    if b<0:
+        b = -b
+        y[1] = -1
     q = int(a / b)
     print "\\begin{array}{ccccccc}"
     print "  r_{i-2} & r_{i-1} & q_i & x_{i-2} & x_{i-1} & y_{i-2} & y_{i-1}\\\\\\hline"
@@ -151,12 +167,27 @@ def euklid_tex(a, b):
 
 # groesster gemeinsamer Teiler 
 def gcd(a, b):
-    return euklid(a, b)[0]
+    if a*b==0:
+        return 0
+    else:
+        return euklid(a, b)[0]
 
 # kleinstes gemeinsames Vielfaches
 def lcm(a, b):
-    return abs(a * b) / gcd(a, b)
-
+    if a*b == 0:
+        return 0
+    else:
+        return abs(a * b) / gcd(a, b)
+
+# chinesischer Restsatz (x equiv a_i mod m_i)
+def chinrest(a,m):
+    M = reduce(lambda x,y:lcm(x,y),m,1)
+    x = 0
+    for i in range(len(a)):
+        s=euklid(m[i],M/m[i])[2]*M/m[i]
+        x = x+a[i]*s
+    return x % M
+   
 #eulersche phi-Funktion
 def phi(n):
     return n * reduce(lambda x, p:x * (1 - 1. / p), factorD(n).keys(), 1)
@@ -261,7 +292,6 @@ def legendre_tex(z, n):
 def jacobi(z, n):
     d = +1
     while 1:
-        print z, "/", n
         if(z > n):
             z = z % n
             continue
@@ -341,7 +371,7 @@ def jacobi_tex(z, n):
 
 
 #Square and Multiply (x,k,*,e)
-def sqmult(x,k,func,e):
+def sqmult(x,k,func=lambda x,y:x*y,e=1):
     y = e;
     while k:
         if(k%2):
@@ -352,7 +382,7 @@ def sqmult(x,k,func,e):
 
 #Square and Multiply (x,k,*,e) mit TeX Ausgabe
 def sqmult_tex(x,k,func,e):
-    y = e;
+    y = e
     print "\\begin{array}{cccc}"
     print "  x & y & k & k \\mod 2\\\\\\hline"
     print " ", x, "&" , y , "&" , k, "&" , k%2, "\\\\"
@@ -364,3 +394,76 @@ def sqmult_tex(x,k,func,e):
         print " ", x, "&" , y , "&" , k, "&" , k%2, "\\\\"
     print "\\end{array}"
     return y
+
+# Berechnet die Fermat'sche Zahl F_n
+def fermat(n):
+    return 2**(2**n)+1
+
+# Berechnet die Mersennsche Zahl M_p
+def mersenne(p):
+    return 2**p-1
+
+# Ueberprueft die Zahl n zur Basis a mit dem Fermat-Test auf Prim
+def fermtest(n,a=3):
+    return sqmult(a,n-1,lambda x,y:(x*y)%n)==1
+
+# Ueberprueft die Zahl n zur Basis a mit dem Solovay-Strassen-Test auf Prim
+def solstratest(n,a=3):
+    return sqmult(a,(n-1)/2,lambda x,y:(x*y)%n)==jacobi(a,n)%n
+
+# Ueberprueft die Zahl n zur Basis a mit dem Rabin-Miller-Test auf Prim
+def milrabtest(n,a=3):
+    s = n-1
+    t = 0
+    while not s%2:
+        s = s/2
+        t = t+1
+    if(sqmult(a,s,lambda x,y:(x*y)%n)==1):
+        return True
+    for i in range(t):
+        if sqmult(a,s*2**i,lambda x,y:(x*y)%n)==n-1:
+            return True
+    return False
+
+# Ueberprueft die FermatZahl n mit dem Pepin-Test auf Prim
+def pepin2(n):
+    return sqmult(3,(2**(2**n))/2,lambda x,y:(x*y)%(2**(2**n)+1),1) == 2**(2**n)
+
+# Ueberprueft die FermatZahl n mit dem Pepin-Test auf Prim
+def pepin(n):
+    return reduce(lambda x,y:(x**2)%2**(2**n)-int((x**2)/2**(2**n)),range(2**n-1),3)==-1
+
+# Ueberprueft die MersennscheZahl n mit dem Lucas-Lehmer-Test auf Prim
+def lucl2(p):
+    if not isprime(p):
+        return False
+    return reduce(lambda s,x:(s*s-2) % (2**p-1),range(p-2),4)==0
+
+# Ueberprueft die MersennscheZahl n mit dem Lucas-Lehmer-Test auf Prim
+def lucl(p):
+    if not isprime(p):
+        return False
+    return reduce(lambda x,y:(x**2-2)%(2**p)+int((x**2-2)/(2**p)),range(p-2),4)==2**p-1
+
+# Ueberprueft die Zahl n mit der Pollardschen rho-Methode auf Prim
+def polrhofak(n,x=3,y=3,e=2,a=1):
+    step = lambda x:sqmult(x,e,lambda x,y:(x*y) % n) +a
+    c = 0
+    while 1:
+        c =c+1
+        x = step(x)
+        y = step(step(y))
+        if gcd(x-y,n)>1:
+            return gcd(x-y,n),c
+ 
+# Ueberprueft die Zahl n mit der Fermat-Methode auf Prim
+def fermatfak(N,k=1):
+    N = k*N
+    u = int(math.ceil(N**0.5))
+    while u<=(N+1)/2:
+        if(int((u**2-N)**0.5)==(u**2-N)**0.5):
+            v = int((u**2-N)**0.5)
+            print v
+            return [u+v,u-v]
+        u = u+1
+
-- 
2.47.3