We present computational results on principal specializations $\mathfrak{S}_w(1^n)$ of Schubert polynomials, which count reduced pipe dreams and reduced bumpless pipe dreams (RBPD). We find the first counterexample, at $n=17$, to the Merzon-Smirnov conjecture (arXiv:1410.6857) that the maximum of $\mathfrak{S}_w(1^n)$ over $S_n$ is attained at a layered permutation. The simulations suggest that $\lim_{n \to \infty} \log(\max_{w\in S_n}\mathfrak{S}_w(1^n))/n^2$ equals the maximal layered permutations' constant from Morales-Pak-Panova (arXiv:1805.04341). We also explore the random permutation drawn from the distribution proportional to $\mathfrak{S}_w(1^n)$, revealing permuton-like asymptotics similar to those for Grothendieck polynomials by Morales-Panova-Petrov-Yeliussizov (arXiv:2407.21653). We implement and compare three recurrences for $\mathfrak{S}_w(1^n)$: the descent formula (Macdonald), transition formula (Lascoux--Schutzenberger), and cotransition formula (Knutson). For sampling uniformly random RBPDs (whose count is $\sum_{w\in S_n} \mathfrak{S}_w(1^n)$), we show that reducedness breaks the sublattice property of the ASM lattice, preventing monotone CFTP and causing false coalescence. We develop an efficient MCMC sampler with macroscopic "droop" updates for connectivity and fast mixing. Our code computes $\mathfrak{S}_w(1^n)$ up to $n\sim 20$ and samples random RBPDs up to $n\sim 60$ on a personal computer ($n\sim 100$ on a cluster).