fork download
copy 
  1. import random
  2. from Queue import Queue
  3.  
  4. class Pow:
  5.     def __init__(self, base, exp):
  6.         self.base = base
  7.         self.exp = exp
  8.  
  9. class Pair:
  10.     def __init__(self, x, y):
  11.         self.x = x
  12.         self.y = y
  13.  
  14. def gcd(a, b):
  15.     while b: a, b = b, a % b
  16.     return a
  17.  
  18. def rabin_miller(p):
  19. 	if p < 2: return False
  20. 	if p != 2 and p % 2 == 0: return False
  21. 	s = p - 1
  22. 	while s % 2 == 0: s >>= 1
  23. 	for i in xrange(10):
  24. 		a = random.randrange(p - 1) + 1
  25. 		temp = s
  26. 		mod = pow(a, temp, p)
  27. 		while temp != p-1 and mod != 1 and mod != p-1:
  28. 			mod = (mod * mod) % p
  29. 			temp = temp * 2
  30. 		if mod != p-1 and temp % 2 == 0:
  31. 			return False
  32. 	return True
  33.  
  34. def brent(n):
  35.     if n % 2 == 0: return 2;
  36.     x, c, m = random.randrange(0,n), random.randrange(1,n), random.randrange(1,n)
  37.     y, r, q = x, 1, 1
  38.     g, ys = 0, 0
  39.     while True:
  40.         x = y
  41.         for i in range(r):
  42.             y, k = (y*y+c)%n, 0
  43.         while True:
  44.             ys = y
  45.             for i in range(min(m,r-k)):
  46.                 y, q = (y*y+c)%n, q*abs(x-y)%n
  47.             g, k = gcd(q,n), k+m
  48.             if k >= r or g > 1: break
  49.         r = 2 * r
  50.         if g > 1: break
  51.     if g == n:
  52.         while True :
  53.             ys, g = (x*x+c)%n, gcd(abs(x-ys),n)
  54.             if g > 1: break
  55.     return g
  56.  
  57. def pollard(n):
  58.     if n % 2 == 0:
  59.         return 2;
  60.     x = random.randrange(2,1000000)
  61.     c = random.randrange(2,1000000)
  62.     y = x
  63.     d = 1
  64.     while d == 1:
  65.         x = (x * x + c) % n
  66.         y = (y * y + c) % n
  67.         y = (y * y + c) % n
  68.         d = gcd(x - y, n)
  69.         if d == n:
  70.             break
  71.     return d
  72.  
  73. #    Prime Factorization by Pollard Rho alrorithm
  74. #    http://c...content-available-to-author-only...e.com/recipes/577037-pollard-rho-prime-factorization/
  75. def factor(n):
  76.     q1 = Queue()
  77.     q2 = []
  78.     q1.put(n)
  79.     while not q1.empty():
  80.         l = q1.get()
  81.         if rabin_miller(l):
  82.             q2.append(l)
  83.             continue
  84.         d = pollard(l)
  85.         if d == l: q1.put(l)
  86.         else:
  87.             q1.put(d)
  88.             q1.put(l / d)
  89.     return q2
  90.  
  91. #   Square root of n in integer.
  92. def sqrt(n):
  93.     if n <= 0: return 0
  94.     s = 1
  95.     t = n;
  96.     while s < t:
  97.         s <<= 1
  98.         t >>= 1
  99.     while True:
  100.         t = s
  101.         s = (n / s + s) >> 1
  102.         if s >= t: break
  103.     return t
  104.  
  105. #    Return (a|n).
  106. #    (a|n) = 0  if a = 0 (mod p)
  107. #          = +1 if x^2 = a (mod p) for some x
  108. #          = -1 if there is no such x
  109. def jacobi(a, n):
  110.     t = 1
  111.     while a != 0:
  112.         while a % 2 == 0:
  113.             a = a / 2
  114.             if n % 8 == 3 or n % 8 == 5: t = -t
  115.         (a, n) = (n, a)
  116.         if a % 4 == 3 and n % 4 == 3: t = -t
  117.         a = a % n
  118.     if n == 1: return t
  119.     return 0
  120.  
  121. #    Return a^n (mod p)
  122. def pow_mod(a, n, p): 
  123.     if n == 0: return 1
  124.     if n % 2 == 0: return pow_mod(a * a % p, n / 2, p)
  125.     return a * pow_mod( a, n - 1, p ) % p
  126.  
  127. #    Find x who satisfies x*x=n (mod prime p) by Shanks-Tonelli algorithm.
  128. #    http://w...content-available-to-author-only...x.com/wiki/Shanks-Tonelli_algorithm
  129. def mod_sqrt_by_shanks_tonelli(prime, arg):
  130.     y = 2
  131.     result = 0
  132.     r = 0
  133.     q = prime - 1
  134.     if jacobi(arg, prime) == -1: return -1
  135.     while jacobi(y, prime) != -1: y += 1
  136.     result = prime - 1
  137.     while q % 2 == 0:
  138.         r += 1
  139.         q /= 2
  140.     result >>= r
  141.     y = pow_mod(y, result, prime)
  142.     result >>= 1
  143.     b = pow_mod(arg, result, prime)
  144.     result = arg * b
  145.     result %= prime
  146.     b *= result
  147.     b %= prime
  148.     while b != 1:   
  149.         t = b * b
  150.         t %= prime
  151.         m = 1
  152.         while t != 1:
  153.             t *= t
  154.             t %= prime
  155.             m += 1
  156.         t = 0
  157.         t |= (1 << (r-m-1))
  158.         t %= prime
  159.         y %= prime
  160.         t = pow_mod(y, t, prime)
  161.         y = t * t
  162.         r = m
  163.         result *= t
  164.         result %= prime
  165.         b *= y % prime
  166.         b %= prime
  167.     return result
  168.  
  169. #    Express prime p as a sum of squared two integers, based on Serret's algorithm.
  170. #    See p.122 of "The higher arithmetic: an introduction to the theory of numbers".
  171. def square_sum_of_prime_by_serret(p):
  172.     if p == 2: return Pair(1, 1)
  173.     if p % 4 != 1: return None
  174.     q = Pair(p, mod_sqrt_by_shanks_tonelli(p, p - 1))
  175.     if 2 * q.y >= p: q.y = p - q.y
  176.     v = []
  177.     while q.y > 0: 
  178.         v.append(q.x / q.y)
  179.         q.x, q.y = q.y, q.x - (q.x / q.y) * q.y
  180.     q.x = v[len(v) / 2]
  181.     q.y = 1
  182.     for i in range(len(v)/2+1, len(v)):
  183.         q.x, q.y = v[i] * q.x + q.y, q.x
  184.     if 2 * q.x >= p:
  185.         q.x = p - q.x
  186.     q.y = sqrt(p - q.x * q.x)
  187.     if q.y < q.x:
  188.         q.x, q.y = q.y, q.x
  189.     return q
  190.  
  191. def expand(a, b, c, d, u):
  192.     q = Pair(0, a * d + b * c)
  193.     if a * c >= b * d:
  194.         q.x = a * c - b * d
  195.     else:
  196.         q.x = b * d - a * c
  197.     if q.y < q.x:
  198.         q.x, q.y = q.y, q.x
  199.     for p in u:
  200.         if p.x == q.x: return False
  201.     u.append(q)
  202.     return u
  203.  
  204. #    Find all pairs x>0 and y>0 which satisfy x^2+y^2=n (x<y).
  205. #    n is the form of exponent expression.
  206. def square_sum_of_integer(v):
  207.     n = 0
  208.     u, v1, v2 = [], [], []
  209.     for p in v:
  210.         if p.base == 2 or p.base % 4 == 1:
  211.             for i in range(p.exp): v1.append(p.base)
  212.             n += p.exp
  213.         else:
  214.             if p.exp % 2 != 0: return False
  215.             for i in range(p.exp / 2): v2.append(p.base)
  216.     if len(v1) == 0:
  217.         q = Pair(1, 0)
  218.         for p in v:
  219.             for i in range(p.exp): q.x *= p.base
  220.         u.append(q)
  221.         return false
  222.     w = []
  223.     r = square_sum_of_prime_by_serret(v1[0])
  224.     if r == None:
  225.         r.x = r.y = v1[0]
  226.     u.append(r)
  227.     for i in range(1, len(v1)):
  228.         r = square_sum_of_prime_by_serret(v1[i])
  229.         if r != None:
  230.             for q in u:
  231.                 w = expand(q.x, q.y, r.x, r.y, w)
  232.                 w = expand(q.x, q.y, r.y, r.x, w)
  233.         u, w = w, u
  234.         w = []
  235.     for i in range(1, len(v2)):
  236.         for q in u:
  237.             q.x *= v2[i]
  238.             q.y *= v2[i]
  239.     return u
  240.  
  241. def solve(n):
  242.     r = 0
  243.     t = 10**(n/2)
  244.     v = []
  245.     u = []
  246.     m = t*t + 1
  247.  
  248.     p = factor(m)
  249.     p.sort()
  250.     q = Pow(p[0], 1)
  251.     for i in range(1, len(p)):
  252.         if p[i] != p[i-1]:
  253.             v.append(q)
  254.             q = Pow(p[i], 1)
  255.         else: q.exp += 1
  256.     v.append(q)
  257.  
  258.     u = square_sum_of_integer(v)
  259.  
  260.     for s in u:
  261.         if s.x % 2 != 0: s.x, s.y = s.y, s.x
  262.         x = t / 2 + s.x / 2
  263.         y = (1 + s.y) / 2
  264.         if x >= t / 10 and x < t and y < t:
  265. #            print "  " + str(t * x + y)
  266.             r += t * x + y
  267.         x = t / 2 - s.x / 2
  268.         if x >= t / 10 and x < t and y < t:
  269. #            print "  " + str(t * x + y)
  270.             r += t * x + y
  271.     print str(n) + ': ' + str(r)
  272.  
  273. for n in range(2, 38, 2):
  274. 	solve(n)
  275. solve(68)
Success #stdin #stdout 2.16s 11840KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
2: 0
4: 10066
6: 990100
8: 94122353
10: 29901173703
12: 2999901309872
14: 681781421599205
16: 147679514888352862
18: 2979998020199999802
20: 699990001048290009817
22: 69614781398394320917605
24: 2882399981176470599995296
26: 299009900990116653047262503
28: 69878645728081601689809338163
30: 62261981581052976186313016462048
32: 699200377406218281243977638802438
34: 66942179745289104121114968938275755
36: 30698564405147901939971608356414351820
68: 694954127022443331439683664029508243184680152400711799691388780508888
https://ideone.com/V56Fbr
「カット・アンド・スクエア」問題 Kawazoeのコード
language:
Python (cpython 2.7.16)
created:
8 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!