尤拉計劃708:你只要2

設f(n)是把一個正整數的質因數全替換成2,所得的積。規定f(1)=1。
例如, 90=2×3×3×5, 然後替換質因數, 2×2×2×2=16, 因此 f(90)=16
設S(N)=∑ f(n),其中 n=1~N。已知S(10⁸)=9613563919。
求S(10¹⁴)

測試程式碼

import math
def prime(n):
 primes=[]
 h  = [True] * int(n+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** (1.0*0.5)) + 1):
    if h[i]:
        h[i*i::i] = [False] * len(h[i*i::i])
 for i, e in enumerate(h):
    if e:
        primes.append(i)
 return primes

def S(n):
 p=prime(n//2)
 m=[1]* int(n+1)
 L=len(p)
 for x1 in range(0,L):
  if p[x1]*p[x1]>n:break
  for a1 in range(1,1+int(math.log(n,p[x1]))):
   t1=p[x1]**a1
   if t1>n:break
   m[t1]=a1
   for x2 in range(x1+1,L):
    if p[x2]>n/t1:break
    for a2 in range(1,1+int(math.log(n,p[x2]))):
     t2=t1*p[x2]**a2
     if t2>n:break
     m[t2]=a1+a2
     for x3 in range(x2+1,L):
      if p[x3]>n/t2:break
      for a3 in range(1,1+int(math.log(n,p[x3]))):
       t3=t2*p[x3]**a3
       if t3>n:break
       m[t3]=a1+a2+a3
       for x4 in range(x3+1,L):
        if p[x4]>n/t3:break
        for a4 in range(1,1+int(math.log(n,p[x4]))):
         t4=t3*p[x4]**a4
         if t4>n:break
         m[t4]=a1+a2+a3+a4
         for x5 in range(x4+1,L):
          if p[x5]>n/t4:break
          for a5 in range(1,1+int(math.log(n,p[x5]))):
           t5=t4*p[x5]**a5
           if t5>n:break
           m[t5]=a1+a2+a3+a4+a5
           for x6 in range(x5+1,L):
            if p[x6]>n/t5:break
            for a6 in range(1,1+int(math.log(n,p[x6]))):
             t6=t5*p[x6]**a6
             if t6>n:break
             m[t6]=a1+a2+a3+a4+a5+a6
             for x7 in range(x6+1,L):
              if p[x7]>n/t6:break
              for a7 in range(1,1+int(math.log(n,p[x7]))):
               t7=t6*p[x7]**a7
               if t7>n:break
               m[t7]=a1+a2+a3+a4+a5+a6+a7
               for x8 in range(x7+1,L):
                if p[x8]>n/t7:break
                for a8 in range(1,1+int(math.log(n,p[x8]))):
                 t8=t7*p[x8]**a8
                 if t8>n:break
                 m[t8]=a1+a2+a3+a4+a5+a6+a7+a8
 #print(m)
 s=0
 for i in m:
  s+=2**i
 print(s-2-1)

import time
t=time.time()
S(10**8)
print(time.time()-t)

>>>> import time
>>>> t=time.time()
>>>> S(10**8)
9613563919
>>>> print(time.time()-t)
42.9400000572