0311 双斜整四边形
* * * *
拉格朗日计划
* * * *
双斜整四边形

ABCD是一个各边长均为整数的凸四边形,且满足$1\le AB < BC < CD < AD$。

设O是BD的中点,若BD和AO的长都是整数且$AO=CO\le BO=DO$,就称ABCD为双斜整四边形。

例如,以下四边形就是双斜整四边形:$AB=19, BC=29, CD=37, AD=43, BD=48, AO=CO=23$。



记$B(N)$是所有满足$AB^2+BC^2+CD^2+AD^2\le N$的双斜整四边形ABCD的数量。

可以验证$B(10^4)=49$以及$B(10^6)=38239$,求$B(10^{10})$。

本题难度:



解答

Apollonius定理可得 $$AB^2+AD^2=2(AO^2+BO^2)=2(CO^2+BO^2)=BC^2+BD^2$$ 又注意到 $$2(AO^2+BO^2)=(AO+BO)^2+(AO-BO)^2$$ 因此本题事实上是求满足 $$a^2+d^2=b^2+c^2=m^2+n^2$$ 的六元组的数量。对此也没有太好的手段,直接枚举平方和并统计频次,若 $$S_n=\{(x,y): x^2+y^2=n\},$$ 那么只要$|S_n|\ge 3$,其中任意三组元素都符合要求。

由于本题数据范围较大,实现时采用分段枚举的方式减少内存开销,最终结果是$2466018557$。

注:以下代码改写自官网论坛。

import math

N=5*10**9
r=0
S=10**8
for i in range(N/S+1):
    L=i*S
    U=(i+1)*S
    counters=[0]*(S/2)
    for x in range(int(math.sqrt(N))+1):
        y=max((int(math.sqrt(L-x*x))&~1)+(x&1)-4,x+2) if x*x < L else x+2
        z=x*x+y*y
        while y <=int(math.sqrt(N-x*x)) and z < U:
            if z>=L:
                counters[(z-L)/2]+=1
            y+=2
            z=x*x+y*y
    r+=sum(c*(c-1)*(c-2)/6 for c in counters)
    print i*2,"percent completed"

print(r)