fork download
copy 
  1. from math import inf, ceil, gcd
  2. from time import clock
  3. from itertools import groupby, accumulate
  4. from operator import mul
  5. from random import randrange
  6.  
  7. def isqrt(n):
  8.     if n < 0:
  9.         raise ValueError("x must be non-negative")
  10.     if n == 0:
  11.         return 0
  12.     x = 2 ** ((n.bit_length() - 1) // 2 + 1)    # x > sqrt(n)
  13.     y = x - 1
  14.     while y < x:
  15.         x = y
  16.         y = (x + n//x) // 2
  17.     return x
  18.  
  19. def iroot(k, n):
  20.     'ínt kth root - uses Newton'
  21.     if k == 2:
  22.         return isqrt(n)
  23.     if n < 0:
  24.         if k % 2:
  25.             return -iroot(k, -n)
  26.         else:
  27.             raise ValueError("n must be non-negative")
  28.     if k < 0:
  29.         raise ValueError("k must be non-negative")
  30.     if n in (0, 1):
  31.         return n
  32.     km1 = k-1
  33.     s = 2 ** ((n.bit_length() - 1)//k + 1)
  34.     u = s - 1
  35.     while u < s:
  36.         s = u
  37.         u = (km1*s + n // s**km1) // k
  38.     return s
  39.  
  40. _small_primes = set((2, 3, 5, 7, 11, 13, 17, 19, 23))
  41.  
  42. def is_prime(n, k=25):
  43.     """ Miller-Rabin
  44.         returns: True if n is probably prime
  45.                     False if n is composite
  46.         k: number of random trial used
  47.     """
  48.     n = abs(n)
  49.     if n < 29:
  50.         return n in _small_primes
  51.     for pi in _small_primes:
  52.         if not n % pi:
  53.             return False
  54.     if pow(2, n-1, n) != 1:
  55.         return False
  56.     d, s = n - 1, 0
  57.     while not d % 2:
  58.         d //= 2
  59.         s += 1
  60.  
  61.     def is_spsp(n, a):
  62.         t = pow(a, d, n)
  63.         if t in (1, n-1):
  64.             return True
  65.         for _ in range(s-1):
  66.             t = (t * t) % n
  67.             if t == n - 1:
  68.                 return True
  69.             if t == 1:
  70.                 return False
  71.         return False
  72.     return all(is_spsp(n, randrange(2, n-1)) for _ in range(k))
  73.  
  74. def pollard_rho(n, c=3):
  75.     t, h, d = 2, 2, 1
  76.     while d == 1:
  77.         t = (t * t + c) % n
  78.         h = (h * h + c) % n
  79.         h = (h * h + c) % n
  80.         d = gcd(t - h, n)
  81.     if d == n:  # cycle detected
  82.         return pollard_rho(n, c + 1)
  83.     return d
  84.  
  85. def rho_factors(n):
  86.     if n < 2:
  87.         return [-1] + rho_factors(-n) if n < 0 else []
  88.     fs = []
  89.     while not n % 2:
  90.         n //= 2
  91.         fs.append(2)
  92.     while not (is_prime(n) or n == 1):
  93.         f = pollard_rho(n)
  94.         flist = [f] if is_prime(f) else rho_factors(f)
  95.         for f in flist:
  96.             while not n % f:
  97.                 n //= f
  98.                 fs.append(f)
  99.     if n > 1:
  100.         fs.append(n)
  101.     return sorted(fs)
  102.  
  103. def divisors_from_factors(factors):
  104.     """FAST - list of all divisors from factors"""
  105.     fact = [(1,) + tuple(accumulate(g, mul)) for _, g in groupby(factors)]
  106.     div = [1]
  107.     for g in fact:
  108.         div = [d * f for d in div for f in g]
  109.     return div
  110.  
  111. def divisors(n, factor_gen=rho_factors):
  112.     """FAST - list of all divisors of n"""
  113.     return divisors_from_factors(factor_gen(n))
  114.  
  115. def bf(n):
  116.     """i <= j <= p"""
  117.     tmin = inf
  118.     fmin = None
  119.     n3 = ceil(n ** (1/3)) + 1
  120.     for i in range(1, n3):
  121.         if n % i == 0:
  122.             m = n // i
  123.             m2 = ceil(m ** (1/2)) + 1
  124.             for j in range(i, m2):
  125.                 if m % j == 0:
  126.                     p = m // j
  127.                     s = i + j + p
  128.                     if s < tmin:
  129.                         tmin = s
  130.                         fmin = i, j, p
  131.     return fmin, tmin
  132.  
  133. def threeway(n):
  134.     divs = sorted(divisors(n), reverse=True)
  135.     n3 = iroot(3, n)
  136.     min_val, min_sol = inf, None
  137.     for x in divs:
  138.         if x > n3:
  139.             continue
  140.         m = n // x
  141.         sqrtm = isqrt(m)
  142.         if x + sqrtm + m // sqrtm >= min_val:
  143.             break
  144.         for y in divs:
  145.             if y < x:
  146.                 break
  147.             if y > sqrtm:
  148.                 continue
  149.             if m % y == 0:
  150.                 z = m // y
  151.                 if x + y + z < min_val:
  152.                     min_sol, min_val = (x, y, z), x + y + z
  153.                 break
  154.     return min_sol, min_val
  155.  
  156. N = 1890
  157. print(N)
  158. print(threeway(N))
  159. M = 10**40
  160. N = randrange(M, 2*M)
  161. print(N)
  162. print(threeway(N))
  163.  
  164. if False:
  165.     for n in range(100000+1, 200000+1):
  166.         if n % 1000 == 0:
  167.             print("n=", n)
  168.         f1, s1 = bf(n)
  169.         f3, s3 = threeway(n)
  170.         assert s1 == s3, "{} {} {} {} {} {}".format(n, s1, s3, facs, f1, f3)
  171.  
  172.     print("Test OK")
  173.  
  174. if False:
  175.     N1 = 10**12
  176.     N2 = N1 + 1
  177.  
  178.     t0 = clock()
  179.     for n in range(N1, N2):
  180.         f1, s1 = threeway(n)
  181.     print(clock() - t0)
  182.  
  183.     t0 = clock()
  184.     for n in range(N1, N2):
  185.         f1, s1 = bf(n)
  186.     print(clock() - t0)
  187.  
Time limit exceeded #stdin #stdout 15s 12064KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/58J5b0
Three way minimum sum
language:
Python 3 nbc (python 3.7.3)
created:
4 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!