CompAlg -- Sage

import random

def to_number(a): n = 0 f = 1 for i in range(len(a)): n += a[i] * f f <<= 64 return n

def to_polynomial(a): poly = [] while a > 0: poly.append(a%(1<<64)) a >>= 64 return poly

def subtract(a, b): c = a[:] carry = 0 for i in range(len(b)): if c[i]< b[i]: c[i] += (1<<64) j = i+1 while c[j] == 0: c[j] = (1<<64)-1 c[j] -= 1 c[i] -= b[i] return c

def add(a, b): modulus = (1<<64) d = [] if len(b) > len(a): temp = b b = a a = temp carry = 0 for i in range(len(a)): t = 0 if i < len(b): t = b[i] c = a[i] + t + carry carry = c >> 64 d.append(c % modulus) if carry != 0: d.append(carry) return d

Het klassieke vermenigvuldigingsalgoritme werkt door elk woord met elk woord te vermenigvuldigen, en toe te voegen aan het juiste woord in het resultaat. Er kan hier mogelijks overflow plaatsvinden, die we dan overbrengen naar het volgende woord.

def classic_multiplication(a, b): modulus = (1<<64) if len(b) > len(a): temp = b b = a a = temp d = [0]*(len(a)+len(b)) for i in range(len(d)): for j in range(max(i - len(a)+1, 0), min(len(b), i+1)): d[i] += b[j] * a[i-j] k = i while d[k] >= modulus: d[k+1] += d[k] >> 64 d[k] = d[k] % modulus k += 1 if d[-1] == 0: d = d[:-1] return d

Het algoritme van Karatsuba herschrijft de vermenigvuldiging van $a*b = (a_1*x^{n/2}+a_0)(b_1*x^{n/2}+b_0) = a_1*b_1*x^n + a_0*b_0 + (a_1*b_0+a_0*b_1)*x^2{n/2} = a_1*b_1*x^n + a_0*b_0 + ((a_1+a_0)(b_1+b_0) -a_1*b_1 a_0*b_0) *x^2{n/2}$. Hierdoor zijn er slechts drie vermenigvuldigingen nodig, en is de tijdscomplexiteit $O(n^{log_2 3})$.

