131. 素数立方数组合(Prime cube partnership)

对于某些素数\(p\)存在正整数\(n\),使得表达式\(n^3+n^2p\)是一个立方数。例如对于\(p=19\)\(8^3+8^2\times19=12^3\)。非常奇特的是,对于满足这个性质的素数,\(n\)的值是唯一的,在小于一百的素数中,只有四个素数满足这个性质。求在小于一百万的素数中,有多少个素数满足这个神奇的性质?

分析:根据题意我们假设\(n^3+n^2p=k^3\)是一个立方数,则两边同时开立方可得:

$$ \sqrt[3]{n^{3}\left( 1+\frac{p}{n}\right)} =n\cdot \sqrt[3]{\frac{n+p}{n}} =n\cdot \frac{\sqrt[3]{n+p}}{\sqrt[3]{n}} =k $$

则显然\(n\)\(n+p\)均为立方数,设\(n=x^3,n+p=y^3\),则有\(p=y^3-x^3\),则有:

$$ p=y^3-x^3=(y-x)(y^2+yx+x^2) $$

\(p\)是一个素数,则只有\(1\)\(p\)两个因子,而显然\((y-x)<(y^2+yx+x^2)\),则有\(y-x=1,y^2+yx+x^2=p\),则:

$$ p=(x+1)^2+(x+1)x+x^2=x^2+2x+1+x^2+x+x^2=3x^2+3x+1 $$

\(p\)应当是一个素数,据题意有\(p=3x^2+3x+1<10^6\Rightarrow x<576.85\),我们只需要逐个尝试\(x\in[1,576]\),检测\(3x^2+3x+1\)是否是一个素数,并对所有符合条件的数统计个数就可以了。代码如下:

# time cost = 8.71 ms ± 878 µs

from sympy import isprime

def main(N=576):
    c = 0
    for i in range(1,N+1):
        if isprime(3*i**2+3*i+1):
            c += 1
    return c