尤拉計劃694:立方滿因數

A positive integer n is considered cube-full, if for every prime p that divides n, so does p3. Note that 1 is considered cube-full.

Let s(n) be the function that counts the number of cube-full divisors of n. For example, 1, 8 and 16 are the three cube-full divisors of 16. Therefore, s(16)=3 Let S(n) represent the summatory function of s(n), that is S(n)=∑i=1~n s(i)

You are given S(16)=19 , S(100)=126 and S(10000)=13344

Find S(10¹⁸).

如果一個正整數n,其所有質因數的立方也是它的因數，則稱其為立方滿數。注意1被當作立方滿數。

設s(n)是n的立方滿因數個數。例如，1、8、16是16的3個立方滿因數。所以s(16)=3。設S(n)表示s(n)的求和，即S(n)=∑_i=1ⁿ s(i)

已知 S(16)=19 , S(100)=126 ， S(10000)=13344

求S(10¹⁸).

考慮

所有質數都不是Cube-full 數
所有質數立方及更高次方都是Cube-full 數比如 8 16 64
所有Cube-full 數乘以另一個Cube-full 數的積也是 8x27=216
每個數都能表示成p1^n1*p2^n2... pm^nm,所以它的因數個數就是n1*n2...nm+1
只有大於1的Cube-full 數才有超過1個以上的Cube-full 因數，其他數n的s(n)都為1

以下程式碼算小於27000（30^3）的n的S(n)

import math
def Sai(n):
 primes=[]
 h  = [True] * int(n**(1.0/3)+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** (1.0/3*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)
 sums=n
 for p in primes:
  for i in range(0,int(math.log(n)/math.log(p))+1):
   sums+=n//(p**(i+3)) #p**3的倍數
 for p1 in primes:
  if (p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2)**3 >=n: break
   if p2>p1:
    for i in range(0,int(math.log(n)/math.log(p1))+1):
     for j in range(0,int(math.log(n/p1**i)/math.log(p2))+1):
      sums+=n//(p1**(i+3)*p2**(j+3)) #(p1**3的倍數*p2**3的倍數
 return sums

在迴圈開頭新增退出條件,可以算出10^18內的了

import math
def Sai(n):
 primes=[]
 h  = [True] * int(n**(1.0/3)+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** (1.0/3*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)
 sums=n
 for p in primes:
  for i in range(0,int(math.log(n)/math.log(p))+1):
   sums+=n//(p**(i+3)) #p**3的倍數
 for p1 in primes:
  if (p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2)**3 >=n: break
   if p2>p1:
    for i in range(0,int(math.log(n)/math.log(p1))+1):
     for j in range(0,int(math.log(n/p1**i)/math.log(p2))+1):
      if (p1**(i+3)*p2**(j+3))>n: break
      sums+=n//(p1**(i+3)*p2**(j+3)) #(p1**3的倍數*p2**3的倍數
 return sums

以上算到2個不同質因數的例子說明,依次把1個質數立方、4次方、5次方...的倍數,2個質數立方、4次方、5次方...之積的倍數進行計數,就可以得出不超過2個質數立方因數的個數,而每次只需要新增當前不同質數的個數即可，因為少於這個數的都已經在前面計數過了，依次類推,只要再羅列下去。就可以得到所有最多7個質數立方組成的10^18（最大的小於10⁶的不同質因數之積2*3*5*7*11*13*17=510510）以內的立方滿因數個數。

例如：三個質數之積(2x3x5)^3=27,000,有8個(=2x2x2)立方滿因數1 8 27 125, 8*27 8*125,27*125, 8*27*125，在前面計算中作為8的倍數和27的倍數和125的倍數各算1個和1個和1個，算上都有的1，再減去前面作為2質數之積的3個，唯一遺漏的，也就是8*27*125。

把上述程式碼一直寫下去很枯燥，也沒有必要，因為隨著不同質因數的增加，每個的次方數範圍會很小，質數本身也很小，到7個因數時，最大的質數只能是31。每個的次方數都只能是1。

2*3*5*7*11*13*31=930930，2*3*5*7*11*17*31=1217370

上述的程式碼存在bug,掃描質數表的上限太大了，其實能減少3個，時間一下變成一半。

import math
def Said(n):
 primes=[]
 h  = [True] * int(n**(1.0/3)+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** (1.0/3*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)
 sums=n
 for p in primes:
  for i in range(0,int(math.log(n)/math.log(p))+1-3):
   sums+=n//(p**(i+3)) #p**3的倍數
 for p1 in primes:
  if (p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2)**3 >=n: break
   if p2>p1:
    for i in range(0,int(math.log(n)/math.log(p1))+1-3):
     for j in range(0,int(math.log(n/p1**i)/math.log(p2))+1-3):
      if (p1**(i+3)*p2**(j+3))>n: break
      sums+=n//(p1**(i+3)*p2**(j+3)) #(p1**3的倍數*p2**3的倍數
 return sums

繼續微調2個質數迴圈的上界，p1^(i+3)*p2^(j+3)<=n,因為p2大於p1,所以p1^(i+3)*p1^3<=n,通過計算可得i<=(log(n)/log(p1)-6,同理j<=log(n/p^(i+3))-3。修改後時間再減一點。加上3個不同質數因數的情況，用時一下擴大了幾倍。

 #-------------------------------p1^n1*p2^n2*p3^n3*x
 for p1 in primes:
  if (p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3)**3 >n: break
     if p3>p2:
      for i in range(0,int(math.log(n,p1)-3-3)+1-3):
       for j in range(0,int(math.log(n/p1**(i+3),p2)-3)+1-3):
        for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3))+1-3):
         if (p1**(i+3)*p2**(j+3)*p3**(k+3))>n: break
         sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3))

