Qmxd1[x_, rho_, p_] :=
Module[{n = Floor[x], b, lambda, f, beta, i, j, k, z},
  b = Sum[i*p[[i]], {i,Length[p]}];
  lambda = rho/b;
  beta = p*lambda*Exp[-lambda*Range[Length[p]]];
  f[1,1] = (1 - rho);
  For[j = 2, j <= n+1, j++,
      f[1,j] = 0;
   ];
  For[i = 2, i <= n+1, i++,
      f[i,1] = Sum[f[(i-1),k], {k,n+1}];
      For[j = 2, j <= n+1, j++,
	  If [Length[p] >= (i-1),
	      f[i,j] = -1/(j-1) Sum[beta[[z]]*f[(i-z),(j-1)],{z,i-1}],
	      f[i,j] = -1/(j-1) Sum[beta[[z]]*f[(i-z),(j-1)],{z,Length[p]}]
	  ];
      ];
  ];
  1 - Exp[lambda x] (f[n+1,1] + Sum[f[n+1,k] (x-n)^(k-1), {k,2,n+1}])
]