题目95：找出元素都不超过一百万的亲和链中最小的元素

欧拉计划

Amicable chains

The proper divisors of a number are all the divisors excluding the number itself. For example, the proper divisors of 28 are 1, 2, 4, 7, and 14. As the sum of these divisors is equal to 28, we call it a perfect number.

Interestingly the sum of the proper divisors of 220 is 284 and the sum of the proper divisors of 284 is 220, forming a chain of two numbers. For this reason, 220 and 284 are called an amicable pair.

Perhaps less well known are longer chains. For example, starting with 12496, we form a chain of five numbers:

12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → ...)

Since this chain returns to its starting point, it is called an amicable chain.

Find the smallest member of the longest amicable chain with no element exceeding one million.

题目：

一个数的真因子是指除了其本身以外的所有因子。例如，28 的真因子是 1,2,4,7 和 14。因为这些因子之和等于 28，我们称 28 为一个完美数。

有趣的是，220 的真因子之和是 284，而 284 的真因子之和是 220，形成了一个两个元素的链。因此 220 和 284 被称为一个亲和对。

可能更长的链就不那么为人所知了，例如，从 12496 开始，我们可以得到一个五个元素的链：

12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → ...)

因为这条链最后回到了开始的数，它被称为一条亲和链。

对于元素都不超过一百万的亲和链，找出最长的亲和链中最小的元素。

愤怒的大头菇

1. import math
   
2. def Hr(x):
   
3.     count =1
   
4.     for i in range(2,int(math.sqrt(x)+1)):
   
5.         if x % i == 0:
   
6.             if i!= x/i:
   
7.                 count += i+x/i
   
8.             else:
   
9.                 count += i
   
10.     return int(count)
    
11. def Prime(x):
    
12.     if x > 1:
    
13.         if x==2:
    
14.             return True
    
15.         if x%2==0:
    
16.             return False
    
17.         for i in range(3,int(math.sqrt(x)+1),2):
    
18.             if x%i==0:
    
19.                 return False
    
20.         return True
    
21.     return False
    
22. i=4
    
23. temp = i
    
24. list1=[]
    
25. list2 = []
    
26. while True:
    
27.     if not Prime(i):
    
28.         list2=[]
    
29.         while temp != 1:
    
30.             temp = Hr(temp)
    
31.             if temp in list2:
    
32.                 break
    
33.             if i not in list2:
    
34.                 list2.append(i)
    
35.             if temp not in list2:
    
36.                 list2.append(temp)
    
37.             if temp > 1000000:
    
38.                 break
    
39.                
    
40.             if temp == i:
    
41.                 break
    
42.         if temp == i and temp<1000000:
    
43.             list1.append(list2)
    
44.     i += 1
    
45.     temp = i
    
46.     if i > 1000000:
    
47.         break
    
48. tmp = 0
    
49. for each in list1:
    
50.     each.sort()
    
51.     if len(each) > tmp:
    
52.         tmp = len(each)
    
53. list3=[]
    
54. for each in list1:
    
55.     if len(each) == tmp:
    
56.         list3.append(each[0])
    
57. list3.sort()
    
58. print(list3[0])

复制代码

效率低得不行，结果：14316

jerryxjr1220

愤怒的大头菇 发表于 2016-11-22 21:17
效率低得不行，结果：14316

我的效率也好低，2分钟才能出结果。。。

1. # -*- coding: utf-8 -*-
   
2. """
   
3. Created on Thu Jan  5 21:21:52 2017
   
4. 
   
5. @author: Jerry Xu
   
6. """
   
7. 
   
8. primes = [True]*1000001
   
9. primes[0], primes[1] = False, False
   
10. for i,prime in enumerate(primes):
    
11.         if prime:
    
12.                 for j in range(i*i,1000001,i):
    
13.                         primes[j] = False
    
14. dic = {}
    
15. chain = {}
    
16. chains = {}
    
17. def ele(n):
    
18.         if primes[n]:
    
19.                 dic[n] = [1]
    
20.         else:
    
21.                 dic[n] = [1]
    
22.                 for i in range(2,int(n**0.5+1)):
    
23.                         if n % i == 0 and i*i != n:
    
