问题提出

2008年11月medie2005在数学研发论坛提出:
\pi(n)表示不大于n的素数个数.
比如,\pi(13)=6. (不大于13的素数有:2,3,5,7,11,13)

13有一个很有意思的性质:\pi(13)=6=1! \times 3!,即:\pi(13)等于13的数字组成(1,3)的阶乘的乘积.
若自然数n满足:\pi(n)等于n的数字组成的阶乘的乘积,我们就称n为PF数.
有同样性质的数还有1512,1520,1521等等.

更难得的,\pi(226130351)=2!\times 2!\times 6!\times 1!\times 3!\times 0!\times 3!\times 5!\times 1!,并且,226130351是素数.
如果n既是PF数,又是素数,则称n为PFP数.

问题:
1):求出10^{18}内的所有PF数.
2): 13,226130351是前2个PFP数,请求出接下来的两个PFP数.

进展

medie2005说问题来源自A066457,
他指出
注意到,数字组成的阶乘的乘积是与顺序无关的。
因此,10^{18}内的的数的数字组成的阶乘的乘积的形式一共只有:C_{27}^9=4686825个。
现在关键是计算\pi(n)
他还给出了一个只包含2829个候选素数的列表,只是可惜这个附件已经无法下载了。

计算\pi(n)

medie2005找到一份源自张一飞O(n^{\frac23+\epsilon})复杂度的计算\pi(n)的算法

#include      
#include      
#include      

const   unsigned   default_n=4000000000;   
const   int   maxlen=1000;   

int   primes[maxlen];   
int   len;   
int   sum;   
unsigned   int   n;   
int   m=7;   
int   Q=1;   
int   phiQ=1;   
int\times    v;   

bool   string2int(unsigned   int&   n,char\times    s);   
void   init(void);   
unsigned   int   phi(unsigned   int   x,int   a);   

void   main(void){   
        char   number_string[80];   
        unsigned   int   time;   

        cout<<"calc   pi(n)\n";   
        cout<<"n(default   =   "<n){   
                                if(!mark[j])++s;   
                                --j;   
                        }   
                        sum+=s;   
                }   
        }   

        for(len3=len2;ilen)m=len;   
        for(i=1;i<=m;++i){   
                Q\times =primes[i];   
                phiQ\times =primes[i]-1;   
        }   

        v=(int\times )GlobalAlloc(GMEM_FIXED,Q\times sizeof(int));   
        for(i=0;i=0;--j){   
                        v[j]-=v[j/primes[i]];   
                }   
        }   

        GlobalFree(mark);   
}   

bool   string2int(unsigned   int&   n,char\times    s){   
        unsigned   int   val=0;   
        for(int   i=0;s[i];++i){   
                if(s[i]>'9'||s[i]<'0')return   false;   
                val\times =10;   
                val+=s[i]-'0';   
                if(val>default_n)return   false;   
        }   
        n=val;   
        return   true;   
}   

unsigned   int   phi(unsigned   int   x,int   a){   
        if(a==m){   
                return   (x/Q)\times phiQ+v[x%Q];   
        }   
        if(x<(unsigned)primes[a]){   
                return   1;   
        }   
        return   phi(x,a-1)-phi(x/primes[a],a-1);   
  }  

成果展示

medie2005于2008年11月20日一下子找到了一组四胞胎数据

PF: 8360755200   209210612202
PF: 8360755200   209210612212
PF: 8360755200   209210612220
PF: 8360755200   209210612221

格式说明:
PF X Y
X表示\pi(Y)
并且他把结果提交到了A066457
2017年3月22日alt找到了
\pi(13030323000581525)=361184624640000=1!\times 3!\times 0!\times 3!\times 0!\times 3!\times 2!\times 3!\times 0!\times 0!\times 0!\times 5!\times 8!\times 1!\times 5!\times 2!\times 5!
而这个结果后来也被收录到A066457