Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Taxicab count cubes

By: ploffie on Nov 12th, 2012  |  syntax: Python  |  size: 3.52 KB  |  hits: 45  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  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
create a new version of this paste RAW Paste Data