最後通過的程式碼，實在是醜陋，執行時間倒是在1分鐘以內（使用pypy3 7.2）。

import math
def Sai7g(n):
 primes=[]
 h  = [True] * int(n**(1.0/3)+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** (1.0/3*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)
 sums=n
 #----------------------------------------p^n*x
 for p in primes:
  for i in range(0,int(math.log(n,p))+1-3):
   sums+=n//(p**(i+3)) 
 #-----------------------------------p1^n1*p2^n2*x
 for p1 in primes:
  if (p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2)**3 >n: break
   if p2>p1:
    for i in range(0,int(math.log(n,p1)-3)+1-3):
     for j in range(0,int(math.log(n/p1**(i+3),p2))+1-3):
      sums+=n//(p1**(i+3)*p2**(j+3)) 
 #-------------------------------p1^n1*p2^n2*p3^n3*x
 for p1 in primes:
  if (p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3)**3 >n: break
     if p3>p2:
      for i in range(0,int(math.log(n,p1)-3-3)+1-3):
       for j in range(0,int(math.log(n/p1**(i+3),p2)-3)+1-3):
        for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3))+1-3):
         sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3)) 
 #-----------------------------------p1^n1*p2^n2*p3^n3*p4^n4*x
 for p1 in primes:
  if (p1*p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3*p3)**3 >=n: break
     if p3>p2:
      for p4 in primes:
       if (p1*p2*p3*p4)**3 >n: break
       if p4>p3:
        for i in range(0,int(math.log(n,p1)-3-3-3)+1-3):
         for j in range(0,int(math.log(n/p1**(i+3),p2)-3-3)+1-3):
          for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3)-3)+1-3):
           for a in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3),p4))+1-3):
            sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3)*p4**(a+3)) 
 #-----------------------------------p1^n1*p2^n2*p3^n3*p4^n4*p5^n5*x
 for p1 in primes:
  if (p1*p1*p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3*p3*p3)**3 >=n: break
     if p3>p2:
      for p4 in primes:
       if (p1*p2*p3*p4*p4)**3 >n: break
       if p4>p3:
        for p5 in primes:
         if (p1*p2*p3*p4*p5)**3 >n: break
         if p5>p4:
          for i in range(0,int(math.log(n,p1)-3-3-3-3)+1-3):
           for j in range(0,int(math.log(n/p1**(i+3),p2)-3-3-3)+1-3):
            for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3)-3-3)+1-3):
             for a in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3),p4)-3)+1-3):
              for b in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3),p5))+1-3):
               sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3)*p4**(a+3)*p5**(b+3)) 
 #-----------------------------------p1^n1*p2^n2*p3^n3*p4^n4*p5^n5*p6^n6*x
 for p1 in primes:
  if (p1*p1*p1*p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2*p2*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3*p3*p3*p3)**3 >=n: break
     if p3>p2:
      for p4 in primes:
       if (p1*p2*p3*p4*p4*p4)**3 >n: break
       if p4>p3:
        for p5 in primes:
         if (p1*p2*p3*p4*p5*p5)**3 >n: break
         if p5>p4:
          for p6 in primes:
           if (p1*p2*p3*p4*p5*p6)**3 >n: break
           if p6>p5:
            for i in range(0,int(math.log(n,p1)-3-3-3-3-3)+1-3):
             for j in range(0,int(math.log(n/p1**(i+3),p2)-3-3-3-3)+1-3):
              for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3)-3-3-3)+1-3):
               for a in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3),p4)-3-3)+1-3):
                for b in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3),p5)-3)+1-3):
                 for c in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3)/p5**(b+3),p6))+1-3):
                  sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3)*p4**(a+3)*p5**(b+3)*p6**(c+3)) 
 #-----------------------------------p1^n1*p2^n2*p3^n3*p4^n4*p5^n5*p6^n6*p7^n7*x
 for p1 in primes:
  if (p1*p1*p1*p1*p1*p1*p1)**3 >=n: break
  for p2 in primes:
   if (p1*p2*p2*p2*p2*p2*p2)**3 >=n: break
   if p2>p1:
    for p3 in primes:
     if (p1*p2*p3*p3*p3*p3*p3)**3 >=n: break
     if p3>p2:
      for p4 in primes:
       if (p1*p2*p3*p4*p4*p4*p4)**3 >n: break
       if p4>p3:
        for p5 in primes:
         if (p1*p2*p3*p4*p5*p5*p5)**3 >n: break
         if p5>p4:
          for p6 in primes:
           if (p1*p2*p3*p4*p5*p6*p6)**3 >n: break
           if p6>p5:
            for p7 in primes:
             if (p1*p2*p3*p4*p5*p6*p7)**3 >n: break
             if p7>p6:
              for i in range(0,int(math.log(n,p1)-3-3-3-3-3-3)+1-3):
               for j in range(0,int(math.log(n/p1**(i+3),p2)-3-3-3-3-3)+1-3):
                for k in range(0,int(math.log(n/p1**(i+3)/p2**(j+3),p3)-3-3-3-3)+1-3):
                 for a in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3),p4)-3-3-3)+1-3):
                  for b in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3),p5)-3-3)+1-3):
                   for c in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3)/p5**(b+3),p6)-3)+1-3):
                    for d in range(0,int(math.log(n/p1**(i+3)/p2**(j+3)/p3**(k+3)/p4**(a+3)/p5**(b+3)/p6**(c+3),p7))+1-3):
                     sums+=n//(p1**(i+3)*p2**(j+3)*p3**(k+3)*p4**(a+3)*p5**(b+3)*p6**(c+3)*p7**(d+3)) 
 return sums

尤拉計劃694:立方滿因數

相關文章