Superabundant Numbers

[medie2005]

一个有意思的问题，也是Knuth出的。

The abundance of an integer n is the sum of the divisors of n (including n itself), divided
by n. Integer n is k-abundant if its abundance is at least k.

For example, the sum of the divisors of 6 is 6+3+2+1=12, and 12/6=2, so 6 is 2-abundant.
As another example, the sum of the divisors of 120 is

120+60+40+30+24+20+15+12+10+8+6+5+4+3+2+1=360

so 120 is 3-abundant. It happens that 6 is the smallest 2-abundant number and 120 is the smallest
3-abundant number. They happen to be exactly 2- and 3-abundant, respectively, but it is generally
possible that the smallest k-abundant number has abundance greater than k.

The task is to write a program that finds the smallest k-abundant number for k=1,2,..
How high can you go?

[mathe]

你是用贪心算法，不能保证全局最优。
前面一次升级的选择对后面的选择是有影响的。
不过贪心算法在本题中还是非常有用的，可以用来先找到一个次优解。然后可以对搜索空间进行更好的裁减。

[mathe]

证明了一个定理，这下应该有很多优化可做了：
i)greedyBound函数里面的二重循环可以改成一重循环了。
ii)find函数里面最后一个循环可以去掉了。
我简单将find里面的后面一个循环去掉，速度就可以提高几十倍了。
看来这个方法改善下去做到20多问题不大。只是现在都是采用浮点数运算，计算精度是个问题。误差会导致结果有偏差

[liangbch]

mark

[zgg___]

窃以为在这个特定问题中，贪心算法有时是可以得到全局最优解。
也就是说：如果37层的方法求出某个整数n1对应于某个k1（k不一定是整数），那么是不存在某个小于n1的整数n2，其所对应的k2大于k1的反例的。感觉上是这样的。
另一方面，可能由于求出的k1略大于题目给出的自然数k，也就是说可能存在小于n1的某个整数n0，它所对应的k0大小位于目标k和k1之间，这样n0才是满足要求的，而n1却不是了。这种情况对应于帖子中的第2个问题，帖子中仅仅说是放宽搜索范围，并没有解决问题的具体方法。

[mathe]

我先写了个程序，估计对于给定的K,最优解最大素因子p(t)的范围
得出
K=2.000000, 2<=t<=5, start searching from t==2
K=3.000000, 3<=t<=5, start searching from t==4
K=4.000000, 4<=t<=5, start searching from t==4
K=5.000000, 6<=t<=9, start searching from t==8
K=6.000000, 9<=t<=15, start searching from t==14
K=7.000000, 14<=t<=20, start searching from t==19
K=8.000000, 22<=t<=28, start searching from t==28
K=9.000000, 35<=t<=46, start searching from t==45
K=10.000000, 55<=t<=69, start searching from t==69
K=11.000000, 89<=t<=107, start searching from t==107
K=12.000000, 142<=t<=165, start searching from t==165
K=13.000000, 230<=t<=260, start searching from t==260
K=14.000000, 373<=t<=407, start searching from t==407
K=15.000000, 609<=t<=651, start searching from t==651
K=16.000000, 996<=t<=1053, start searching from t==1053
K=17.000000, 1637<=t<=1707, start searching from t==1707
K=18.000000, 2698<=t<=2790, start searching from t==2790
K=19.000000, 4461<=t<=4577, start searching from t==4577
K=20.000000, 7398<=t<=7544, start searching from t==7544
K=21.000000, 12301<=t<=12301, start searching from t==12301
K=22.000000, 20503<=t<=20503, start searching from t==20503
K=23.000000, 34253<=t<=34254, start searching from t==34254
K=24.000000, 57348<=t<=57348, start searching from t==57348
K=25.000000, 96198<=t<=96198, start searching from t==96198
K=26.000000, 161659<=t<=161659, start searching from t==161659
K=27.000000, 272124<=t<=272125, start searching from t==272125
K=28.000000, 458789<=t<=458789, start searching from t==458789
K=29.000000, 774616<=t<=774616, start searching from t==774616
好像对于比较大的素数，最后一个值是用一种贪心算法找出最大可能的t,不过好像找的不是很好。

