 # Taxicab count cubes

Nov 12th, 2012
96
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import itertools
2. import bisect
3.
4. cubs = [i**3 for i in xrange(1000000)] # needed for cubes_bisect
5.
6. def invpow(x, n):
7.     """ x , n > 0 integers
8.        return y, such that y ** n <= x and (y+1) ** n > x
9.    """
10.     return int((x + 0.5) ** (1.0/3))
11.
12.
13. def cubes_bisect(n):
14.     """ using bisection to find x,y pairs - a little slower than cubes
15.    """
16.     y = invpow(n, 3)
17.     x = 1
18.     while x <= y:
19.         y3 = cubs[y]
20.         x = bisect.bisect_left(cubs, n - y3, x, y + 1)
21.         if cubs[x] + y3 == n:
22.             yield x, y
23.         x += 1
24.         y -= 1
25.
26. def cubes(n):
27.     """generator for all sums of 3rd powers to get n
28.       start with pair x, y, where y is the highest number, such that y**3<=n
29.       If n is odd, (x,y) can only be an odd-even pair. If n is even, (x,y) can
30.       only be odd-odd or even-even. While x^3+y^3 is smaller than n, x can
31.       advance with 2 steps i.s.o one step, as a one step increase cannot be a
32.       solution. If y is decreased x must increase to make sure that the sum can
33.       be a solution (is odd or even, like n).
34.       This version is about a factor 2 faster than cubes0
35.    """
36.     y = invpow(n, 3)
37.     y3 = y ** 3
38.     x = 0 + (n - y) & 1 # make x odd if n - y^3 is odd
39.     while x <= y:
40.         test = x ** 3 + y3
41.         if test < n:
42.             x += 2
43.         else:
44.             if test == n:
45.                 yield x, y
46.             x += 1
47.             y -= 1
48.             y3 = y ** 3
49.
50. def check_cubes(n, ncubes):
51.     """check if number of cubes is ncubes for number n"""
52.     cub = list(itertools.islice(cubes(n), ncubes+1))
53.     return cub if len(cub) == ncubes else False
54.
55. def list_numbers(maxnumber, ncubes):
56.     """ list all numbers with ncubes ways of writing it as the sum of 2 cubes
57.    """
58.     for i in xrange(2, maxnumber):
59.         res = check_cubes(i, ncubes)
60.         if res:
61.             print i, res
62. import time
63. t0 = time.clock()
64. list_numbers(100000, 2)
65. print time.clock() - t0
66. """
67. 1729 [(1, 12), (9, 10)]
68. 4104 [(2, 16), (9, 15)]
69. 13832 [(2, 24), (18, 20)]
70. 20683 [(10, 27), (19, 24)]
71. 32832 [(4, 32), (18, 30)]
72. 39312 [(2, 34), (15, 33)]
73. 40033 [(9, 34), (16, 33)]
74. 46683 [(3, 36), (27, 30)]
75. 64232 [(17, 39), (26, 36)]
76. 65728 [(12, 40), (31, 33)]
77. """
78.
79. import sys
80. from time import clock
81. t0 = clock()
82. def taxicab(ncubes):
83.     """return lowest number with ncubes cubes"""
84.     n = 0
85.     while 1:
86.         n += 1
87.         if not n % 1000000:
88.             print n, clock() - t0
89.             sys.stdout.flush()
90.         res = check_cubes(n, ncubes)
91.         if res:
92.             return n, res
93.
94. #print taxicab(3) # this takes 5268 sec -> 87539319
95. #print clock() - t0
96.
97. def test_invpow():
98.     """test for function invpow"""
99.     i = 2
100.     while i < 100000000:
101.         a = invpow(i, 3)
102.         assert a ** 3 <= i
103.         assert (a+1)**3 > i
104.         i += 1234
105.
106. def test_cubes():
107.     """test if functions cubes and cubes0 give same solution"""
108.     i = 2
109.     while i < 1000000:
110.         for e1,e2 in itertools.izip(cubes(i), cubes0(i)):
111.             assert e1 == e2
112.         i += 1
113.
114. def cubes0(n):
115.     """generator for all sums of 3rd powers to get n
116.       cubes is a faster version
117.    """
118.     x = 0
119.     y = invpow(n, 3)
120.     while x <= y:
121.         x3 = x ** 3
122.         y3 = y ** 3
123.         if x3 + y3 > n:
124.             y -= 1
125.         elif x3 + y3 < n:
126.             x += 1
127.         else:
128.             yield x, y
129.             x += 1
130.             y -= 1
RAW Paste Data