# Taxicab count for large values

Nov 12th, 2012
72
0
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
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