ploffie

Taxicab count for large values

Nov 12th, 2012
48
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import ma.primee as PR
  2. import itertools as IT
  3. import operator
  4. import math
  5. import cubes as CC
  6. from time import clock
  7. combinations = IT.combinations
  8. mul = operator.mul
  9.  
  10. def solveq(a, b, c):
  11.     """ solves ax^2 + bx + c = 0"""
  12.     y = b ** 2  - 4 * a * c
  13.     if y < 0:
  14.         return False
  15.     sqrty = math.sqrt(y)
  16.     return ((-b - sqrty) / (2*a), (-b + sqrty) / (2*a))
  17.  
  18. def check_sol(sol, p):
  19.     """checks if solution of quadratic equation is integer and in (0,p)"""
  20.     res= []
  21.     if sol:
  22.         a = sol[0]
  23.         inta = int(a)
  24.         if abs(inta - a) > 1e-15:
  25.             return 0
  26.         if 0 < inta < p:
  27.             res.append(inta)
  28.         inta = sol[1]
  29.         if 0 < inta < p:
  30.             res.append(inta)
  31.     return res
  32.  
  33. def posint(a, x1, x2):
  34.     """check if a is an integer in (x1,x2)"""
  35.     if not x1 < a < x2:
  36.         return False
  37.     if abs(int(a) - a) < 1e-15:
  38.         return int(a)
  39.    
  40. def combs(n):
  41.     """decompose n in (p,q) with p < q and pq = n
  42.        combs returns all possible values for p
  43.    """
  44.     # search through primefactors
  45.     # make sure that no duplications exist:
  46.     #    go from low to higher factors and keep adding factors until p > sqrt(n)
  47.     sqrtn = math.sqrt(n)
  48.     fac =  [(k, len(list(g))) for k, g in IT.groupby(PR.rho_factors(n))]
  49.     pvalues = []
  50.     Q = [(0, 0, 1)]
  51.     while Q:
  52.         # j: current prime, m: exponent of current prime, p: factor of n
  53.         j, m, p = Q.pop()    
  54.         for i, (fi, ni) in enumerate(fac[j:], j):
  55.             mi = m if i == j else 0
  56.             if mi < ni:
  57.                 pi = p * fi
  58.                 if pi <= sqrtn:
  59.                     pvalues.append(pi)
  60.                     Q.append((i, mi + 1, pi))
  61.     return pvalues
  62.  
  63. def count_cubes(n):
  64.     """method: https://cs.uwaterloo.ca/journals/JIS/wilson10.html#HW54"""
  65.     sol = set()
  66.     for p in combs(n):
  67.         aa = check_sol(solveq(3, -3*p, p**2 - n // p), p)
  68.         if aa:
  69.             for a in aa:
  70.                 b = p - a
  71.                 if a > b:
  72.                     a, b = b, a
  73.                 sol.add((a, p - a))
  74.     return sol
  75.  
  76. if __name__ == '__main__':              
  77.     t0 = clock()    
  78.     print count_cubes(19242016141146624)
  79.     print clock() - t0
  80.     t0 = clock()
  81.     print CC.check_cubes(19242016141146624,4)
  82.     print clock() - t0
  83.    
  84.     t0 = clock()    
  85.     print count_cubes(87539319)
  86.     print clock() - t0
  87.     t0 = clock()
  88.     print CC.check_cubes(87539319 ,3)
  89.     print clock() - t0
RAW Paste Data