@Carl Love Here are some more-primitive methods for the Cartesian product of an arbitrary number of lists or vectors. The first uses modular arithmetic and represents each entry of the product as a multi-radix integer with the radices being the number of entries in the lists and the digits being the list entries:

CartProd_MultiRadix:= proc(V::seq({Vector, list}))
local K, i, v, B:= seq(numelems(v), v= V);
   [seq([seq(V[i][1+irem(K, B[i], 'K')], i= 1..nargs)], K= 0..mul([B])-1)]
end  proc:

A slight variation of the above should work in just about any programming language. Note that the 3-argument form of irem (integer remainder) used above returns the quotient in the third-argument, which is why it's in quotes ('K').

My next method corresponds to how I personally do Cartesian products by hand: Lay out all the first entries, then all the second entries, etc.:

CartProd_byRows:= proc(V::seq({Vector, list}))
local 
   i, j, k, v, 
   B:= [seq(numelems(v), v= V)], n:= nops(B), P:= seq(mul(B[1..k]), k= 0..n), N:= P[-1]
;
   convert(
      rtable([seq([seq(seq(V[k][irem(j,B[k])+1], i= 1..P[k]), j= 0..N/P[k]-1)], k= 1..n)])^%T,
      list, nested
   )
end proc:

The ^%T is the transpose operator (switching rows and columns of a matrix). The "-1" entry of an indexed structure (as in P[-1]) means the last entry.

My third method builds the appropriate nested seq command (to an arbritrary nesting depth) and then evaluates it:

CartProd_NestedSeq:= proc(V::seq({Vector, list}))
local i,k,S;
   eval(subs(S= seq, [foldl(S, [seq(i[k], k= 1..nargs)], seq(i[k]= V[k], k= 1..nargs))]))
end proc:

The command foldl (see ?foldl) takes a function (or function symbol, S in this case) and builds a nested invocation of that function. A variation of this method should be possible in most languages that support functional programming style. The eval and subs are only needed due to idiosyncracies of Maple (because seq behaves differently from most other functions).

If you're a serious student of computer programming, I recommend that you study these methods.

@vv Sorry. It is not my purpose to befuddle or obfuscate. On the other hand, I'm only interested in solutions that work for an arbitrary number of factors. In this case, that necessitates some complexity, and writing some code for this is an exercise that I recommend for any student of programming. (And I have three more completely different solutions that I'll post as a Reply to my Answer below.)

So, here's a clearer variation of the code above (the method based on Array or rtable indexing):

CartProd1A:= proc(V::seq({Vector, list}))
local v, i;
   [seq(
      rtable(
         seq(1..numelems(v), v= V), (k::seq(posint))-> [seq(V[i][k[i]], i= 1..nargs)]
       )
   )]
end proc:

@Matt C Anderson Thanks. Maple has a vast variety of predefined types, several tools for defining your own, and many commands that use types (such as the very important commands indets and subsindets). Types are at the heart of symbolic computation. A type is a boolean predicate that is applied to an expression. 

@acer No offense taken, even if you had meant that I had missed something, because I often do.

Here's some more-elaborate (non-uniform) examples, as requested. But this procedure cannot handle the first Example at ?DiscreteValueMap, because the support in that case is infinite.

I added an internal remember table to the procedure. It's just for efficiency; the functionality is the same.

restart
:
St:= Statistics
:
CombineDiscrete:= proc(R::list(RandomVariable), f)
uses S= Statistics, It= Iterator;
local 
   r, X, t, P:= table(sparse), nR:= nops(R), fX,
   Pr:= proc(e::(RandomVariable=realcons)) option remember; S:-Probability(e) end proc
;
   for X in 
      It:-CartesianProduct(
         seq([solve~(op~(indets(S:-PDF(r,t), specfunc(Dirac))), t)[]], r= R)
      )
   do
      P[(fX:= f(seq(X)))]:= P[fX] + mul(Pr(R[r]=X[r]), r= 1..nR)
   od;
   P:= [entries(P, 'pairs')];
   S:-RandomVariable(EmpiricalDistribution(lhs~(P), 'probabilities'= rhs~(P)))
