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 = [1,-p[0]] else: mu2=mobi(p[0:length-1]) mu=mu2[:] for i in mu2: mu.append(-1*p[length-1]*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=[i for i in xrange(1,rrl+1)]#定义整数表 j=2 while j<=rrrl: if prime[j-1] != 0: m=int(rrl/j) for i in xrange(j,m+1): prime[i*j-1]=0 j=j+1 else: j=j+1 prime.sort() #从小到大重新排列(让0位于前面) z=prime.count(0) #统计0的个数 prime=prime[z+1:rl] #素数表生成完毕 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[i] <= rn: if n%prime[i]==0: n=n//prime[i] if n%prime[i] != 0: p.append(prime[i]) 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)