1. 程式人生 > >程式設計師程式設計藝術:第六章、求解500萬以內的親和數

程式設計師程式設計藝術:第六章、求解500萬以內的親和數

   第六章、親和數問題--求解500萬以內的親和數


前奏
    本章陸續開始,除了繼續保持原有的字串、陣列等面試題之外,會有意識的間斷性節選一些有關數字趣味小而巧的面試題目,重在突出思路的“巧”,和“妙”。本章親和數問題之關鍵字,“500萬”,“線性複雜度”。

第一節、親和數問題
題目描述:
求500萬以內的所有親和數
如果兩個數a和b,a的所有真因數之和等於b,b的所有真因數之和等於a,則稱a,b是一對親和數。
例如220和284,1184和1210,2620和2924。

分析:
    首先得明確到底是什麼是親和數?

親和數問題最早是由畢達哥拉斯學派發現和研究的。他們在研究數字的規律的時候發現有以下性質特點的兩個數:
220的真因子是:1、2、4、5、10、11、20、22、44、55、110;
284的真因子是:1、2、4、71、142。
而這兩個數恰恰等於對方的真因子各自加起來的和(sum[i]表示數i 的各個真因子的和),即
220=1+2+4+71+142=sum[284],
284=1+2+4+5+10+11+20+22+44+55+110=sum[220]。
得284的真因子之和sum[284]=220,且220的真因子之和sum[220]=284,即有sum[220]=sum[sum[284]]=284。

如此,是否已看出絲毫端倪?

