0351 六边形果园
* * * *
拉格朗日计划
* * * *
六边形果园

n阶六边形果园是由边长为n的正六边形分割而成的正三角形网格。例如下图是5阶六边形果园:



站在中间向外看,图中用绿色标出的点都会被更近的点挡住。可以看出,对于5阶六边形果园,从中间向外看时有30个点会被挡住。

记在n阶六边形果园的中间向外看时会被挡住的点的数目为$H(n)$。

可以验证$H(5)=30, H(10)=138, H(1000)=1177848$。

求$H(10^8)$。

本题难度:



解答

显然六边形果园就是由$x=(1,0)$和$y=(\cos\frac{2\pi}{3}, \sin\frac{2\pi}{3})$生成的格,即果园中每个点都可以唯一地写成$ax+by$的形式,其中$a,b\in\mathbb Z$。

若$d=\gcd(a,b)>1$,那么$(a,b)$就会被$(a/d,b/d)$挡住。

对每个固定的a,与a互素的数有$\varphi(a)$个,其中$\varphi$是欧拉Totient函数。

a取遍1到n之间的值,汇总$\varphi(a)$,再乘以6(对称性)即可得: $$H(10^8)=6\cdot\frac{10^8\cdot (10^8+1)}{2}-6\sum_{a=1}^{10^8}\varphi(a)$$ 其中欧拉Totient函数的和$\sum_{a=1}^{10^8}\varphi(a)$在OEIS A002088给出,此处使用的是该页面上PROG栏目中的代码。

最终结果是$11762187201804552$。

注:因用到记忆化搜素的cache装饰器,以下代码为Python 3。

from functools import *

@cache
def f(n):
    if n==0:
        return 0
    c,j=0,2
    m=n//j
    while m>1:
        k=n//m+1
        c+=(k-j)*(2*f(m)-1)
        j,m=k,n//k
    return (n*(n-1)-c+j)//2

print(3*(10**8)*(10**8+1)-6*f(10**8))