fork(4) download
  1. # 100011322432435545
  2.  
  3. def isqrt(n): # newton
  4. x = n; y = (x + 1) // 2
  5. while y < x:
  6. x = y; y = (x + n // x) // 2
  7. return x
  8.  
  9. def isSquare(n):
  10. # def q(n):
  11. # from sets import Set
  12. # s, sum = Set(), 0
  13. # for x in xrange(0,n):
  14. # t = pow(x,2,n)
  15. # if t not in s:
  16. # s.add(t)
  17. # sum += pow(2,t)
  18. # return sum
  19. # q(32) => 33751571
  20. # q(27) => 38348435
  21. # q(25) => 19483219
  22. # q(19) => 199411
  23. # q(17) => 107287
  24. # q(13) => 5659
  25. # q(11) => 571
  26. # q(7) => 23
  27. # 99.82% of non-squares caught by bloom
  28. # filters before square root calculation
  29. if 33751571>>(n%32)&1==0: return False
  30. if 38348435>>(n%27)&1==0: return False
  31. if 19483219>>(n%25)&1==0: return False
  32. if 199411>>(n%19)&1==0: return False
  33. if 107287>>(n%17)&1==0: return False
  34. if 5659>>(n%13)&1==0: return False
  35. if 571>>(n%11)&1==0: return False
  36. if 23>>(n% 7)&1==0: return False
  37. s = isqrt(n)
  38. if s * s == n: return s
  39. return False
  40.  
  41. def isPrime(n, k=5): # miller-rabin
  42. from random import randint
  43. ps = [2,3,5,7,11,13,17,19,23,29]
  44. if n < 2: return False
  45. for p in ps:
  46. if n%p == 0: return n==p
  47. s, d = 0, n-1
  48. while d%2 == 0:
  49. s, d = s+1, d/2
  50. for i in xrange(k):
  51. x = pow(randint(2, n-1), d, n)
  52. if x == 1 or x == n-1: continue
  53. for r in xrange(1, s):
  54. x = (x * x) % n
  55. if x == 1: return False
  56. if x == n-1: break
  57. else: return False
  58. return True
  59.  
  60. def factors(n):
  61. def gcd(a, b): # euclid
  62. if b == 0: return a
  63. return gcd(b, a%b)
  64. # 2,3,5-wheel to cube root of n
  65. wheel = [1,2,2,4,2,4,2,4,6,2,6]
  66. w, f, fs = 0, 2, []
  67. while f*f*f <= n:
  68. # print f, n, fs
  69. while n % f == 0:
  70. fs.append(f); n /= f
  71. if n == 1: return fs
  72. if isPrime(n):
  73. fs.append(n)
  74. return fs
  75. f, w = f+wheel[w], w+1
  76. if w == 11: w = 3
  77. # n must be semi-prime, so use
  78. # william hart's fermat variant
  79. if isSquare(n):
  80. # print "square", n
  81. f = isqrt(n);
  82. fs.append(f); fs.append(f)
  83. return fs
  84. for i in xrange(1,n):
  85. s = isqrt(n*i)
  86. if s*s <= n*i: s += 1
  87. m = pow(s,2,n)
  88. # print i, s, m
  89. if isSquare(m):
  90. t = isqrt(m); f = gcd(n, s-t)
  91. fs.append(f); fs.append(n/f)
  92. return sorted(fs)
  93.  
  94. print factors(100011322432435545)
Success #stdin #stdout 0.02s 9840KB
stdin
Standard input is empty
stdout
[3, 3, 5, 911, 1051, 48179, 48179]