如上所示,考慮到1是每個整數的因子,把出去整數本身之外的所有因子叫做這個數的“真因子”。如果兩個整數,其中每一個真因子的和都恰好等於另一個數,那麼這兩個數,就構成一對“親和數”(有關親和數的更多討論,可參考這:http://t.cn/hesH09)。

求解:
    瞭解了什麼是親和數,接下來咱們一步一步來解決上面提出的問題(以下內容大部引自水的原話,同時水哥有一句原話,“在你真正弄弄懂這個範例之前,你不配說你懂資料結構和演算法”)。

  1. 看到這個問題後,第一想法是什麼?模擬搜尋+剪枝?回溯?時間複雜度有多大?其中bn為an的偽親和數,即bn是an的真因數之和大約是多少?至少是10^13(@iicup:N^1.5 對於5*10^6 , 次數大致 10^10 而不是 10^13.)的數量級的。那麼對於每秒千萬次運算的計算機來說,大概在1000多天也就是3年內就可以搞定了(iicup的計算: 10^13 / 10^7 =1000000(秒) 大約 278 小時. )。如果是基於這個基數在優化,你無法在一天內得到結果的。
  2. 一個不錯的演算法應該在半小時之內搞定這個問題,當然這樣的演算法有很多。節約時間的做法是可以生成伴隨陣列,也就是空間換時間,但是那樣,空間代價太大,因為資料規模龐大。
  3. 在稍後的演算法中,依然使用的伴隨陣列,只不過,因為題目的特殊性,只是它方便和巧妙地利用了下標作為伴隨陣列,來節約時間。同時,將回溯的思想換成遞推的思想(預處理陣列的時間複雜度為logN(調和級數)*N,掃描陣列的時間複雜度為線性O(N)。所以,總的時間複雜度為O(N*logN+N)(其中logN為調和級數)  )。


第二節、伴隨陣列線性遍歷
依據上文中的第3點思路,編寫如下程式碼:

  1. //求解親和數問題
  2. //第一個for和第二個for迴圈是logn(調和級數)*N次遍歷,第三個for迴圈掃描O(N)。
  3. //所以總的時間複雜度為 O(n*logn)+O(n)=O(N*logN)(其中logN為調和級數)。
  4. //關於第一個for和第二個for尋找中,調和級數的說明:
  5. //比如給2的倍數加2,那麼應該是  n/2次,3的倍數加3 應該是 n/3次,...
  6. //那麼其實就是n*(1+1/2+1/3+1/4+...1/(n/2))=n*(調和級數)=n*logn。
  7. //[email protected] 上善若水
  8. //July、updated,2011.05.24。
  9. #include<stdio.h>
  10. int sum[5000010];   //為防越界
  11. int main()   
  12. {  
  13.     int i, j;  
  14.     for (i = 0; i <= 5000000; i++)   
  15.         sum[i] = 1;  //1是所有數的真因數所以全部置1
  16.     for (i = 2; i + i <= 5000000; i++)  //預處理,預處理是logN(調和級數)*N。
  17.         //@litaoye:調和級數1/2 + 1/3 + 1/4......的和近似為ln(n),
  18.         //因此O(n *(1/2 + 1/3 + 1/4......)) = O(n * ln(n)) = O(N*log(N))。
  19.     {    
  20.         //5000000以下最大的真因數是不超過它的一半的
  21.         j = i + i;  //因為真因數,所以不能算本身,所以從它的2倍開始
  22.         while (j <= 5000000)   
  23.         {    
  24.             //將所有i的倍數的位置上加i
  25.             sum[j] += i;    
  26.             j += i;       
  27.         }  
  28.     }  
  29.     for (i = 220; i <= 5000000; i++)   //掃描,O(N)。
  30.     {  
  31.         // 一次遍歷,因為知道最小是220和284因此從220開始
  32.         if (sum[i] > i && sum[i] <= 5000000 && sum[sum[i]] == i)  
  33.         {  
  34.             //去重,不越界,滿足親和
  35.             printf("%d %d/n",i,sum[i]);  
  36.         }  
  37.     }  
  38.     return 0;  
  39. }  

執行結果:

    @上善若水:

    1、可能大家理解的還不是很清晰,我們建立一個5 000 000 的陣列,從1到2 500 000 開始,在每一個下標是i的倍數的位置上加上i,那麼在迴圈結束之後,我們得到的是什麼?是 類似埃斯托拉晒求素數的陣列(當然裡面有真的親和數),然後只需要一次遍歷就可以輕鬆找到所有的親和數了。時間複雜度,線性。

    2、我們可以清晰的發現連續資料的對映可以通過陣列結構本身的特點替代,用來節約空間,這是資料結構的藝術。在大規模連續資料的回溯處理上,可以通過轉化為遞推生成的方法,逆向思維操作,這是演算法的藝術。

    3、把最簡單的東西運用的最巧妙的人,要比用複雜方法解決複雜問題的人要頭腦清晰。


第三節、程式的構造與解釋
    我再來具體解釋下上述程式的原理,ok,舉個例子,假設是求10以內的親和數,求解步驟如下:

因為所有數的真因數都包含1,所以,先在各個數的下方全部置1

  1. 然後取i=2,3,4,5(i<=10/2),j依次對應的位置為j=(4、6、8、10),(6、9),(8),(10)各數所對應的位置。
  2. 依據j所找到的位置,在j所指的各個數的下面加上各個真因子i(i=2、3、4、5)。
    整個過程,即如下圖所示(如sum[6]=1+2+3=6,sum[10]=1+2+5=8.):
    1  2  3  4  5  6  7  8  9  10
    1  1  1  1  1  1  1  1  1  1
               2      2      2      2
                       3          3 
                               4
                                       5
  3. 然後一次遍歷i從220開始到5000000,i每遍歷一個數後,
    將i對應的數下面的各個真因子加起來得到一個和sum[i],如果這個和sum[i]==某個i’,且sum[i‘]=i,
    那麼這兩個數i和i’,即為一對親和數。
  4. i=2;sum[4]+=2,sum[6]+=2,sum[8]+=2,sum[10]+=2,sum[12]+=2...
    i=3,sum[6]+=3,sum[9]+=3...
    ......
  5. i=220時,sum[220]=284,i=284時,sum[284]=220;即sum[220]=sum[sum[284]]=284,
    得出220與284是一對親和數。所以,最終輸出220、284,...

特別鳴謝

  1. //求解親和數問題
  2. //[email protected] litaoye
  3. //July、胡濱,updated,2011.05.26。
  4. using System;  
  5. using System.Collections.Generic;  
  6. namespace CSharpTest  
  7. {  
  8.     class Program  
  9.     {  
  10.         publicstaticvoid Main()  
  11.         {  
  12.             int max = 5000000;  
  13.             DateTime start = DateTime.Now;  
  14.             int[] counter = CreateCounter(max);  
  15.             for (int i = 0; i < counter.Length; i++)  
  16.             {  
  17.                 int num = counter[i] - i;  
  18.                 //if (num < counter.Length && num > i && counter[num] == counter[i])
  19.                 // Console.WriteLine("{0} {1}", i, num);
  20.             }  
  21.             Console.WriteLine((DateTime.Now - start).TotalSeconds);  
  22.             Console.ReadKey();  
  23.         }  
  24.         staticint[] CreateCounter(int n)  
  25.         {  
  26.             List<int> primes = new List<int>();  
  27.             int[] counter = newint[n + 1];  
  28.             counter[1] = 1;  
  29.             for (int i = 2; i <= n; i++)  
  30.             {  
  31.                 if (counter[i] == 0)  
  32.                 {  
  33.                     counter[i] = i + 1;  
  34.                     primes.Add(i);  
  35.                 }  
  36.                 for (int j = 0; j < primes.Count; j++)  
  37.                 {  
  38.                     if (primes[j] * i > n)   
  39.                         break;  
  40.                     if (i % primes[j] == 0)  
  41.                     {  
  42.                         int k = i;  
  43.                         int l = primes[j] * primes[j];  
  44.                         while (k % primes[j] == 0)  
  45.                         {  
  46.                             l *= primes[j];  
  47.                             k /= primes[j];  
  48.                         }  
  49.                         counter[primes[j] * i] = counter[k] * (l - 1) / (primes[j] - 1);  
  50.                         break;  
  51.                     }  
  52.                     else
  53.                         counter[primes[j] * i] = counter[i] * (primes[j] + 1);  
  54.                 }  
  55.             }  
  56.             return counter;  
  57.         }  
  58.     }  
  59. }  
  60. /* 
  61. 測試結果: 
  62. 0.484375   
  63. 0.484375 
  64. 0.46875 
  65. 單位second。 
  66. */

本章完。

版權所有,本人對本blog內所有任何內容享有版權及著作權。實要轉載,請以連結形式註明出處。