Would someone be able and willing to help to address the problem (which is reverse to the approximation solution described in one of the answers to https://math.stackexchange.com/questions/5030561 question) by writing the program, which will create four files for the four groups, containing in total 1000 OEIS A223881 prime numbers terms, where terms in the each group are prime numbers found just by using approximation formula y = C ∗ x^m, where “y” is the function real number value (which is the closest value to some prime in the OEIS A223881) for the integer argument “x”, iterated from 1 to 1000, “m”- is the slope (all four curves have the same slope: m=1.197311990); “C” gives the displacement for the each of the curve. C1 =6.86845 C2 =3.420580 C3 =2.28335 C4 =1.70460 It appears to me that coding such a program is not a trivial task both mathematically and programmatically and is beyond my skills in those areas…
I am attaching 4 text files each containing prime numbers, belonging to
one of the 4 groups. Those four groups (files) were generated by Kristian Pontoppidan Larsen using unknown to me method. The task of the program which I am trying to write is to recreate generation of those 4 groups (files).
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:
Please correct me if my interpretation of your results is wrong but as the follow up (in attempt to connect your findings with the previously found approximation formulas), should C1, C2, C3 and C4 values in the approximation formulas be, at least, readjusted (taking in the account that per your findings the ratios should be exact natural numbers: 2, 3, 4, … ) using a simple average to determine the base value as:
C_1’ = 6.844515
C_2’ = 3.4222575
C_3’ = 2.281505
C_4’ = 1.71112875
and then to have the general (for all distinct curves) approximation formula as:
y = B * (p(m - 1) + 1)/L ∗ x^D
where
B ~ = 1.71112875,
D ~= 1.197311990
L is the largest prime factor of the
p(m - 1) + 1
?
Could we somehow figure out what is the real mathematical meaning concealed behind the empirically found B and D approximate values?