[mathe]

关于这个问题，我现在描写了部分数学结论，放在一个word文件中，可以去
http://download.csdn.net/source/328992
下载。

[mathe]

提供了一个更新的文档：
http://download.csdn.net/source/329898
应该可以更具里面的结论写出一个不错的代码了

[mathe]

我还是用上面的程序，在我的计算机上经过长时间计算，算出了K=19和K=20的情况（不过由于采用浮点计算，结果无法保证）
KK=19
p1^6 p2 p3^3 p4 p5 p6 p9 p16 p4485
P[1]^16 P[2]^10 P[3]^9 P[4]^6 P[5]^5 P[6-8]^4 P[9-15]^3 P[16-61]^2 P[62-4485]
19.000029

real 48m46.056s
user 48m44.560s
sys 0m0.040s

KK=20
p1^8 p2^4 p4 p5 p6 p9 p18 p7429
P[1]^18 P[2-3]^10 P[4]^6 P[5]^5 P[6-8]^4 P[9-17]^3 P[18-76]^2 P[77-7429]
20.000003
real 215m46.801s
user 215m38.650s
sys 0m0.060s
所以如果花些时间，将算法优化一下，计算到K=24是非常可能的。不过更大的估计就很难了，毕竟用我现在的计算结果来估计，K每增加1，大概要花费4倍的时间

[medie2005]

mahte,我怎么没法下载你那个更新的文档啊？

[mathe]

不过写文档时，又发现了一些过去的错误，比如上面的程序中
greedyInit中
if(2*the_prime[i]*the_prime[i]<the_prime[t]){
可以改为：
if(the_prime[i]*the_prime[i]<the_prime[t]){
这个也可以提高一些速度

[mathe]

晕，我要重新上传说已经上传过了，不能上传了，我还是上载到gxqcn的BBS上吧

[mathe]

http://bbs.emath.ac.cn/viewthread.php?tid=34&extra=page%3D1&frombbs=1

[medie2005]

还是没怎么懂，呵呵，
没办法，大学没怎么学数学啊，我是全凭高中的数学底子。

[mathe]

那没办法了。其实除了拉格朗日极值法和琴生不等式，其实都是中学内容差不多可以覆盖了

[medie2005]

拉格朗日极值法和琴生不等式，这些我还是会的。但是，我反而不明白整个算法的流程了。

[mathe]

其实其中使用拉格朗日极值法计算极值的部分也可以改成采用琴生不等式计算。
如文件中定义：
那没办法了。其实除了拉格朗日极值法和琴生不等式，其实都是中学内容差不多可以覆盖了
g(x)=log(1-exp(x))
那么这个函数的二阶导数小于0，所以函数是上凸函数（或者说凹函数），
二阶导数小于0的几何意义就是这个函数在x点的切线随着x增大斜率变小，
上凸函数就是说，这个函数就向一个开口向下的抛物线（比如正弦函数0到Pi之间部分等等）
琴生不等式是说，对于任意x1,x2,...,xn
g((x1+...+xn)/n)>=(g(x1)+...+g(xn))/n,这个从几何意义上面应该比较容易看出来

[medie2005]

虽然我还是不太懂，但是还是谢谢mathe。
其实能得到k=19和k=20，我已经很满足了。（十分感谢mathe花费宝贵的4小时来算这些数）

[mathe]

我在gxqcn的bbs上贴了个更新版本，那个版本不使用拉格朗日极值法了。
而琴生不等式几何意义还是很明显的。

[mathe]

呵呵，计算机花费的4个小时不算什么，不过是机器而已。倒是我自己花费了很多时间去推导设计这个问题。
不过我比较喜欢帮助追求上进的小伙子们:)

[medie2005]

mathe真是太好了，这么帮助我这种菜鸟，感动...