I have a sampling algorithm for you written in sage / python.

Instead of $\underline{123}$ avoiding permutations of length $n$ I consider words $w_1 w_2 \dots w_n$ satisfying the following properties:

- $0\le w_k \le n-k$
- the word is weakly $\underline{123}$ avoiding, which means there is no index $i\in [1,n-2]$ such that $w_i\le w_{i+1} \le w_{i+2}$.

I call these words $\underline{123}$ avoiding indexwords.
They are in bijection with $\underline{123}$ avoiding permutations using the following bijection: Let $S=\{1,2,\dots,n\}$ then for a given indexword $w_1 w_2 \dots w_n$ we define $\pi_i = w_i\text{-th smallest entry of } S\setminus\{\pi_1,\dots,\pi_{i-1}\}$, where I start indexing the entries of the set with $0$.
Put differently: If we keep track of the unused letters to build up the permutation in an ordered list, then $\pi_i$ is precisely the $w_i\text{-th}$ entry in that list (explaining my name for these words).

These words will simplify the following algorithm for random uniform sampling.
The algorithm comes in two steps: First count the number of $\underline{123}$ avoiding indexwords with a given prefix. Then create a random word iteratively by appending letters to the already constructed word using certain probabilities using the precomputed values.

For this let $\alpha_n^{w_1 w_2 \dots w_l}$ be the number of $\underline{123}$ avoiding indexwords of length $n$ starting with $w_1 w_2 \dots w_l$.
For calculating these values, we have two rules:
Let $w_1 w_2 \dots w_l$ be weakly $\underline{123}$ avoiding, then

- Expansion rule
$$
\alpha_n^{w_1 w_2 \dots w_l} = \sum_{0\le x\le n-l-1 \text{ and } \neg(w_{l-1}\le w_l\le x)} \alpha_n^{w_1 w_2 \dots w_l x}
$$
- Reduction rule
$$
\alpha_n^{w_1 w_2 \dots w_l} = \alpha_{n-l+2}^{w_{l-1}w_l}
$$
- Start cases
$$
\alpha_2^{00}=\alpha_2^{10}=1
$$

Combining the rules we obtain
$$
\alpha_n^{xy}=\sum_{0\le z\le n-3 \text{ and } \neg(x\le y\le z)}\alpha_{n-1}^{yz}
$$

Note that you can obtain $a_n$ by summing $\alpha_n^{xy}$ over all pairs $(x,y)\in [0,n-1]\times [0,n-2]$

Now for the sampling part.
Randomly choose a pair $(x,y)\in [0,n-1]\times [0,n-2]$ with probability $\frac{\alpha_n^{xy}}{a_n}$.
Assume now, we have already created a sequence $w_1 \dots w_l$ with $l<n$, we append an appendable letter $x\in [0,n-l-1]$ with probabilty
$\frac{\alpha_{n+1-l}^{w_l x}}{\alpha_{n+2-l}^{w_{l-1} w_{l}}}$.
(A letter is appendable, if it does not create the $\underline{123}$ pattern).
Repeat this until you have a $\underline{123}$ avoiding indexword of length $n$.

Here is the sage code:

```
@cached_function
def alpha(n,x,y):
assert n>1, "n needs to be at least 2"
if n == 2:
return 1
return sum(alpha(n-1,y,z) for z in range(n-2) if not (x<=y<=z))
def random_indexword(n):
assert n>=2, "Algorithm only works for length n>=2"
pairs = [(x,y) for x in range(n) for y in range(n-1)]
P = [alpha(n,x,y) for x,y in pairs]
s1,s2 = pairs[GeneralDiscreteDistribution(P).get_random_element()]
indexword = [s1,s2]
while len(indexword)<n:
appendable = [x for x in range(n-len(indexword)) if not (indexword[-2] <= indexword[-1] <= x)]
P = [alpha(n+1-len(indexword), indexword[-1],x) for x in appendable]
indexword.append(appendable[GeneralDiscreteDistribution(P).get_random_element()])
return indexword
def indexword_to_permutation(word):
letters = [i for i in range(1,len(word)+1)]
p = []
for i in word:
p.append(letters[i])
letters.remove(p[-1])
return p
# for testing
def is_123_avoiding(perm):
return not any(perm[i]<perm[i+1]<perm[i+2] for i in range(len(perm)-2))
w = random_indexword(10)
perm = indexword_to_permutation(w)
perm,is_123_avoiding(perm)
```

Edit: Updated typo in bijection between index words and permutations.