趣味编程之计算相亲数(上)

    一直想写这篇关于算法的文章，但是由于看到园子里众多研究算法的高手让我一直没有决心写下来，但高手归高手，不是高手也可以写出来让高手们拍砖，所以今天就在这里献丑了。相亲数和完全数作为数学问题的确是顶级难题，可是拿来做娱乐就不同了，从刚接触编程时C语言书上的课后习题到两年前的Intel多核编程大赛，这个题目一直伴随着我们，让我们来娱乐一下吧。

    简单说一下概念，相亲数是指两个正整数中，彼此的全部约数之和（本身除外）与另一方相等。举例来说：

220的全部约数(除掉本身)相加是：1+2+4+5+10+11+20+22+44+55+110=284

284的全部约数(除掉本身)相加的和是：1+2+4+71+142=220

    所以220和284就是一对相亲数。

    那什么是完全数呢？即它所有的真因子（即除了自身以外的约数）的和恰好等于它本身。例如：

第一个完全数是6，它有约数1、2、3、6，除去它本身6外，其余3个数相加，1+2+3=6

第二个完全数是28，它有约数1、2、4、7、14、28，除去它本身28外，其余5个数相加，1+2+4+7+14=28

    概念不多说了，直接上算法。

    算法一：直接计算约数之和

    这是最直接的方式了，循环计算所有可能成为约数的数字然后加和，直接到不用做过多的解释，只要会写程序，用任何语言都能实现！

  1: /// <summary>
  2: /// 直接计算约数之和(串行)
  3: /// </summary>
  4: public class Algorithm1
  5: {
  6:     private int GetSum(int num)
  7:     {
  8:         int sum = 1;
  9:         int limit = (int)Math.Sqrt(num);
 10:         for (int i = 2; i <= limit; i++)
 11:             if (num % i == 0) sum += i + num / i;
 12:         return sum;
 13:     }
 14: 
 15:     public void Run(int from, int to)
 16:     {
 17:         int perfertCount = 0;
 18:         int amicablePairCount = 0;
 19:         for (int num = from; num <= to; num++)
 20:         {
 21:             int sum1 = this.GetSum(num);
 22:             if (sum1 == num)
 23:             {
 24:                 Console.WriteLine("{0}是完全数", num);
 25:                 perfertCount++;
 26:             }
 27:             if (sum1 > num)
 28:             {
 29:                 int sum2 = this.GetSum(sum1);
 30:                 if (sum2 == num)
 31:                 {
 32:                     Console.WriteLine("{0}和{1}是一对相亲数", sum1, sum2);
 33:                     amicablePairCount++;
 34:                 }
 35:             }
 36:         }
 37:         Console.WriteLine("在{0}到{1}中共有{2}个完全数和{3}对相亲数", from, to, perfertCount, amicablePairCount);
 38:     }
 39: }

    测试代码，从2计算到5000000：

  1: static void Main(string[] args)
  2: {
  3:     var stopwatch = Stopwatch.StartNew();
  4:     Algorithm1 algorithm = new Algorithm1();
  5:     algorithm.Run(2, 5000000);
  6:     stopwatch.Stop();
  7:     Console.WriteLine("计算完成共花费{0}秒", stopwatch.Elapsed.TotalSeconds);
  8:     Console.ReadKey();
  9: }

    在我的ThinkPad R400上测试运行时间大概在51秒左右，速度能不能再提高呢，让我们看看.Net4.0为我们带来的并行计算的新特性表现如何。

  1: /// <summary>
  2: /// 直接计算约数之和(并行)
  3: /// </summary>
  4: public class Algorithm2
  5: {
  6:     private int GetSum(int num)
  7:     {
  8:         int sum = 1;
  9:         int limit = (int)Math.Sqrt(num);
 10:         for (int i = 2; i <= limit; i++)
 11:             if (num % i == 0) sum += i + num / i;
 12:         return sum;
 13:     }
 14: 
 15:     public void Run(int from, int to)
 16:     {
 17:         int perfertCount = 0;
 18:         int amicablePairCount = 0;
 19:         Parallel.For(from, to, num =>
 20:         {
 21:             int sum1 = this.GetSum(num);
 22:             if (sum1 == num)
 23:             {
 24:                 Console.WriteLine("{0}是完全数", num);
 25:                 perfertCount++;
 26:             }
 27:             if (sum1 > num)
 28:             {
 29:                 int sum2 = this.GetSum(sum1);
 30:                 if (sum2 == num)
 31:                 {
 32:                     Console.WriteLine("{0}和{1}是一对相亲数", sum1, sum2);
 33:                     amicablePairCount++;
 34:                 }
 35:             }
 36:         });
 37:         Console.WriteLine("在{0}到{1}中共有{2}个完全数和{3}对相亲数", from, to, perfertCount, amicablePairCount);
 38:     }
 39: }

    注意第19行，我们使用System.Threading.Tasks下的Parallel类取代传统的for循环，由于在该算法中每一次计算都是独立的，所以很适合并行，废话不多说，直接运行看结果，运行时间在26秒左右，由于我的机器是双核，所以同样是从2计算到5000000，并行的时间差不多是之前的(51秒)一半，看来Parallel真是不错的新工具啊！当然，这个是从技术上达到了速度的提升，算法本质还没有变，那能不能从算法本身提高计算效率呢？答案当然是肯定的！

    算法二：通过计算所有质因数来计算约数之和

    先说一下原理：记得小学的奥赛有一个题型是计算一个自然数的约数的个数(现在还有没有不清楚)，截法很简单，先分解质因数，然后把所有质因子的幂加一再相乘就是约数的个数，例如数字36=2^2×3^2，那么36就有(2+1)×(2+1)=9个约数(1,2,3,4,6,9,12,18,36)。其实能算出有9个约数，自然可以计算9个约数是什么，是2个2和2个3排列组合的结果，剩下的就不用我废话了，该算法的思路就是先分解质因数然后在逐个计算约数并加和，上代码：

  1: /// <summary>
  2: /// 通过计算所有质因数来计算约数之和（串行）
  3: /// </summary>
  4: public class Algorithm3
  5: {
  6:     private int GetNextFactor(int num)
  7:     {
  8:         int limit = (int)(Math.Sqrt(num));
  9:         for (int i = 2; i <= limit; i++)
 10:             if (num % i == 0)
 11:                 return i;
 12:         return num;
 13:     }
 14:     private List<int> Decomposition(int num)
 15:     {
 16:         var factors = new List<int>();
 17:         while (true)
 18:         {
 19:             var divisor = GetNextFactor(num);
 20:             factors.Add(divisor);
 21:             if (divisor == num) break;
 22:             num = num / divisor;
 23:         }
 24:         return factors;
 25:     }
 26:     private int Sum(List<int> divisors)
 27:     {
 28:         int sum = 0;
 29:         for (int i = 0, count = divisors.Count - 1; i < count; i++)
 30:             sum += divisors[i];
 31:         return sum;
 32:     }
 33:     private int GetSum(List<int> factors)
 34:     {
 35:         if (factors.Count == 1) return 1;//质数
 36:         var divisors = new List<int>() { 1 };
 37:         var factorPows = new List<List<int>>() { new List<int>() { factors[0] } };
 38:         for (int i = 1, count = factors.Count; i < count; i++)
 39:         {
 40:             var length = factorPows.Count;
 41:             if (factors[i] == factorPows[length - 1][0])
 42:                 factorPows[length - 1].Add(Convert.ToInt32(Math.Pow(Convert.ToDouble(factors[i]), Convert.ToDouble(factorPows[length - 1].Count + 1))));
 43:             else
 44:                 factorPows.Add(new List<int>() { factors[i] });
 45:         }
 46:         for (int f = 0, fCount = factorPows.Count; f < fCount; f++)
 47:             for (int d = 0, dCount = divisors.Count; d < dCount; d++)
 48:                 for (int p = 0, pCount = factorPows[f].Count; p < pCount; p++)
 49:                     divisors.Add(divisors[d] * factorPows[f][p]);
 50:         return Sum(divisors);
 51:     }
 52:     public void Run(int from, int to)
 53:     {
 54:         int perfertCount = 0;
 55:         int amicablePairCount = 0;
 56:         for (var num = from; num < to; num++)
 57:         {
 58:             int sum1 = this.GetSum(Decomposition(num));
 59:             if (sum1 == num)
 60:             {
 61:                 Console.WriteLine("{0}是完全数", num);
 62:                 perfertCount++;
 63:             }
 64:             if (sum1 > num)
 65:             {
 66:                 int sum2 = this.GetSum(Decomposition(sum1));
 67:                 if (sum2 == num)
 68:                 {
 69:                     Console.WriteLine("{0}和{1}是一对相亲数", sum1, sum2);
 70:                     amicablePairCount++;
 71:                 }
 72:             }
 73:         }
 74:         Console.WriteLine("在{0}到{1}中共有{2}个完全数和{3}对相亲数", from, to, perfertCount, amicablePairCount);
 75:     }
 76: }

    先看速度如何，是否比算法一更快，从2计算到5000000花费近27秒，几乎与算法一的并行版本接近，果然快很多，那么快在哪里了呢？算法一做了大量的取模和除法操作，相比于加、减、乘等操作这些操作都是很耗时的，而算法二先计算质因数，减少了很多取模和除法操作，取而代之的是很多乘法和指数运算，性能自然得到提升，但还算细心的我并没有就此下结论，我把计算范围缩小到2到100000，此时，算法一花费时间是0.18秒而算法而则是0.36秒，算法一反而取胜了，其实道理很简单，随着数字的增大，算法一计算取模和除法的操作增加也不断增加，而算法二随着数字增大计算次数却增加不明显，因为分解质因数其实是找质数的过程，但找到第一个质因数后，数字就变小了，计算量并没增加多少，相反，找到的质因数越多，计算约数的计算量就越大，所以在数字不大(相对不大)的领域里，试商次数不多所以算法一很快完成了计算，而算法二相对于算法一的却没什么太大优势，但随着数字的变大，算法二的优势会越来越明显！例如我从5000000计算到5500000，算法一就豪无优势了，落后算法二一半都不止，这让我想起了一个古老的但却众所周知的梵塔问题，算法一就是这样一个梵塔……

    当然我也没忘记把这个算法改成并行版本，我就不贴代码了，只要改第56行的for就可以了，从2计算到5000000花费在16秒左右，这样我们已经从最初的51秒降低到16秒，还能更快吗？我绞尽脑汁暂时还没什么结果，因为我发现最后所花的时间都是在寻找质数上了，难怪数学界的顶级难题都跟质数有关或者围绕质数展开，还有我们程序员熟悉的加密算法RSA，也是基于大整数难以分解的“特性”上，如果不幸被我找到了快速分解算法我就不用在这写博客啦，扯远了，还是回归娱乐吧，我们还有没有办法让它再快点呢？

    算法二SP1：之所以叫SP1是因为它并没有本质上的更改，只是在外围让它显得更快。话说算法界有两大秘籍，九阴真经和九阳神功——都不是，开个玩笑，其实没什么秘籍，而是方法论，无非就是时间换空间和空间换时间，根据我们的需要，时间和空间在我们的代码中转换来转换去，既然我追求的是速度，那自然是用空间来换时间喽。

    算法一我还没想到哪里可以用空间来换取速度，倒是在对算法二的研究过程中我意识到大量的重复计算，最典型的是计算质数，如果缓存这些质数速度会不会快一些呢，其实是废话当然会快，花了空间速度还没提升的事情谁会愿意做呢，但仅仅缓存质数远远不够，因为最大量的计算根本不在这里，如果连续的看待分解的过程，你就会发现它是一个递归的过程，之前的计算结果对后面很有用，比如我们已经分解了36=2^2×3^2。那么当我们分解72的时候是怎样的过程呢，先找到了第一个因子2，然后得到36，继续分解，36的分解过程又做一次，重复量可想而知，有人说，你把2^2×3^2记录下来，在计算72的时候直接利用36的计算结果，说的没错，但我记录的信息是不是太大了，空间也不是这么挥霍的啊，于是我权衡再三，采用了下面的算法，先看代码：

  1: /// <summary>
  2: /// 通过计算所有质因数来计算约数之和（空间算法）
  3: /// </summary>
  4: public class Algorithm5
  5: {
  6:     public List<int> primeList = new List<int>();
  7:     public int[] firstFactorList;
  8:     public int[] remainingList;
  9:     public int[] resultList;
 10:     public int GetNextFactor(int num)
 11:     {
 12:         var max = (int)Math.Sqrt(num);
 13:         for (int i = 0; i < primeList.Count; i++)
 14:         {
 15:             var p = primeList[i];
 16:             if (p > max) break;
 17:             if (num % p == 0)
 18:                 return p;
 19:         }
 20:         primeList.Add(num);
 21:         return num;
 22:     }
 23:     public List<int> Decomposition(int num)
 24:     {
 25:         var divisors = new List<int>();
 26:         var factor = firstFactorList[num] = GetNextFactor(num);
 27:         if (factor == num)
 28:             remainingList[num] = 1;
 29:         else
 30:             remainingList[num] = num / firstFactorList[num];
 31:         while (true)
 32:         {
 33:             divisors.Add(firstFactorList[num]);
 34:             if (remainingList[num] == 1) break;
 35:             num = remainingList[num];
 36:         }
 37:         return divisors;
 38:     }
 39:     private int Sum(List<int> divisors)
 40:     {
 41:         int sum = 0;
 42:         for (int i = 0, count = divisors.Count - 1; i < count; i++)
 43:             sum += divisors[i];
 44:         return sum;
 45:     }
 46:     public int GetSum(List<int> factors)
 47:     {
 48:         if (factors.Count == 1) return 1;//质数
 49:         var divisors = new List<int>() { 1 };
 50:         var factorPows = new List<List<int>>() { new List<int>() { factors[0] } };
 51:         for (int i = 1, count = factors.Count; i < count; i++)
 52:         {
 53:             var length = factorPows.Count;
 54:             if (factors[i] == factorPows[length - 1][0])
 55:                 factorPows[length - 1].Add(Convert.ToInt32(Math.Pow(Convert.ToDouble(factors[i]), Convert.ToDouble(factorPows[length - 1].Count + 1))));
 56:             else
 57:                 factorPows.Add(new List<int>() { factors[i] });
 58:         }
 59:         for (int f = 0, fCount = factorPows.Count; f < fCount; f++)
 60:             for (int d = 0, dCount = divisors.Count; d < dCount; d++)
 61:                 for (int p = 0, pCount = factorPows[f].Count; p < pCount; p++)
 62:                     divisors.Add(divisors[d] * factorPows[f][p]);
 63:         return Sum(divisors);
 64:     }
 65:     public void Run(int limit)
 66:     {
 67:         firstFactorList = new int[limit];
 68:         remainingList = new int[limit];
 69:         resultList = new int[limit];
 70:         int perfertCount = 0;
 71:         int amicablePairCount = 0;
 72:         for (var num = 2; num < limit; num++)
 73:         {
 74:             var result = resultList[num] = this.GetSum(Decomposition(num));
 75:             if (result == num)
 76:             {
 77:                 Console.WriteLine("{0}是完全数", num);
 78:                 perfertCount++;
 79:             }
 80:             else if (result < num && resultList[result] == num)
 81:             {
 82:                 Console.WriteLine("{0}和{1}是一对相亲数", result, num);
 83:                 amicablePairCount++;
 84:             }
 85:         }
 86:         Console.WriteLine("在{0}到{1}中至少有{2}个完全数和{3}对相亲数", 2, limit, perfertCount, amicablePairCount);
 87:     }
 88: }

    我缓存了质数，每个数字的第一个因子和剩余因子的乘积以及每个数字的约数和，代码之所以没有注释是因为我写不清楚，给大家举个例子，大家再对照代码看就好：比如分解72，先找到它的第一个因子2，和剩余因子的乘积36(其实就是72/2)，然后缓存2和36做为72对应的缓存变量，然后在缓存列表中找到36的第一个因子(因为之前已经计算过)，也是2，然后看看36剩余因子的乘积是18有找到了18的第一个因子9……就这样利用了原来的结果，把取模和除法变成了查(缓存)表，这样无疑有个弊端，我们不能从中间开始算了，必须从2开始计算，先不管了，看看速度再说，从2计算到5000000花费了13秒左右，注意这里计算次数跟之前不一样，之前的算法是算到5000000，但可以超过5000000，而现在是只算到5000000，不管怎样，速度比并行版本的算法二还要快一些(前提是从2开始计算)，缓存的效果不错哦，不过内存也被吃去大片，空间换时间的代价……

    算法二SP2：本来以上就是我要记录的全部内容，但由于迟迟未成文，所以最近我又发现一个“超级”快速的空间换时间的算法，比SP1快很多，请允许我先卖一个关子，先公布结果，从2计算到5000000只需要2.8秒，具体实现且听下回分解……(无耻的淫笑)

    后记：郑重声明，本贴纯属娱乐帖，很多代码并没有在细节性能上过多纠缠(比如Math.Sqrt的性能还有提升空间，List申请空间处的性能提升等等……)，另外卖关子不是吊胃口，而是要说明这个SP2算法所花费的篇幅可能比较长，实在写不动了，就卖了一个关子，小弟会尽快完成下一篇，敬请期待！不过通过把这篇文章写出来(很多代码多年前就写好了)还是感触颇深的，我们面对一个问题，可以从多个角度去分析优化解决，可以从根本上想办法(算法本身)，也可以找工具(并行)，还可以在外围和结构等方面想办法(缓存……)，只要我们不停的思考，从多个角度想办法，总会有些手段给我们利用。最后，不知大家是不是也研究过此类问题，或许您有更好的算法，无论是出于兴趣还是其它，欢迎一起讨论and拍砖。

    希望本文对您有帮助(完整代码下篇一起提供)。