This is not a rigorous proof, but since
$$\frac{m!}{p_{m-1} + 1}$$ will only be a fraction when
$p_{m-1} +1$ has a prime factor that cannot be cancelled by
$m!$, and
$m!$ has enough prime factors to cancel any prime factors of
$p_{m-1}+1$ less than
$m$, the values of your OEIS sequence can be found by identifying where the largest prime factor of
$p_{m-1}+1$ is greater than
$m$. This will always be prime since this largest factor can only divide
$p_{m-1}+1$ once since
$p_{m-1}+1 < m^2$
Therefore the sequence can be found by using the LargestPrimeFactor resource function, and identifiying where Prime[m - 1] + 1 has a largest prime factor greater than m (there is one exception when m=3):
ClearAll["`*"]
lpf = ResourceFunction["LargestPrimeFactor"]
sequenceIndex = 1;
dat = Table[
denom = Prime[m - 1] + 1;
largest = lpf[denom];
If[largest > m || m == 3,
{denom/largest, sequenceIndex++, largest}
,
Nothing]
, {m, 2, 6000}
];
Note that I keep track also of denom/largest, the value of
$p_{m-1}+1$ divided by its largest factor. If we gather the sequence terms by denom/largest, we see we recover the distinct families:
gathered = #[[All, 2 ;; All]] & /@ GatherBy[dat, First];
ListLinePlot[gathered]
These families correspond to when
$p_{m-1}+1 = 2L , 4L, 6L.... $ where
$L$ is the largest prime factor of
$p_{m-1}+1$.
