Yes - it is called a binomial distribution. The probability of a trial having a probability of success of $p$ being successful $k$ times out of $n$ trials is
$$P(k) = \binom{n}{k} p^k (1-p)^{n-k}$$
For your numbers, the probability is
$$P(K \ge N) = \sum_{k=N}^M \binom{M}{k} p^k (1-p)^{M-k}$$
For your die example:
$$P(K \ge 75) = \sum_{k=75}^{100} \binom{100}{k} \left (\frac13\right)^k \left ( \frac{2}{3}\right )^{100-k} \approx 1.89 \cdot 10^{-17}$$
Here is the full distribution over all rolls, which explains why this probability is so small:
ADDENDUM
You may also have noticed that the distribution looks quite normal. In fact, for large $N$, it is well-known that binomial distributions (among others) approximate a normal distribution very well. In this case, $\mu = N p = 100/3$ and $\sigma^2 = N p (1-p) = 200/9$.
UPDATE: Seeing again the same ugly counting problem we can try a different approach using a ordinary generating function.
Let $s:=\{s_1,s_2,s_3,s_4\}$ the setup that we want to count where $s_k$ is the number of $k$-valued sides in the throw (by example if $s_3=5$ it means that the throw must contain exactly five $3's$.)
Because we are throwing $n$ dices then for each valid setup $s$ there are $\ell:=n-\sum_k s_k$ dice who values are $5$ or $6$. Then there is a maximum number $\ell$ of "wildcards" (that is, of possible $6's$). Let define the order relation $s\ge t$ if and only if $s_k\ge t_k$ for $k=1,2,3,4$. Then if $s\ge t$ then it is enough to have a valid throw when $\ell\ge s_6\ge\sum_k(s_k-t_k)=:|s-t|_1$.
Then for some $s\ge t$ there are $\ell+1-|s-t|_1$ possible valid setups. Now the ordinary generating function defined by
$$g_s(x):=\prod_{k=1}^4\left(\sum_{j=0}^{s_k}x^j\right)$$
is counting all $s\ge t$ for each $|t|_1$, that is $[x^{|t|_1}]g_s(x)$ is the number of $t\le s$ with $|t|_1:=\sum_k t_k$. And because for every $|t|_1$ there are $\ell+1-|s-t|_1=n+1-2|s|_1+|t|_1$ valid setups to count them all it is enough to make the dot product $(v_s,r_s)$, where $r_s$ is the vector defined by the coefficients of the polynomial $g_s$, that is
$$r_s:=([x^0]g_s(x),[x^1]g_s(x),\ldots,[x^{|s|_1}]g_s(x))$$
and $v_s:=(n+1-|s|_1,\ldots,2,1,0,\ldots,0)$.
Then the probability $p_s$ that some setup equivalent to $s$ appear in this game will be $p_s=(v_s,r_s)/6^n$. From this notation it could be easier to write efficient programs that evaluate these probabilities.
However if $n$ is too big probably an approximation could be better, by example from some Montecarlo simulation. In any case I doubt that a "closed" formula can be found.
Best Answer
The chance that $1$ comes up exactly $33$ times in $100$ comes from the binomial distribution. The chance of success is $\frac 16$ and failure is $\frac 56$ so it is ${100 \choose 33}(\frac 16)^{33}(\frac 56)^{67}\approx 0.00003$ If we sum from $33$ to $100$ we get the chance of at least $33\ 1$s, which is about $0.0005$ per Alpha. You can multiply these by $6$ to get the chance for any number, as it is very unlikely we doublecount by having at least $33$ of two different numbers. So the chance of a "lucky number" happening by chance is about $0.0003$ or one in $3300$. Pretty unlikely, but rarer things happen all the time.