Refining List Comprehension in Haskell

Question

So, I'm trying to create a list of taxi numbers in Haskell. Taxi numbers are numbers that can be written as the sum of two different cubes in two different ways: the smallest number is 1729 = 1^3 + 12^3 = 9^3 + 10^3 .

Currently, I'm just worried about creating four numbers that "make up" a taxi number, for example. (1,12,9,10), and they were asked to use an understanding of the list (which I am not very familiar with). This function will generate all 4-tuples, where the largest number is at most n:

taxi n = [(a,b,c,d) | a <- [1..n], b <- [1..n], c <- [1..n], d <- [1..n], a^3 + b^3 == c^3 + d^3, a < b, a < c, c < d]

However, there may be several reasons:

• The domain for a, b, c, d is identical, but I can’t decide how to simplify this code, so [1..n] written only once.
• I want an endless list with no upper limit, so I can just leave the program running as long as possible and stop it when I want. It is clear that if I simply set a <- [1..] etc., then the program will never finish evaluating anything.
• The program is very slow: only taxi 50 takes 19 seconds.

Any speed optimization would also be good, but if not the most naive method that I use is enough.

Answer 1

Your limitations imply a < c < d < b . So, let b start from the outside and let the rest work in the corresponding lower ranges:

 taxi n = [ (a,b,c,d) | b <- [1..n], d <- [1..b-1], c <- [1..d-1], a <- [1..c-1], a^3 + b^3 == c^3 + d^3 ] 

To go endlessly, just use b <- [1..] .

---

Another big improvement is to calculate one of four variables from the other three:

 taxi = [ (a,b,c,d) | b <- [1..], c <- [1..b-1], a <- [1..c-1], let d3 = a^3 + b^3 - c^3, let d = round(fromIntegral(d3)**(1/3)), c < d, d^3 == d3 ] 

---

Benchmarking taxi 50 in GHCi with :set +s , as you did:

 Yours: (16.49 secs, 17,672,510,704 bytes) My first: (0.65 secs, 658,537,184 bytes) My second: (0.09 secs, 66,229,376 bytes) (modified to use b <- [1..n] again) Daniel first: (1.94 secs, 2,016,810,312 bytes) Daniel second: (2.87 secs, 2,434,309,440 bytes) 

Answer 2

The domain for a, b, c, d is identical, but I can’t figure out how to simplify this code, so [1..n] written only once.

Use [a,b,c,d] <- replicateM 4 [1..n] .

The program is very slow: only a taxi 50 takes 19 seconds.

One of the cheap improvements is to bake the conditions a<b , a<c and c<d into understanding.

 taxi n = [(a,b,c,d) | a <- [1..n], b <- [a+1..n], c <- [a+1..n], d <- [c+1..n], a^3 + b^3 == c^3 + d^3] 

This greatly speeds up the work on my machine.

Or, for better compatibility with the next (and previous) part of the answer, relate b , c and d to offsets.

 taxi n = [ (a,b,c,d) | a <- [1..n] , b_ <- [1..n], let b = a+b_, b<=n , c_ <- [1..n], let c = a+c_, c<=n , d_ <- [1..n], let d = c+d_, d<=n , a^3 + b^3 == c^3 + d^3 ] 

I want an endless list with no upper limit.

See my answer to the Cartesian product of 2 lists in Haskell for a hint. tl; dr use choices .

Answer 3

Raising Stefan is his excellent answer. Given a^3 + b^3 == c^3 + d^3 , we only need to look at integers for which 0 < a < c < b . We introduce this (infinite) iterative structure

 -- generates all integers x, y and z for which holds 0 < x < y < z triplets = [(x, y, z) | z <- [3 .. ], y <- [2 .. z - 1], x <- [1 .. y - 1]] 

This will give us triplets that are easily accessible from our understanding of the list, which we will present here after. For people with a Python background, this should be equivalent to Python yield .

 1 2 3 1 2 4 1 3 4 2 3 4 1 2 5 1 3 5 2 3 5 1 4 5 2 4 5 3 4 5 1 2 6 1 3 6 2 3 6 1 4 6 2 4 6 3 4 6 1 5 6 2 5 6 3 5 6 4 5 6 

Next, we need something (quickly) to find the largest cube and check that the integers are a cube, also called an integer cube root. There is this package Math.NumberTheory.Powers.Cubes , which has functions for these tasks. Or just use these

 -- given integer x >= 0 find the largest integer r such that r^3 <= x largestCube :: Integral a => a -> a largestCube x = let powers_of_two = iterate ((*) 2) 1 upper = head [j | j <- powers_of_two, x < j ^ 3] in largestCubeSub 0 upper x largestCubeSub :: Integral a => a -> a -> a -> a largestCubeSub lower upper x | lower + 1 == upper = lower | b ^ 3 <= x = largestCubeSub b upper x | otherwise = largestCubeSub lower bx where b = div (lower + upper) 2 -- test if an integer x >= 0 is a cube isCube :: Integral a => a -> Bool isCube x = (largestCube x) ^ 3 == x 

Now your compact list comprehension for the first 50 taxi numbers looks like

 *Main> condition = \abc -> and [isCube (a^3 + b^3 - c^3), a^3 + b^3 - c^3 > c^3] *Main> taxi = [(a, b, c, largestCube (a^3 + b^3 - c^3)) | (a, c, b) <- triplets, condition abc] *Main> first50 = take 50 taxi 

Print them using

 *Main> single_line = \(x, y, z, u) -> unwords [show i | i <- [x, y, z, u]] *Main> putStrLn $ unlines $ map single_line first50 

Gives

 1 12 9 10 2 16 9 15 2 24 18 20 10 27 19 24 4 32 18 30 2 34 15 33 9 34 16 33 3 36 27 30 17 39 26 36 12 40 31 33 6 48 27 45 4 48 36 40 12 51 38 43 8 53 29 50 20 54 38 48 17 55 24 54 9 58 22 57 3 60 22 59 5 60 45 50 8 64 36 60 30 67 51 58 4 68 30 66 18 68 32 66 42 69 56 61 6 72 54 60 17 76 38 73 5 76 48 69 34 78 52 72 10 80 45 75 15 80 54 71 24 80 62 66 30 81 57 72 51 82 64 75 7 84 63 70 2 89 41 86 11 93 30 92 23 94 63 84 12 96 54 90 50 96 59 93 8 96 72 80 20 97 33 96 47 97 66 90 35 98 59 92 24 98 63 89 29 99 60 92 6 102 45 99 27 102 48 99 23 102 60 95 24 102 76 86 1 103 64 94 

Within a few seconds, he returns the first 50 taxi numbers.