24.                                 dic[n] += [i,n//i]
    
25.                         elif i*i == n:
    
26.                                 dic[n] += [i]
    
27.         chain[n] = sum(dic[n])
    
28. def find(x):
    
29.         tmp = [x]
    
30.         y = chain[x]
    
31.         while True:
    
32.                 if y > 1000000:
    
33.                         tmp = []
    
34.                         chains[y] = 0
    
35.                         break
    
36.                 else:
    
37.                         if y in tmp:
    
38.                                 index = tmp.index(y)
    
39.                                 if len(tmp)-index > chains.get(y,0):
    
40.                                         chains[y] = len(tmp)-index
    
41.                                 break
    
42.                         else:
    
43.                                 tmp.append(y)
    
44.                                 y = chain[y]
    
45.         return y,chains[y]
    
46. for i in range(1,1000001):
    
47.         ele(i)
    
48. maxlength = 0
    
49. for i in range(1,1000001):
    
50.         if find(i)[1] > maxlength:
    
51.                 maxlength = find(i)[1]
    
52.                 maxi = find(i)[0]
    
53. k = maxi
    
54. longest = []
    
55. for i in range(maxlength):
    
56.         k = chain[k]
    
57.         longest.append(k)
    
58. print (min(longest))

复制代码

结果和你是一样的：
14316
[Finished in 123.9s]

芒果加黄桃

1. # encoding:utf-8
   
2. # 找出不超过100万的亲和链中的最小元素
   
3. from time import time
   
4. N = 1000000
   
5. l_sums, l_lens = [1] * (2 * N), [0] * (2 * N)
   
6. def getFacSum():
   
7.     for i in range(2, N // 2 + 1):
   
8.         for t in range(2 * i, N + 1, i):
   
9.             l_sums[t] += i
   
10. def euler095():
    
11.     getFacSum()
    
12.     max_len_nums = set()
    
13.     for i in range (1, N + 1):
    
14.         s_t , t, flag = list(), i, False
    
15.         s_t.append(t)
    
16.         while True:
    
17.             if l_sums[t] > N + 1:
    
18.                 break
    
19.             if l_sums[t] not in s_t:
    
20.                 s_t.append(l_sums[t])
    
21.                 t = l_sums[t]
    
22.             else:
    
23.                 s_t = s_t[s_t.index(l_sums[t]):]
    
24.                 flag = True
    
25.                 break
    
26.         if flag:
    
27.             l_lens[s_t[0]] = len(s_t)
    
28.             if len(s_t) > len(max_len_nums):
    
29.                 max_len_nums = s_t.copy()
    
30.     print(max_len_nums)
    
31.     print(min(max_len_nums))
    
32. if __name__ == '__main__':
    
33.     start = time()
    
34.     euler095()
    
35.     print('cost %.6f sec' % (time() - start))
    

复制代码

[14316, 19116, 31704, 47616, 83328, 177792, 295488, 629072, 589786, 294896, 358336, 418904, 366556, 274924, 275444, 243760, 376736, 381028, 285778, 152990, 122410, 97946, 48976, 45946, 22976, 22744, 19916, 17716]
14316
cost 16.423642 sec

fc1735

1. bound = int(1e6)
   
2. 
   
3. count = [0] * bound
   
4. for i in range(1, bound // 2):
   
5.   for j in range(2 * i, bound, i):
   
6.     count[j] += i
   
7. 
   
8. flag = [0] * bound
   
9. def find(n):
   
10.   if flag[n] == 0:
    
11.     s = set()
    
12.     p = n
    
13.     while p not in s:
    
14.       s.add(p)
    
15.       flag[p] = -1
    
16.       p = count[p]
    
17.       if p > bound:
    
18.         return -1
    
19.     l = list()
    
20.     while p in s:
    
21.       l.append(p)
    
22.       s.remove(p)
    
23.       p = count[p]
    
24.     for i in l:
    
25.       flag[i] = len(l)
    
26.   return flag[n]
    
27. 
    
28. print(max(range(1, bound), key = find))

复制代码

https://github.com/devinizz/project_euler/blob/master/page02/95.py

永恒的蓝色梦想

C++ 151ms

1. #include<iostream>
   
2. #include<vector>
   
3. #include<unordered_map>
   
4. 
   
5. 
   
6. 
   
7. int main() {
   
8.     using namespace std;
   
9.     ios::sync_with_stdio(false);
   
10. 
    
11.     constexpr static unsigned int UPPER_BOUND = 1000000;
    
12.     static unsigned int prime_divisor_sum[UPPER_BOUND + 1];
    
13.     unordered_map<unsigned int, unsigned int> indecies;
    
14.     vector<unsigned int> chain;
    
15.     unsigned int i, copy_i, temp, min_of_longest = 0, max_length = 0;
    
16.     unsigned int & index = copy_i, & length = temp;//别名
    
17. 
    
18. 
    
19.     for (i = 1; i <= UPPER_BOUND / 2; i++) {
    
20.         for (temp = i << 1; temp <= UPPER_BOUND; temp += i) {
    
21.             prime_divisor_sum[temp] += i;
    
22.         }
    
23.     }
    
24. 
    
25. 
    
26.     for (i = 1; i <= UPPER_BOUND; i++) {
    
27.         if (prime_divisor_sum[i]) {
    
28.             chain.push_back(i);
    
29.             indecies[i] = 1;
    
30.             copy_i = prime_divisor_sum[i];
    
31.             prime_divisor_sum[i] = 0;
    
32. 
    
33.             for (;;) {
    
34.                 if (copy_i > UPPER_BOUND) {
    
35.                     break;
    
36.                 }
    
37.                 else if (prime_divisor_sum[copy_i]) {
    
38.                     chain.push_back(copy_i);
    
39.                     indecies[copy_i] = chain.size();
    
40.                     temp = prime_divisor_sum[copy_i];
    
41.                     prime_divisor_sum[copy_i] = 0;
    
42.                     copy_i = temp;
    
43.                 }
    
44.                 else {
    
45.                     if (index = indecies[copy_i]) {
    
46.                         length = chain.size() - index;//比真的长度少1
    
47. 
    
48.                         if (length > max_length) {
    
49.                             max_length = length;
    
50. 
    
51.                             for (min_of_longest = chain[index - 1]; index < chain.size(); index++) {
    
52.                                 if (chain[index] < min_of_longest) {
    
53.                                     min_of_longest = chain[index];
    
54.                                 }
    
55.                             }
    
56.                         }
    
57.                     }
    
58. 
    
59.                     break;
    
60.                 }
    
61.             }
    
62. 
    
63.             indecies.clear();
    
64.             chain.clear();
    
65.         }
    
66.     }
    
67. 
    
68. 
    
69.     cout << min_of_longest << endl;
    
70.     return 0;
    
71. }

复制代码

guosl

1. /*
   
2. 答案：14316
   
3. 耗时：2.58101秒（单核）
   
4.       0.287085秒（八核）
   
5. */
   
6. #include <iostream>
   
7. #include <algorithm>
   
8. #include <vector>
   
9. #include <cmath>
   
10. #include <omp.h>
    
11. using namespace std;
    
12. 
    
13. const int N = 1000000;
    
14. int a[N + 1];
    
15. 
    
16. int getAmi(int n)
    
17. {
    
18.   int r = 1;
    
19.   for (int i = 2; i <= (int)sqrt((double)n); ++i)
    
20.   {
    
21.     if (n%i == 0)
    
22.     {
    
23.       int m = n / i;
    
24.       if (m != i)
    
25.         r += m + i;
    
26.       else
    
27.         r += i;
    
28.     }
    
29.   }
    
30.   if (r <= N)
    
31.     return r;
    
32.   else
    
33.     return -1;
    
34. }
    
35. 
    
36. int main(void)
    
37. {
    
38.   double t = omp_get_wtime();
    
39.   int ps = omp_get_num_procs();
    
40.   int *nMaxLens = new int[ps], *nMinElements = new int[ps];
    
41. #pragma omp parallel shared(nMaxLens,nMinElements)
    
42.   {
    
43. #pragma omp for
    
44.     for (int i = 0; i < ps; ++i)
    
45.     {
    
46.       nMaxLens[i] = 0;
    
47.       nMinElements[i] = N;
    
48.     }
    
49. #pragma omp for schedule(guided,4)
    
50.     for (int i = 2; i <= N; ++i)
    
51.       a[i] = getAmi(i);
    
52. #pragma omp for schedule(dynamic)
    
53.     for (int i = 2; i <= N; ++i)
    
54.     {
    
55.       int pt = omp_get_thread_num();
    
56.       vector<int> v;
    
57.       int k = i;
    
58.       bool bFind = false;
    
59.       while (true)
    
60.       {
    
61.         v.push_back(k);
    
62.         k = a[k];
    
63.         if (k <= 1)
    
64.           break;
    
65.         auto itr = find(v.begin(), v.end(), k);
    
66.         if (itr != v.end())
    
67.         {
    
68.           int nc = int(v.end() - itr);
    
69.           if (nc > nMaxLens[pt])
    
70.           {
    
71.             nMaxLens[pt] = nc;
    
72.             nMinElements[pt] = *min_element(itr, v.end());
    
73.           }
    
74.           break;
    
75.         }
    
76.       }
    
77.     }
    
78.   }
    
79.   for (int i = 1; i < ps; ++i)
    
80.   {
    
81.     if (nMaxLens[i] > nMaxLens[0])
    
82.     {
    
83.       nMaxLens[0] = nMaxLens[i];
    
84.       nMinElements[0] = nMinElements[i];
    
85.     }
    
86.   }
    
87.   t = omp_get_wtime() - t;
    
88.   cout << nMinElements[0] << endl << t << endl;;
    
89.   delete[] nMaxLens;
    
90.   delete[] nMinElements;
    
91.   return 0;
    
92. }

复制代码