import warnings
import random
warnings.simplefilter('ignore', DeprecationWarning)
"""
Algorytm Euklidesa
"""
def gcd(a,b):
while b: a, b = b, a % b
return a
"""
Zwraca rozwiazania rownania x**2 = n mod p
"""
def solve_xx_n_mod_p(n, p):
for x in range(p):
if x*x % p == n % p:
return [x, p-x]
return []
"""
Zwraca listę liczb B-gładkich postaci x**2 - n
"""
def Sieving_B_smooth_in_x2_n(B, n, K):
numbers = [[x*x - n, x] for x in range(ceil(sqrt(n)), ceil(sqrt(n)) + K)]
for i in range(len(numbers)):
x = numbers[i]
while x[0] % 2 == 0: x[0] /= 2
for p in list(primes(3,B)):
if p > K: break
if n % p == 0:
print "---> Znaleziono dzielnik %s: %s" % (n, p)
s = solve_xx_n_mod_p(n, p) if n % p != 0 else [p*x for x in solve_xx_n_mod_p(n/p, p)]
if s == None or len(s) < 2: continue
for i in [0] + range(1, K/p + 1):
for _s in s:
_i = i*p + (_s - ceil(sqrt(n))) % p
try:
x = numbers[_i]
while x[0] % p == 0: x[0] /= p
except IndexError:
pass
return [y*y - n for [x,y] in numbers if x == 1]
"""
Tworzy drzewo z n w węźle i o poddrzewach left i right.
Dla takiego drzewa T obsługa jest prosta:
* wartość korzenia to T[0]
* lewe poddrzewo to T[1] o ile len(T) > 1
* prawe poddrzewo to T[2] o ile len(T) > 2
"""
def makeTree(n, left, right):
return [n, left, right]
"""
Zwraca listę liści drzewa tree
"""
def leaves(tree):
if len(tree) == 1: return tree
return leaves(tree[1]) + ([] if len(tree) == 2 else leaves(tree[2]))
"""
Tworzy drzewo produktowe [z ograniczeniem bound] z listy liczb numbers
"""
def productTree(numbers, bound = None):
n = len(numbers)
if n == 1: return numbers
T = productTree(numbers[0:int(n/2)], bound)
U = productTree(numbers[int(n/2):n], bound)
return makeTree(T[0]*U[0] if not bound else min(bound, T[0]*U[0]), T, U)
"""
Tworzy drzewo reszt modulo x z drzewa produktowego product_tree
"""
def remainderTree(x, product_tree, parent = None):
product_tree[0] = \
x if product_tree[0] == x else \
x % (product_tree[0] if not parent else parent)
if len(product_tree) > 1: remainderTree(x, product_tree[1], parent)
if len(product_tree) > 2: remainderTree(x, product_tree[2], parent)
"""
Testuje liczby z listy numbers pod kątem B-gładkości
"""
from math import log, ceil
def BernsteinBatchSmoothnessTest(numbers, B):
primes_product_tree = productTree(list(primes(B)))
P = primes_product_tree[0]
e = ceil(log(log(max(numbers), 2), 2))
numbers_product_tree = productTree(numbers)
remainderTree(P, numbers_product_tree)
remainders = leaves(numbers_product_tree)
assert len(numbers) == len(remainders)
smooth = []
for i in range(len(numbers)):
x, r = numbers[i], remainders[i]
if gcd(int(r**(2**e) % x), x) == x:
smooth.append(x)
return smooth
B = 50
K = 100
P_list = list(primes(B*2, B*4))
p1 = P_list[randrange(0, len(P_list))]
p2 = P_list[randrange(0, len(P_list))]
assert p1 != p2
n = p1*p2
print "Następujące liczby są %s-gładkie: \n" % B
print " -według zmodyfikowanego sita Eratostenesa: \n\n %s\n" % Sieving_B_smooth_in_x2_n(B, n, K)
print " -według testu Bernsteina: \n\n %s\n" % \
BernsteinBatchSmoothnessTest([x*x - n for x in range(int(ceil(sqrt(n))), int(ceil(sqrt(n))) + K)], B)
Następujące liczby są 50-gładkie:
-według zmodyfikowanego sita Eratostenesa:
[152, 435, 720, 1296, 1587, 1880, 2175, 3072, 3375, 3680, 4920, 6192, 6840, 8160, 8832, 9512, 10200, 11600, 12312, 13395, 13760, 19475, 19872, 20672, 22707, 23120, 25215, 26496, 27360, 30000, 31347, 32712, 34560, 35496, 37392]
-według testu Bernsteina:
[152, 435, 720, 1296, 1587, 1880, 2175, 3072, 3375, 3680, 4920, 6192, 6840, 8160, 8832, 9512, 10200, 11600, 12312, 13395, 13760, 19475, 19872, 20672, 22707, 23120, 25215, 26496, 27360, 30000, 31347, 32712, 34560, 35496, 37392]
