ploffie

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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×