Lwins_G 发表于 2014-3-21 12:01
直接计算\[ \sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \le ...
$\psi(x,y)$是$O(\sqrt(x))$的，那$\mu(d)$呢，我怎么感觉它的计算量也不低呀
本帖最后由 sunwukong 于 2014-6-26 10:54 编辑
由 13#，
\(n^2=st,s \lt n \lt t\le L-n\)
设 \(s=k*a^2\),其中，\(k\) 无平方素因子，由整数分解定理可知，这个分解是唯一的。
由\(n^2=k*a^2*t\)
得，\(t=k*b^2,n=k*a*b\)，由 \(s\lt t\)，得 \(a\lt b\)，
由\(t\le L-n\)，得\(k*b*(a+b)\le L\)
所以，原问题转化为
给定 \(L\) ，求满足
\(\mu(k)\not=0\)
\(k*b*(a+b)\le L\)
\(a\lt b\)
的解 \(\left(k,a,b \right)\) 的个数。
因为 \(b*(a+b)>=2*(1+2)=6\)，所以 \(k\le L/6\)
所以
\[
\begin{split}
f(L)=\sum_{\begin{subarray}{cc} k=1 \\ \mu(k)\not=0 \end{subarray}}^{L/6} \sum_{b=2}^{\sqrt{L/k}} \min\left( b-1,\lfloor \frac{L}{kb} \rfloor-b \right)
\end{split}
\]
注： \(\mu(x)\) 是莫比乌斯函数 http://zh.wikipedia.org/wiki/默比乌斯函数
本帖最后由 282842712474 于 2014-6-26 17:14 编辑
Lwins_G 发表于 2014-3-21 12:01
直接计算\[ \sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \le ...
刚刚写完了$O(L^{3/4})$的公式的代码，运行测试$f(10^9)$的时间为24s左右（由上次的190s到了24s）。
但是$f(10^{10})$运行了半个钟都没有结果，而运行$f(2\times 10^9)$需要75s
我觉得由于$d|n$的存在，需要分解n，两个公式其实都不止$O(N)$
$O(L^{3/4})$算法的代码如下（还是python)
import time
import os
start=time.clock()
import math
#主要计算函数（最后一个求和号）
def psi(l, n, d):
sum=0
for s in xrange(n//d+1, (2*n-1)//d + 1):
sum = sum + l//(n*d*s)
return sum
#生成一个素数表的莫比乌斯表
def mobi(p):
length = len(p)
if length == 1:
mu = ]
else:
mu2=mobi(p)
mu=mu2[:]
for i in mu2:
mu.append(-1*p*i)
return mu
#开始
l=2000000000 #上限
s=0 #统计个数
rl=int(math.sqrt(l))
rrl=int(math.sqrt(int(math.sqrt(l))))
rrrl=int(math.sqrt(int(math.sqrt(int(math.sqrt(l))))))
#生成rrl内的素数表
prime=#定义整数表
j=2
while j<=rrrl:
if prime != 0:
m=int(rrl/j)
for i in xrange(j,m+1):
prime=0
j=j+1
else:
j=j+1
prime.sort() #从小到大重新排列（让0位于前面）
z=prime.count(0) #统计0的个数
prime=prime
#素数表生成完毕
prime.append(rl) #往素数表里边添加一个大数，避免溢出错误
for n in xrange(2,rl+1):
p=[] #n的素因子表
mu=[] #莫比乌斯表
nn=n
i=0
rn=int(math.sqrt(n))
#开始生成n的素因子表
while prime <= rn:
if n%prime==0:
n=n//prime
if n%prime != 0:
p.append(prime)
i=i+1
else:i=i+1
#下面这句的意思是，如果n的一个因子不能被前面的素数整除，那么这个因子就是素数。
if n>rn:
p.append(n)
#n的素因子表p生成完毕
mu=mobi(p) #生成莫比乌斯表
for d in mu:
s=s+d//abs(d)*psi(l,nn,abs(d))
print(s)
end=time.clock()
print(end-start)
我发觉主要的问题是，我想不出$psi(x,y)$的$O(\sqrt{x})$算法...求指教。
本帖最后由 282842712474 于 2014-6-28 00:33 编辑
终于做出来了。$f(10^{12})=5435004633092$，用时101秒！
稍后会整理代码出来。
=========================
http://spaces.ac.cn/index.php/archives/2665/
import time
start=time.clock()
import math
#主要计算函数（最后一个求和号）
def psi(x, y1, y2):
s=0
rx=int(math.sqrt(x))
if y2<= rx:
for d in xrange(y1, y2+1):
s=s+x//d
elif y1<=rx:
for d in xrange(y1, rx+1):
s=s+x//d
p=x//rx-1
q=x//y2
d=0
q1=x//(p-d+1)
while d < p-q:
q2=q1
q1=x//(p-d)
s=s+(q1-q2)*(p-d)
d=d+1
s=s+(y2-x//(q+1))*q
else:
p=x//y1-1
q=x//y2
if p<q:
s=(1+y2-y1)*q
else:
d=0
q1=x//(p-d+1)
while d < p-q:
q2=q1
q1=x//(p-d)
s=s+(q1-q2)*(p-d)
d=d+1
s=s+(y2-x//(q+1))*q
s=s+(x//(p+1)-y1+1)*(p+1)
return s
#生成一个素数表的莫比乌斯表
def mobi(p):
length = len(p)
mu=]
for i in xrange(1,length):
mu2=mu[:]
for j in mu2:
mu.append(-1*j*p)
return mu
#开始
l=1000000000000 #上限
s=0 #统计个数
rl=int(math.sqrt(l))
rrl=int(math.sqrt(int(math.sqrt(l))))
#生成rl内的素数表
prime=#定义整数表
j=2
while j<=rrl:
if prime != 0:
m=int(rl/j)
for i in range(j,m+1):
prime=0
j=j+1
else:
j=j+1
prime.sort() #从小到大重新排列（让0位于前面）
z=prime.count(0) #统计0的个数
prime=prime
#素数表生成完毕
prime.append(rl) #往素数表里边添加一个大数，避免溢出错误
for n in xrange(2,rl+1):
p=[] #n的素因子表
mu=[] #莫比乌斯表
nn=n
i=0
rn=int(math.sqrt(n))
#开始生成n的素因子表
while prime <= rn:
if n%prime==0:
n=n//prime
if n%prime != 0:
p.append(prime)
i=i+1
else:i=i+1
#下面这句的意思是，如果n的一个因子不能被前面的素数整除，那么这个因子就是素数。
if n>rn:
p.append(n)
#n的素因子表p生成完毕
mu=mobi(p) #生成莫比乌斯表
for d in mu:
ll=l//(nn*abs(d))
s=s+d//abs(d)*(psi(ll,nn//abs(d)+1,(2*nn-1)//abs(d)))
print(s)
end=time.clock()
print(end-start)