-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy patheuler343.py
More file actions
54 lines (49 loc) · 1.48 KB
/
euler343.py
File metadata and controls
54 lines (49 loc) · 1.48 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
from eulertools import primes
def is_square(n, p):
return pow(n, (p-1) // 2, p) == 1
def exp_in_fp2(l, exponent, k, p):
if exponent == 0:
return (1, 0)
elif exponent == 1:
return l
elif exponent % 2 == 0:
x, y = l
return exp_in_fp2(((x**2+y**2*k) % p , (2*x*y) % p), exponent//2, k, p)
else:
x, y = l
z, w = exp_in_fp2(((x**2+y**2*k) % p , (2*x*y) % p), (exponent-1)//2, k, p)
return ((x*z+y*w*k) % p, (x*w+y*z) % p)
def find_square_roots(n, p):
"""Implementing Cipolli's algorithm"""
a = 1
while is_square((a**2-n) % p, p):
a += 1
result = exp_in_fp2((a, 1), (p+1)//2, (a**2-n) % p, p)
return result[0], p - result[0]
def main(n):
prime_list = primes(n+2)
l = [(1, k**3+1) for k in xrange(n+1)]
for p in prime_list:
remainders = set([])
if p == 2:
remainders.add(1)
elif p == 3:
remainders.add(2)
else:
if is_square(p-3, p):
a, b = find_square_roots(p-3, p)
remainders.add((1+a)*((p+1)//2) % p)
remainders.add((1+b)*((p+1)//2) % p)
remainders.add(p-1)
for rem in remainders:
for i in xrange(rem, n+1, p):
a, b = l[i]
while b % p == 0:
b //= p
l[i] = (p, b)
result = 0
for i in l[1:]:
result += max(i) - 1
return result
print main(100)
print main(2*10**6)