Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import itertools
- import bisect
- cubs = [i**3 for i in xrange(1000000)] # needed for cubes_bisect
- def invpow(x, n):
- """ x , n > 0 integers
- return y, such that y ** n <= x and (y+1) ** n > x
- """
- return int((x + 0.5) ** (1.0/3))
- def cubes_bisect(n):
- """ using bisection to find x,y pairs - a little slower than cubes
- """
- y = invpow(n, 3)
- x = 1
- while x <= y:
- y3 = cubs[y]
- x = bisect.bisect_left(cubs, n - y3, x, y + 1)
- if cubs[x] + y3 == n:
- yield x, y
- x += 1
- y -= 1
- def cubes(n):
- """generator for all sums of 3rd powers to get n
- start with pair x, y, where y is the highest number, such that y**3<=n
- If n is odd, (x,y) can only be an odd-even pair. If n is even, (x,y) can
- only be odd-odd or even-even. While x^3+y^3 is smaller than n, x can
- advance with 2 steps i.s.o one step, as a one step increase cannot be a
- solution. If y is decreased x must increase to make sure that the sum can
- be a solution (is odd or even, like n).
- This version is about a factor 2 faster than cubes0
- """
- y = invpow(n, 3)
- y3 = y ** 3
- x = 0 + (n - y) & 1 # make x odd if n - y^3 is odd
- while x <= y:
- test = x ** 3 + y3
- if test < n:
- x += 2
- else:
- if test == n:
- yield x, y
- x += 1
- y -= 1
- y3 = y ** 3
- def check_cubes(n, ncubes):
- """check if number of cubes is ncubes for number n"""
- cub = list(itertools.islice(cubes(n), ncubes+1))
- return cub if len(cub) == ncubes else False
- def list_numbers(maxnumber, ncubes):
- """ list all numbers with ncubes ways of writing it as the sum of 2 cubes
- """
- for i in xrange(2, maxnumber):
- res = check_cubes(i, ncubes)
- if res:
- print i, res
- import time
- t0 = time.clock()
- list_numbers(100000, 2)
- print time.clock() - t0
- """
- 1729 [(1, 12), (9, 10)]
- 4104 [(2, 16), (9, 15)]
- 13832 [(2, 24), (18, 20)]
- 20683 [(10, 27), (19, 24)]
- 32832 [(4, 32), (18, 30)]
- 39312 [(2, 34), (15, 33)]
- 40033 [(9, 34), (16, 33)]
- 46683 [(3, 36), (27, 30)]
- 64232 [(17, 39), (26, 36)]
- 65728 [(12, 40), (31, 33)]
- """
- import sys
- from time import clock
- t0 = clock()
- def taxicab(ncubes):
- """return lowest number with ncubes cubes"""
- n = 0
- while 1:
- n += 1
- if not n % 1000000:
- print n, clock() - t0
- sys.stdout.flush()
- res = check_cubes(n, ncubes)
- if res:
- return n, res
- #print taxicab(3) # this takes 5268 sec -> 87539319
- #print clock() - t0
- def test_invpow():
- """test for function invpow"""
- i = 2
- while i < 100000000:
- a = invpow(i, 3)
- assert a ** 3 <= i
- assert (a+1)**3 > i
- i += 1234
- def test_cubes():
- """test if functions cubes and cubes0 give same solution"""
- i = 2
- while i < 1000000:
- for e1,e2 in itertools.izip(cubes(i), cubes0(i)):
- assert e1 == e2
- i += 1
- def cubes0(n):
- """generator for all sums of 3rd powers to get n
- cubes is a faster version
- """
- x = 0
- y = invpow(n, 3)
- while x <= y:
- x3 = x ** 3
- y3 = y ** 3
- if x3 + y3 > n:
- y -= 1
- elif x3 + y3 < n:
- x += 1
- else:
- yield x, y
- x += 1
- y -= 1
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement