fork(5) download
copy 
  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
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[3, 3, 5, 911, 1051, 48179, 48179]
https://ideone.com/hTmzjd
language:
Python (cpython 2.7.16)
created:
2 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!