This is effectively partitioning a multiset of prime factors.

The simple algorithm to partition an integer is via a recursive function which takes as parameters the value remaining to partition, the maximum part, and suitable accumulator(s). E.g. if the handling of the final partition is done with a visit function then

_partitionInteger(int n, int maxPart, int[] accum)
    if (n == 0) visit(accum)
    else for (int m = 1; m <= maxPart; m++) _partitionInteger(n - m, m, accum :: m)

This is then wrapped for nice presentation as

partitionInteger(int n)
    _partitionInteger(n, n, [])

Exactly the same approach works here. The difference is that instead of dealing with integers we’re dealing with multisets. Instead of n == 0 we want to test t == {}, and instead of iterating from 1 to maxPart we need to iterate over multisets in lexicographic order. It will be easiest to represent them in frequency notation with an order of parts: for your application, that means that each multiset is an ordered sequence of prime powers. If the argument n has its first non-zero value at position i then the multisets m over which we iterate must have a non-zero value at position i: otherwise we’ll never terminate. (This is analogous to starting the iteration at m = 1 instead of m = 0).

So to take your example of 18 = 2^1 3^2, the recursive calls are:

_partition({2:1, 3:2}, {2:1, 3:2}, [])
    _partition({2:0, 3:2}, {2:1, 3:0}, [{2:1, 3:0}])
        _partition({2:0, 3:1}, {2:0, 3:1}, [{2:1, 3:0}, {2:0, 3:1}])
            _partition({2:0, 3:0}, {2:0, 3:1}, [{2:1, 3:0}, {2:0, 3:1}, {2:0, 3:1}])
                visit([{2:1, 3:0}, {2:0, 3:1}, {2:0, 3:1}])
        _partition({2:0, 3:0}, {2:0, 3:2}, [{2:1, 3:0}, {2:0, 3:2}])
            visit([{2:1, 3:0}, {2:0, 3:2}])
    _partition({2:0, 3:1}, {2:1, 3:1}, [{2:1, 3:1}])
        _partition({2:0, 3:0}, {2:0, 3:1}, [{2:1, 3:1}, {2:0, 3:1}])
            visit([{2:1, 3:1}, {2:0, 3:1}])
    _partition({2:0, 3:0}, {2:1, 3:2}, [{2:1, 3:2}])
        visit([{2:1, 3:2}])

Actually writing out the details gets a bit messier, but e.g. in Python we get:

def partitionFactors(n):
    primes, powers, p = [], [], 2
    while p * p <= n:
        pow = 0
        while n % p == 0:
            n /= p
            pow += 1
        if pow > 0:
            primes.append(p)
            powers.append(pow)
        p += 1 + (p & 1)
    if n > 1:
        primes.append(n)
        powers.append(1)
    return partitionFactorsImpl(primes, powers, list(powers), [])

def partitionFactorsImpl(primes, powers, max, accum):
    curr, firstIdx = [0] * len(primes), 0
    while firstIdx < len(primes) and powers[firstIdx] == 0:
        firstIdx += 1

    if firstIdx == len(primes):
        yield list(accum)
        return

    curr[firstIdx] = 1
    while True:
        i, factor = firstIdx, 1
        while i < len(powers):
            factor *= int(primes[i] ** curr[i])
            powers[i] -= curr[i]
            i += 1

        accum.append(factor)
        for partition in partitionFactorsImpl(primes, powers, curr, accum):
            yield partition
        accum.pop()

        i = firstIdx
        while i < len(powers):
            powers[i] += curr[i];
            i += 1

        # Advance curr. This is effectively counting in a variable base to avoid exceeding the available powers.
        idx = len(powers) - 1
        while idx >= firstIdx:
            curr[idx] += 1
            if curr[idx] > powers[idx]:
                curr[idx] = 0
                idx -= 1
            else:
                break
        # There are two ways of running out of space
        # Note that the test `curr > max` can probably be optimised
        if idx < firstIdx or curr > max:
            break

Sample output:

for x in partitionFactors(4 * 9 * 17):
    print (x)

gives

[2, 2, 3, 3, 17]
[2, 2, 51, 3]
[2, 2, 9, 17]
[2, 2, 153]
[34, 2, 3, 3]
[34, 2, 9]
[6, 2, 3, 17]
[6, 2, 51]
[6, 34, 3]
[6, 6, 17]
[102, 2, 3]
[102, 6]
[18, 2, 17]
[18, 34]
[306, 2]
[4, 3, 3, 17]
[4, 51, 3]
[4, 9, 17]
[4, 153]
[68, 3, 3]
[68, 9]
[12, 3, 17]
[12, 51]
[204, 3]
[36, 17]
[612]

(Note that the C# code in an earlier revision of this answer was buggy: the error is in the way max is handled).

In very general terms:

1. write out all prime factors of the number of interest in a long array or list (i. e. instead of [(2, 3) (5,2)], use [2, 2, 2, 5, 5])
2. map a function over this that returns all possible permutations of this list.
3. (to increase performance, filter out duplicates in the resulting list-of-lists)
4. Map over each of these permuted lists a function that returns all different ways the list can be split. (after the first elem, the second,… the n-1th)
5. For each split in each permutation, take the product of all factors in the left split, take the product of all factors in the right split.
6. Flatten the final list of list of number-pairs to end up with a list of number pairs that together will form your required number.

Edit: Above algorithm works if you want splits with two numbers. By introducing multiple splits in step 4, you can create the rest of the results.

I will ignore factor 1 because you could add it as many times as you want.

To solve the problem you need a procedure that produces the list of factorizations of N with the additions requirement that all factor are below or equal some max M.

To produce the factorizations of N

• Build the sorted list of divisors.
• Get the largest divisor D.
• produce the list of factorizations of N with factors <= D (see below)

To produce the factorization of N with factors <= M:

• For each divisor D if N with 1 < D <= M (use the sorted list create before)
  • If D=N, return the factorization (D).
  • else
    • Generate the factorization of N/D with factors <= D (recursive call)
    • Add factor D at the end of each returned factorization
• collect all factorizations from the loop into a single list and return that.

That could work. I haven’t actually tested it.