Generate a stream of all possible numbers `x^k*y^l` for given `k` and `l` in Haskell

> Generate the function generateExponents k l, which for given k and l generates a stream of all unique possible numbers x^k*y^l in increasing order. For example generateExponents 2 3 = [1,4,8,9,16,25,27...]

For obvious reasons this doesn't work:

generateExponents k l = sort [x^k*y^l | x <- [1..], y <- [1..]]

Then I tried this, which doesn't work either:

generateExponents k l = [n | n <- [1 ..], n `elem` products n]
  where
    xs n = takeWhile (\x -> x ^ k <= n) [1 ..]
    ys n = takeWhile (\y -> y ^ l <= n) [1 ..]
    products n = liftA2 (*) (xs n) (ys n)

What am I doing wrong?

Your algorithm is pretty slow -- it checks every number, and for every number it searches for an appropriate factorization! You can do better by producing an infinite table of answers, and then collapsing the table appropriately. For example, for x^2*y^3, the table looks like:

 x  1    2    3    4    5
y
1   1    4    9   16   25
2   8   32   72  128  200
3  27  108  243  432  675
4  64  256  576 1024 1600
5 125  500 1125 2000 3125

Note two nice features of this table: each row is sorted, and the rows themselves are sorted. This means we can merge them efficiently by simply taking the top-left value, then re-inserting the tail of the first row in its new sorted position. For example, the table above, after emitting 1, would look like:

  4    9   16   25   36
  8   32   72  128  200
 27  108  243  432  675
 64  256  576 1024 1600
125  500 1125 2000 3125

Then, after emitting the top-left value 4:

  8   32   72  128  200
  9   16   25   36   49
 27  108  243  432  675
 64  256  576 1024 1600
125  500 1125 2000 3125

Note how the top row has now become the second row to keep the doubly-sorted property.

This is an efficient way to construct all the right numbers in the right order. Then, the only remaining trick needed is to deduplicate, and for that you can deploy the standard trick map head . group, since duplicates are guaranteed to be next to each other. Here's the full code:

import Data.List

generateExponents' k l = map head . group . go $ [[x^k*y^l | x <- [1..]] | y <- [1..]] where
    go ((x:xs):xss) = x:go (insert xs xss)

It's much, much faster. Compare:

> sum . take 400 $ generateExponents 2 3
5994260
(8.26 secs, 23,596,249,112 bytes)
> sum . take 400 $ generateExponents' 2 3
5994260
(0.01 secs, 1,172,864 bytes)
> sum . take 1000000 {- a million -} $ generateExponents' 2 3
72001360441854395
(6.99 secs, 13,460,753,616 bytes)

I think you just forgot to map the actual function over the xs and ys:

generateExponents k l = [n | n <- [1 ..], n `elem` products n]
  where
    xs n = takeWhile (<= n) $ map (^ k) [1 ..]
    ys n = takeWhile (<= n) $ map (^ l) [1 ..]
    products n = liftA2 (*) (xs n) (ys n)