n = 3850178868941
# The starting point
s0 = 709609764751
# A constant for the step updates
d = 56
# We want to accumulate the product of 5 differences into S, so we need to do five iterations
B = 5
def iterate(s_i):
s_iplus1 = s_i^2 + d
return int(mod(s_iplus1, n))
def iterate_double(s_2i):
s_2_iplus1 = (s_2i^2+d)^2 + d
return int(mod(s_2_iplus1, n))
def initialize_walks(starting_point = s0):
# slow_walk = starting_point
# fast_walk = starting_point
# return slow_walk, fast_walk
# This function returns a set of the initial points for both walks.
# Thus, we return the starting point for both the slow walk and the fast walk.
return starting_point, starting_point
# First, re-initialize the walks to the starting point (i.e. rho_0)
slow_walk, fast_walk = initialize_walks()
s = 1
print("Slow walk: " + str(slow_walk) + ", fast walk: " + str(fast_walk))
for i in range(B):
slow_walk = iterate(slow_walk)
fast_walk = iterate_double(fast_walk)
print("Slow walk: " + str(slow_walk) + ", fast walk: " + str(fast_walk))
s = int(mod(s*(fast_walk-slow_walk), n))
p = gcd(s, n)
p
Slow walk: 709609764751, fast walk: 709609764751 Slow walk: 2106196996297, fast walk: 685190823535 Slow walk: 685190823535, fast walk: 3837605587870 Slow walk: 2688451017762, fast walk: 3528230340708 Slow walk: 3837605587870, fast walk: 1277042349402 Slow walk: 3525830075524, fast walk: 1388605140877
79
s = lcm([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
s
720720
list(range(1,16+1))
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
lcm(range(1,16+1))
720720
s.factor()
2^4 * 3^2 * 5 * 7 * 11 * 13
c = p
nprime = n/c
nprime
48736441379
a = 1364612
# Do not forget the __mod n__!
b = Mod(a ^ s, nprime)
b
8965317773
q1 = int(gcd(b-1, int(nprime)))
q1
1995841
q2 = n / (p*q1)
q2
24419
24419*1995841*79 == n
True
# The naming in the exercise gives this
p = q1
q = q2
(p-1).factor()
2^6 * 3^4 * 5 * 7 * 11
(q-1).factor()
2 * 29 * 421
s.factor()
2^4 * 3^2 * 5 * 7 * 11 * 13
GF(p)(a).multiplicative_order().factor()
2^4 * 3^2 * 5 * 7 * 11
GF(q)(a).multiplicative_order().factor()
29 * 421
ord_a_inFp_divides_s = 0
ord_a_inFp_not_divides_s = 0
for a in range(1, GF(p).cardinality()):
order_of_a = GF(p)(a).multiplicative_order()
if order_of_a.divides(s):
ord_a_inFp_divides_s += 1
else:
ord_a_inFp_not_divides_s += 1
ord_a_inFp_divides_s/(ord_a_inFp_divides_s + ord_a_inFp_not_divides_s)
1/36
ord_a_inFq_divides_s = 0
ord_a_inFq_not_divides_s = 0
for a in range(1, GF(q).cardinality()):
order_of_a = GF(q)(a).multiplicative_order()
if order_of_a.divides(s):
ord_a_inFq_divides_s += 1
else:
ord_a_inFq_not_divides_s += 1
ord_a_inFq_divides_s/(ord_a_inFq_divides_s + ord_a_inFq_not_divides_s)
1/12209