end proc
:  
Dice:= 'St:-RandomVariable(DiscreteUniform(1,6))' $ 2: #quotes make them independent
Game:= CombineDiscrete([Dice], (x,y)-> piecewise(x=1, -4, y=1, -4, abs(x-y))): 
`2D6`:= CombineDiscrete([Dice], `+`): #sum of two six-sided dice

We add a "double-or-nothing" kicker to the game: If a second roll of a pair
of dice has an even sum, the payoff (either positive or negative) is 
doubled; if odd, the payoff is zero.
`2or0`:= CombineDiscrete([Game, `2D6`], (g, `2d6`)-> piecewise(`2d6`::even, 2*g, 0)):
St:-Mean(`2D6`);
                               7
St:-StandardDeviation(`2D6`, numeric);
                          2.415229458
St:-Mean(`2or0`);
                               -1
                               --
                               9 
St:-Variance(`2or0`);
                              1241
                              ----
                               81 

Now we consider the sum of 4 plays of the double-or-nothing game:
`4G`:= CombineDiscrete(['`2or0`'$4], `+`):
S:= St:-Sample(`4G`, 2^17):
St:-Histogram(S, gridlines= false);

@acer I said "with finite support", which is a subset of "from discrete distributions". I'm using support as per its usual mathematical definition rather than as it's used by Statistics. I don't see a way to handle a discrete distribution with infinite support.

@acer Could the OP be thinking about the CUDA package (?CUDA)?

I know that here you're just experimenting with MutableSets, but is any type of set, mutable or not, appropriate for this task? I mention this because in some of your other Questions I'm sure that you've used sets for tasks for which they were they were ill suited, where lists would've been better.

A set is a good choice over a list if the slightly increased time to build it is more than compensated by the substantial savings in element look-up times. A desire to keep numbers in order is not a good reason to use a set, IMO.

@vv I think that this is a substantially simpler inverse Fibonacci computation:

InvFib:= proc(x::nonnegative)
uses C= combinat;
local r, rr; 
   if x=0 then return 0 fi;
   r:= evalf(ln(sqrt(5)*x)/ln((1+sqrt(5))/2));
   if x::posint and x = C:-fibonacci((rr:= round(r))) then return rr fi;
   r
end proc
: 

It returns nonnegint if the result is exact, float otherwise. This simplification is possible because the second root of the characteristic equation (or, equivalently, the second eigenvalue of the matrix) has magnitude less than 1, and therefore it can be ignored for integer computation.

These Fibonacci methods are ad hoc to the specifc matrix from the original Question, and are difficult or impossible to generalize to some larger class (such as square matrices of nonnegint with an eigenvalue greater than 1). That is why I didn't address them in my original Answer even though I recognized the matrix as the generator for that sequence.

@izhammulya I made your follow-up question into a new Question with the title "Combining lists."

@tomleslie There are at least four reasons to not do it that way:

• Its average-case time is several times more than that of my procedure
• Its best-case time is several times more than that of my procedure
• It behaves ungracefully when the upperbound is greater than the sum of all the elements
• It's not a procedure

All of these things are obvious without needing to run tests, but feel free to do so anyway.

@tomleslie Unfortunately Dice1() + Dice1() does automatically simplify to 2*Dice1(), and this is a common pitfall for generating samples and performing other operations with multivariate r.v. expressions. You can use unevaluation quotes to get independent draws:

'Dice1()' + 'Dice1()'

@vv I would expect and hope that Maple would consider the complex case as the default. But here's a counterexample to that:

limit(exp(-1/x^2), x= 0);
                                                    0

Compare with

limit(exp(-1/x^2), x= 0, complex);
                                           undefined

@Kitonum Numbers of the same general "type" are ordered the usual way within a set, but they get grouped by type. Observe:

{0, 1/2, 1};
                                         {0, 1, 1/2}

Regardless, one should not rely on the order because it may change in different versions.

@Kitonum It seems to me to be mere happenstance that each of the last three elements (as they're listed above) is greater than each of the first three elements.

@tomleslie Your randomize should be randomize().