def karatsuba(a, b): if len(b) > len(a): temp = b b = a a = temp if len(b) == 1: return classic_multiplication(a, b) a1 = a[len(a)//2:] a0 = a[:len(a)//2] b1 = [0] b0 = b if len(b) > len(a)//2: b1 = b[len(a)//2:] b0 = b[:len(a)//2] a12 = add(a1, a0) b12 = add(b1, b0) f = karatsuba(a1, b1) g = karatsuba(a0, b0) h = karatsuba(a12, b12) h = subtract(h, add(f,g)) return add(add(([0]*(len(a)//2*2) + f), ([0]*(len(a)//2)+h)), g)

for k in range(1,16): i = 2**k a = [0]*i b = [0]*i for j in range(i): a[j] = floor(random.random()*(100)) b[j] = floor(random.random()*(100)) time1 = cputime() res1 = classic_multiplication(a, b) time1 = cputime() - time1 time2 = cputime() res2 = karatsuba(a, b) time2 = cputime() - time2 print(res1==res2, i, time1, time2)

(True, 2, 0.00012300000003051537, 0.00019199999996999395)
(True, 4, 8.500000012645614e-05, 0.00042500000017753337)
(True, 8, 0.00016699999991942605, 0.0013279999998303538)
(True, 16, 0.0003979999999046413, 0.004126999999925829)
(True, 32, 0.0011730000001080043, 0.012536000000181957)
(True, 64, 0.0038919999999507127, 0.03792399999997542)
(True, 128, 0.013846000000057757, 0.11475199999995311)
(True, 256, 0.05231499999990774, 0.3351199999999608)
(True, 512, 0.19084700000007615, 0.5938400000000001)
(True, 1024, 0.42531700000017736, 1.5718580000000202)
(True, 2048, 1.696279000000004, 4.728260000000091)
(True, 4096, 6.679397999999992, 14.18733599999996)
(True, 8192, 26.168273, 42.64646799999991)
(True, 16384, 103.19446300000004, 128.25840799999992)
(True, 32768, 409.68479, 384.0914730000004)

(True, 2, 0.00012300000003051537, 0.00019199999996999395)
(True, 4, 8.500000012645614e-05, 0.00042500000017753337)
(True, 8, 0.00016699999991942605, 0.0013279999998303538)
(True, 16, 0.0003979999999046413, 0.004126999999925829)
(True, 32, 0.0011730000001080043, 0.012536000000181957)
(True, 64, 0.0038919999999507127, 0.03792399999997542)
(True, 128, 0.013846000000057757, 0.11475199999995311)
(True, 256, 0.05231499999990774, 0.3351199999999608)
(True, 512, 0.19084700000007615, 0.5938400000000001)
(True, 1024, 0.42531700000017736, 1.5718580000000202)
(True, 2048, 1.696279000000004, 4.728260000000091)
(True, 4096, 6.679397999999992, 14.18733599999996)
(True, 8192, 26.168273, 42.64646799999991)
(True, 16384, 103.19446300000004, 128.25840799999992)
(True, 32768, 409.68479, 384.0914730000004)

We zien duidelijk dat het klassieke vermenigvuldiginsalgoritme sneller is voor kleine inputs, maar dat het kwadratisch stijgt en voor grotere inputs trager wordt.

Onderstaande functie kijkt na of de orde van een element gelijk is aan een parameter orde. Factors wordt meegegeven zodat het niet voor elke oproep opnieuw berekent moet worden.

Het kijkt na of $element^{order}$ gelijk is aan $1$, zodat de orde van het element zeker een deler van order is.

Daarna kijkt het of de orde niet kleiner is dan order, door $order/factor$ als macht te proberen voor elke factor van order.

def check_order(element, order, factors): if element**(order) != 1: return False for (f,_exp) in factors: if element ** (order/f) == 1: return False return True

De Lucas Test of N-1 test, test of n een priemgetal is door te kijken of er een getal $g$ bestaat zodat de orde van $g$ over $\mathbb{Z}/n\mathbb{Z}$ gelijk is aan $n-1$.

def lucasTest(n, maxIter=None): if n == 2: return True ring = IntegerModRing(n) factors = list((n-1).factor()) iter = 0 while iter != maxIter: iter += 1 g = ring.random_element() if check_order(g, n-1, factors): return True return False

De N+1 test, test of n een priemgetal is door te proberen een $t$ te vinden zodat $x$ orde $n+1$ heeft over $R_t = (\mathbb{Z}/n\mathbb{Z})[x]/(x^2-tx+1)$.

def nplus1test(n, maxIter = None): if n == 2: return True nz = IntegerModRing(n) ring = PolynomialRing(nz, 'x', sparse=True) R.<x> = ring factors = list((n+1).factor()) iter = 0 while iter != maxIter: iter += 1 t = nz.random_element() I = x**2 - t*x+1 Q = R.quotient(I) if check_order(Q(x), n+1, factors): return True return False

nplus1test((1<<1279)-1)

True

True

# Test multiplication algorithms for i in range(1, 100): a = [0]*i b = [0]*i for j in range(i): a[j] = floor(random.random()*(100)) b[j] = floor(random.random()*(100)) res1 = classic_multiplication(a, b) res2 = karatsuba(a, b) res = to_number(a) * to_number(b) assert to_number(res1)==res, "Classic multiplication failed. %d * %d" % (to_number(a), to_number(b)) assert to_number(res2)==res, "Karatsuba multiplication failed. %d * %d" % (to_number(a), to_number(b))

# Test primality test algorithms for i in range(2,101): assert lucasTest(i, 100) == is_prime(i), "Lucas Test failed for number %d." % (i) assert nplus1test(i, 100) == is_prime(i), "N+1 test failed for number %d." % (i) if i % 10 == 0: print("Tested first %d numbers succesfully." % (i))

Tested first 10 numbers succesfully.
Tested first 20 numbers succesfully.
Tested first 30 numbers succesfully.
Tested first 40 numbers succesfully.
Tested first 50 numbers succesfully.
Tested first 60 numbers succesfully.
Tested first 70 numbers succesfully.
Tested first 80 numbers succesfully.
Tested first 90 numbers succesfully.
Tested first 100 numbers succesfully.

Tested first 10 numbers succesfully.
Tested first 20 numbers succesfully.
Tested first 30 numbers succesfully.
Tested first 40 numbers succesfully.
Tested first 50 numbers succesfully.
Tested first 60 numbers succesfully.
Tested first 70 numbers succesfully.
Tested first 80 numbers succesfully.
Tested first 90 numbers succesfully.
Tested first 100 numbers succesfully.

CompAlg

1808 days ago by David.Vandorpe