0428 圆环项链
* * * *
拉格朗日计划
* * * *
圆环项链

在一条直线上依次取W, X, Y, Z四点,使得$|WX|=a, |XY|=b, |YZ|=c, |WZ|=a+b+c$且$a,b,c\in\mathbb N$。

考虑分别以XY和WZ为直径的两个圆C和C',若存在另外k个圆($k\ge 3$) $C_1,\ldots, C_k$满足:

所有这k个圆都同时与C和C'相切,这k个圆中的任意两个都不相交,但$C_1$与$C_2$相切,$C_2$与$C_3$相切,。。。,$C_k$与$C_1$相切,

就称$(a,b,c)$为项链三元组。例如,(5, 5, 5)和(4, 3, 21)都是项链三元组,而(2, 2, 5)则不是。



记$T(n)$为满足$b\le n$的项链三元组(a, b, c)的数量,可以验证$T(1)=9, T(20)=732, T(3000)=438106$。

求$T(10^9)$。

本题难度:



解答

本问题的圆环项链即Steiner圆族(可以参考例如此贴),以下解法改编自官方论坛:

根据文中的inversion可以得到如下关系: $$\sin(\frac{\pi}{n})=\sqrt{\frac{ac}{(a+b)(b+c)}},$$ 从而 $$\frac{b(a+b+c)-ac}{b(a+b+c)+ac}=0,\pm\frac{1}{2},\pm1.$$ 接下来的算法似乎出自该论文:Diophantine Steiner Triples and Pythagorean-Type Triangles

根据其中的公式,只需找出每个n的素因子的数量,再根据n模2、3、6的余数分别带入sub1、sub2、sub3、sub4计算。

由于数据范围很大,因此按每组$10^7$个数的方式分组运用筛法计算素因子的数量,最终结果是$747215561862$。(代码中打印了进度信息)

N=1000000000
window=10**7

target=31622
d=[0]*target
k=2
while k < target:
    for i in range(k*k,target,k):
        d[i]=d[i]+1
    k+=1
    while k < target and d[k]>0:
        k+=1
primeList=[k for k in range(2,target) if d[k]==0]

res=0
end=1
while end < N:
    start=end
    end+=window
    f=[0]*window
    d=list(range(start,end))
    for p in primeList:
        a=(p-start%p) if start%p>0 else 0
        for i in range(start+a,end,p):
            f[i-start]+=1
            while d[i-start]%p==0:
                d[i-start]//=p
  
    for i in range(start,end):
        x=i-start
        if d[x]>1:
            f[x]+=1
        sub1=1+(i%3>0)
        sub2=1+(i%3==2)
        sub3=0
        sub4=0+(i%3 < 2)
        if i%2==1:
            sub1*=3
            sub2*=2
            sub4*=2
            if i%3==1:
                sub3+=1
            else:
                sub4+=1
        total=(sub1+sub2+sub3)*(N//i)+sub4*(N//(3*i))
        res+=total*(1 << (f[x]-1)) if f[x]>0 else total//(1 << (1-f[x]))
    if (end-1)%window==0:
        print end//window, "percent completed, current sum:",res

print 2*res