0378 三角三元组
* * * *
拉格朗日计划
* * * *
三角三元组

记$T(n)=n(n+1)/2$为第n个三角数,并记$dT(n)$为$T(n)$的约数数量。

例如:$T(7)=28$,因此$dT(7)=6$。

再记$Tr(n)$为满足$1\le i < j < k\le n$以及$dT(i)>dT(j)>dT(k)$的三元组$(i, j, k)$的数量。

可以验证$Tr(20)=14, Tr(100)=5772$以及$Tr(1000)=11174776$。

求$Tr(60000000)$的最后18位数字。

本题难度:



解答

本题的思路很简单,但由于数据量很大,因此时间和空间开销都很高。

首先注意到n与n+1互素,因此$n(n+1)$的约数数量就是两者约数数量之积。

因此先用筛法找出不超过$\sqrt{60000000}$,再用筛法找出每个n的质因数分解,并计算其约数数量。

若$n=p_1^{a_1}\ldots p_k^{a_k}$是n的质因数分解,那么n有$(1+a_1)(1+a_2)\ldots(1+a_k)$个约数。由于计算$T(n)$时还需要除2,因此当$p_1=2$时取$a_1(1+a_2)\ldots(1+a_k)$。

统计三元组总数的方法是遍历中间的数,若其前面有a个约数数量大于它的数,后面有b个约数数量小于它的数,那么以该数为中间数的三元组共有ab个。

实现时先正向遍历一次,用有序列表记录每个数的约数数量,这样就可以通过二分查找来确定当前数之前多少个有约数数量大于其约数数量的数,用一个数组保存这一信息。

再逆向遍历一次,同样用有序列表记录每个数的约数数量,并通过二分查找来确定当前数之后多少个有约数数量小于其约数数量的数,并汇总更新。

最终结果是$147534623725724718$。

注:因需要使用sortedcontainers库,以下代码为Python 3,且代码中打印了进度信息。

注2:本题内存开销非常大,因此最后一步逆向遍历的同时还需要不断弹出约数列表和正向遍历所得的列表这两者末端不再使用的数据。

from sortedcontainers import SortedList

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

print("primes generated")

target=60000000

nums=list(range(target+2))
d=[1]*(target+2)
for j,p in enumerate(primeList):
    for n in range(p,target+2,p):
        k=0
        while nums[n]%p==0:
            nums[n]//=p
            k+=1
        d[n]*=(k if p==2 else k+1)
    if j%10==0:
        print("generating divisors", 1+j//10, "percent completed")

for n in range(2,target+2):
    if nums[n]>1:
        d[n]*=2
nums.clear()
print("divisors generated")

mod=10**18
infty=99999999
tick=600000

s=SortedList([0])
a=[0,0]
r=0
for n in range(2,target):
    c=d[n]*d[n+1]
    i=s.bisect_right(c)
    a.append(n-1-i)
    s.add(c)
    if n%tick==0:
        print("traversing", n//tick, "percent completed")
print("traversing completed")

s=SortedList([d[target]*d[target+1]])
d.pop()
r=0
for n in range(target-1,2,-1):
    c=d[n]*d[n+1]
    d.pop()
    i=s.bisect_left(c)
    r=(r+a.pop()*i)%mod
    s.add(c)
    if n%tick==0:
        print("computing", 100-n//tick, "percent completed")

print(r)