simply code it using a table just as I coded power. Hence you have to know upper and lower bounds on m (btw, you talk about gamma distributions etc. The code I give you here implicity assume that you have models that besides these two complicating operators are simlpe linear/quadraticx programming implementable. YALMIP can not be used for arbitrarily complicated nonlinear operators etc, it is a language for very structured problems. before you start trying to implement these things, you should post your complete problem so I can judge if it has any chance of being implemented in YALMIP and wether there is any solver available for the problem
Constraints = ;
mfactorial = sdpvar(1);
for i 1:mupper
Constraints = [Constraints,implies(m==i,mfactorial==factorial(i))];
As I said regarding the previous operator, if it enters in a convex way, it can be represented without implies, i.e. additional binary variables. Hence, it is crucial that you post the complete model.