0157 分数不定方程
* * * *
拉格朗日计划
* * * *
分数不定方程

考虑不定方程 $$\frac{1}{a}+\frac{1}{b}=\frac{p}{10^n}$$ 其中a、b、p、n均为自然数且$a\le b$。若$n=1$,则该方程有20组解: \begin{align*} \frac{1}{1}+\frac{1}{1}&=\frac{20}{10} \\ \frac{1}{1}+\frac{1}{2}&=\frac{15}{10} \\ \frac{1}{1}+\frac{1}{5}&=\frac{12}{10} \\ \frac{1}{1}+\frac{1}{10}&=\frac{11}{10} \\ \frac{1}{2}+\frac{1}{2}&=\frac{10}{10} \\ \frac{1}{2}+\frac{1}{5}&=\frac{7}{10} \\ \frac{1}{2}+\frac{1}{10}&=\frac{6}{10} \\ \frac{1}{3}+\frac{1}{6}&=\frac{5}{10} \\ \frac{1}{3}+\frac{1}{15}&=\frac{4}{10} \\ \frac{1}{4}+\frac{1}{4}&=\frac{5}{10} \\ \frac{1}{4}+\frac{1}{20}&=\frac{3}{10} \\ \frac{1}{5}+\frac{1}{5}&=\frac{4}{10} \\ \frac{1}{5}+\frac{1}{10}&=\frac{3}{10} \\ \frac{1}{6}+\frac{1}{30}&=\frac{2}{10} \\ \frac{1}{10}+\frac{1}{10}&=\frac{2}{10} \\ \frac{1}{11}+\frac{1}{110}&=\frac{1}{10} \\ \frac{1}{12}+\frac{1}{60}&=\frac{1}{10} \\ \frac{1}{14}+\frac{1}{35}&=\frac{1}{10} \\ \frac{1}{15}+\frac{1}{30}&=\frac{1}{10} \\ \frac{1}{20}+\frac{1}{20}&=\frac{1}{10} \\ \end{align*} 对$1\le n \le 9$,该方程一共有多少组解?

本题难度:



解答

将原方程写做 $$\frac{1}{ap}+\frac{1}{bp}=\frac{1}{10^n},$$ 则根据第108题第110题的分析有 $$ap=10^n+x,\quad bp=10^n+y,$$ 其中 $$10^{2n}=xy.$$ 因此找出$10^{2n}$的约数对$x,y$,并计算 $$d=\gcd(10^n+x,10^n+y).$$ d的每个约数都可以作为p。

用筛法生成不超过45000的素数以计算d的约数的个数,汇总后即得结果$53490$。
import math

target=45000

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
p=[k for k in range(2,target) if d[k]==0]

r=0
for n in range(1,10):
    m=10**n
    m2=m*m
    for i in range(2*n+1):
        for j in range(2*n+1):
            x=(5**i)<< j
            y=m2//x
            if x<=y:
                d=math.gcd(m+x,m+y)
                k=1
                u=0
                while u < len(p) and d>1 and d>=p[u]:
                    if d%p[u]==0:
                        z=0
                        while d%p[u]==0:
                            d//=p[u]
                            z+=1
                        k*=(z+1)
                    u+=1
                if d>1:
                    k*=2
                r+=k
    print(n,r)

print(r)