Superabundant Numbers
一个有意思的问题，也是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?
[解决办法]
我觉得可以用排列组合堆素数的办法来找
2^n差一点就是个2-abundant, 3^n则差多一些,考虑同样的数量级上用其他质数来放大
另外,已经找到的最小k-abundance可以参考作为k+1-abundance的步长,因为其间的数都不足够大
[解决办法]
对于整数
p1^a1*...*pt^at
其最小abundance为（不一定是整数）：
- Assembly code
p1-p1^(-a1) pt - pt^(-at)-----------*...* -------------- p1-1 pt - 1
[解决办法]
结果就变成求
Min{a1*log(p1)+a2*log(p2)+...+at*log(pt)}
在给定条件
- Assembly code
p1-p1^(-a1) pt - pt^(-at)-----------*...* -------------->=k p1-1 pt - 1
[解决办法]
这个题目用整数规划做挺有意思
可以首先对充分大的t解方程
[code=Assemble]
p1-p1^(-a1) pt - pt^(-at)
-----------*...* -------------->=k
p1-1 pt - 1
[/code]
得到一组浮点解(a1',a2',...,at')
这样，我们可以得到一个最小值一个下界
a1'*log(p1)+a2'*log(p2)+...+at'*log(pt)
然后对于其中一个浮点数ai,将整个解空间拆分成（先拆分小的下标）
ai<=[ai'] 和ai>=[ai']+1两种情况
分别在两种情况下面求级值，直到找到一组整数解。
不过这道题目数据本身挺规律，应该还有更强的结论
[解决办法]
为什么6是 smallest 2-abundant,而不是3或2？
2+2= 4, 4/2=2
3+2+1 = 6, 6/3 = 2
[解决办法]
做些简单分析，
假设结果是
p1^a1*p2^a2*...*pt^at
其中a1,a2,...,at是正整数,p1,p2,...,pt是素数而且p1<p2<...<pt
我们可以容易得出一下两个结论
i)p1,p2,...,pt正好是前t个素数
ii)a1>=a2>=...>=at
[解决办法]
接下去提供一个非常有用的结论:
对于整数X,我们记abd(X)为它所有正因子的和除以X,也就是X的abundant值(注意这个是一个浮点数)。
假设X0是使得abd(X)=M的最小整数X,那么除了X0=4以外，对于所有的其他这样的X0,其最大的素因子指数必然是1。
比如对于M=7/4,X0=4. 对于M=2,X0=6=2*3,最大素因子3的指数为1，对于M=3,X0=120=8*3*5,最大素因子5的指数也是1.
[解决办法]
现在已经可以给出一种解法了，只是的最大素数的估计范围太大了，如果能够改善一下会好很多。
首先对于方程：
[code=Assemble]
p1-p1^(-a1) pt - pt^(-at)
-----------*...* -------------- = K
p1-1 pt - 1
[/code]
求非负实数(a1,a2,...,at)使得
a1*log(p1)+a2*log(p2)+...+at*log(pt)
的解很简单
令W=power(K*(1-1/p1)(1-1/p2)*...*(1-1/pt),1/t)
其中W不许大于1，不然无解。
ai=-log(1-W)/log(pi)-1即可
也就是取级值时，pi^(ai+1)对于所有的i相等。
对于给定的k,我们需要估计t的上界和下界。
t的下界很简单，因为我们要求K*(1-1/p1)(1-1/p2)...(1-1/pt)<1。通过这个方法可以得到t的下界
比如对于K=2,2*(1-1/2)=1不满足小于1，所以t>=2 (至少有两个素数2,3)
对于K=3, 3*(1-1/2)(1-1/3)=1不满足小于1，所以t>=3(至少有三个素数2,3,5)
而对于t的上界，我估计的不是很好
我们可以让a1=a2=...=at=1得到t的上界，也就是
(1+1/p1)(1+1/p2)...(1+1/pt)>=k的第一个数字t可以作为t的上界
比如对于k=2,(1+1/2)(1+1/3)=2,所以对于k=2，最多指需要前面两个素数
对于k=3,(1+1/2)(1+1/3)(1+1/5)(1+1/7)(1+1/11)(1+1/13)>3,所以对于k=3,最多需要6个素数
可以看出这种方法估计t的上界增长很快。我一直在考虑如何降低这个t的上界的问题。
[解决办法]
有了t的上界以后
我们可以对于这个上界先求出一个浮点解(a1,a2,...,at)
然后我们可以通过从a1开始，每个数字取其上取整，得到整数列(b1,b2,...,bt).然后我们去掉后面一部分数据，使得前面部分对应abundant数还是大于k.这时我们可以得到一个
次优整数解(b1,b2,...,bh)并且得到最终结果的一个上界B=b1*log(p1)+b2*log(p2)+...+bh*log(ph)
然后接下去我们就可以从(a1,a2,...,at)放入栈中采用整数规划的方法求解：
每一次
我们从栈中取一组浮点解(a1',a2',...,at')中找到最后一个不是整数的下标a(s)'(也就是a(s+1)',...,a(t)'都是整数)
我们接下去搜索下面的问题，分别设a(s)=a(s+1)'+1,...,[a(s)']+1.
而a(s+1)=a(s+1)',a(s+2)=a(s+2)',...,a(t)=a(t)'，求最优的浮点解对应的(a(1),a(2),...,a(s-1))
这个浮点问题同样可以用上一楼的方法去计算（如果s==1时，我们改为直接计算对应的最优整数解)
这时对于这个对应的浮点解都可以得到一个浮点的a1*log(p1)+a2*log(p2)+...+at*log(pt)的下界。如果这个值大于B,我们知道这个分支没有最优整数解，裁减；
不然我们将这个浮点解放入栈中。
而对应到某个s==1时（这时得到一组新的整数解），如果这个解得到一个更小的B,那么更新B.
对上面的搜索过程要采用深度优先的方法。
[解决办法]
现在拿k=3作为例子
首先，取6个素数p1=2,p2=3,p3=5,p4=7,p5=11,p6=13
可以解得浮点解:
A[0]=2.506464
A[1]=1.212332
A[2]=0.510152
A[3]=0.249028
A[4]=0.013595
A[5]=-0.052420
由于A[5]<0,这个说明总是只能取A[5]=0
同时我们可以得到一个A[0]=3,A[1]=2,A[2]=1,A[3]=A[4]=A[5]=0时对应的次优解(<K,M>=<3.250000,5.886104>), B=5.886104
list={A[5]=0}
然后对于A[5]=0,可以解得:
A[0]=2.470868
A[1]=1.189874
A[2]=0.494822
A[3]=0.236348
A[4]=0.003306,<K,M>=<3.000000,4.284111> 4.284111<5.886104
list={A[4]=0,A[5]=0; A[4]=1,A[5]=0}
然后对于A[4]=0,A[5]=0,可以解得：
A[0]=2.473743
A[1]=1.191688
A[2]=0.496060
A[3]=0.237372,<K,M>=<3.000000,4.284154>
list={A[3]=0,A[4]=0,A[5]=0;A[3]=1,A[4]=0,A[5]=0;A[4]=1,A[5]=0}
解:A[3]=0,A[4]=0,A[5]=0
A[0]=2.802241
A[1]=1.398947
A[2]=0.637536, <K,M>=<3.000000,4.505340>
list={A[2]=0,A[3]=0,A[4]=0,A[5]=0;A[2]=1,A[3]=0,A[4]=0,A[5]=0;A[3]=1,A[4]=0,A[5]=0;A[4]=1,A[5]=0}
解：A[2]=0,A[3]=0,A[4]=0,A[5]=0得无解淘汰
解: A[2]=1,A[3]=0,A[4]=0,A[5]=0
A[0]=2.520702
A[1]=1.221316 <K,M>=<3.000000,4.698408>
list={A[1]=1,A[2]=1,A[3]=0,A[4]=0,A[5]=0;A[1]=2,A[2]=1,A[3]=0,A[4]=0,A[5]=0;A[3]=1,A[4]=0,A[5]=0;A[4]=1,A[5]=0}
解:A[1]=1,A[2]=1,A[3]=0,A[4]=0,A[5]=0得
A[0]=3,<K,M>=<3.000000,4.787492>,设B=4.787492
list={A[1]=2,A[2]=1,A[3]=0,A[4]=0,A[5]=0;A[3]=1,A[4]=0,A[5]=0;A[4]=1,A[5]=0}
解A[1]=2,A[2]=1,A[3]=0,A[4]=0,A[5]=0得
A[0]=2,<K,M>=<3.033333,5.192957>,M>B淘汰
list={A[3]=1,A[4]=0,A[5]=0;A[4]=1,A[5]=0}
解A[3]=1,A[4]=0,A[5]=0
A[0]=2.157193
A[1]=0.991967
A[2]=0.359729，<K,M>=<3.000000,5.109912>，M>B淘汰
list={A[4]=1,A[5]=0}
解A[4]=1,A[5]=0，得
A[0]=2.189775
A[1]=1.012524
A[2]=0.373761
A[3]=0.136221，<K,M>=<3.000000,5.894722> M>B淘汰
所以得到最优解对应B=4.787492,这时A[0]=3,A[1]=1,A[2]=1,A[3]=0,A[4]=0,A[5]=0
对应整数为exp(B)=120
可以看出，上面的算法非常有效，几乎看不到一步多余的判断
[解决办法]
不过从上面的计算中，我们得到一个估计t的比较好的方法。
因为我们可以发现t太大时，解出的浮点最优解最后一项会变成负数。
现在问题简单了，只要找到第一个使得出现负数最优解的t就可以了，这个非常容易了。
[解决办法]
记 q 为素数
若a= q1^n1
则所有分解除以a ，结果为（1 + 1/q1 + 1/q1^2 +...)
若a= q1^n1 * q2^n2
则所有分解除以a 结果记为 x=（1 + 1/q1 + 1/q1^2 +...1/q1^n1)（1 + 1/q2 + 1/q2^2 +...1/q2^n2)
记 x1=（1 + 1/q1 + 1/q1^2 +...)
x2=（1 + 1/q2 + 1/q2^2 +...)
若此时增加 n1' = n1+1 则dx= x/x1 * 1/q1^(n1+1) a=q1*a
若此时增加 n2' = n2+1 则dx= x/x2 * 1/q2^(n2+1) a=q2*a
若增加q3
a= q1^n1 * q2^n2 * q3
则dx= x * 1/q3 a=q3*a
取 dx/qn最大的即为增加速度最快的方向
在找到相应的数例如 abundance(a)=4.2 同理，可按这个方法取 qn/dx 最大的 来去除因子 使其abundance(a)=4 即可
[解决办法]
- Assembly code
记p(t)为第t个素数,q(t)=1/p(t)集合N[n;u(1),u(2),...,u(s)] (s<=n)表示那些所有素因子全部为前n个素数，而且后面s个素因子的指数分别为u(1),u(2),...,u(s)的整数构成的集合比如N[2;]表示所有素因子都是2或3的整数集合，比如1,2,3,4,6,8,9,...而N[3;1]表示所有素因子都是2,3或5而且5的指数必须是1的整数，比如5,10,15,20,30,...而对于记号N[n;u(1),u(2),...,u(s)],有时如过u(1),u(2),...,u(s)中有些连续数据，我们可以用*加上次数来缩写比如N[5;1*3]表示N[5;1,1,1], N[5;2*2,1*2] 表示N[5;2,2,1,1]对于整数x,记abd(x)为x的最小abundant数，也就是如果x=p(1)^a(1)*p(2)^a(2)*...*p(t)^a(t),那么abd(x)=(p(1)-p(1)^(-a(1)))/(p(1)-1)*(p(2)-p(2)^(-a(2)))/(p(2)-1)*...*(p(t)-p(t)^(-a(t)))/(p(t)-1)或者abd(x)=(1-q(1)^(a(1)+1))/(1-q(1))*(1-q(2)^(a(2)+1))/(1-q(2))*...*(1-q(t)^(a(t)+1))/(1-q(t))对于给定的实数K,对于集合N[n;u(1),u(2),...,u(s)]记b[K,n;u(1),u(2),...,u(n-s)]为集合N[n;u(1),u(2),...,u(n-s)]中abd(x)>=K的整数x中最小的整数记 K*(1-q(s+1))*(1-q(s+2))*...*(1-q(n))x1=---------------------------- (1-q(s+1)^(u(1)+1))(1-q(s+2)^(u(2)+1))....(1-q(n)^(u(n-s)+1))p(h)为使得 （1-q(h)^2)^h----------------------------*(1+q(h+1))(1+q(h+2))*...*(1+q(s))<=x1的最小的h(1-q(1))(1-q(2))...(1-q(h)) x1*(1-q(1))(1-q(2))...(1-q(h))x2=1-pow(------------------------------------, 1/h) = q(h)^(d+1) （这里实数d>=1) (1+q(h+1))(1+q(h+2))*...*(1+q(s))那么记 p(h)^((d+1)*h)m[K,n;u(1),u(2),...,u(n-s)]=-------------------*p(h+1)p(h+2)...p(s)*p(s+1)^u(1)*p(s+2)^u(2)*...*p(n)^u(n-s) p(1)p(2)....p(h)可以理论证明,b[k,n;u(1),u(2),...,u(n-s)]<=m[k,n;u(1),u(2),...,u(n-s)]也就是对于给定的集合N[n;u(1),u(2),...,u(s)]，我们得到了这个集合中最小abundant K数的一个下界。
[解决办法]
我是这么想的：
首先，某数的约数和与其本身的比例是积性函数，也就是满足：当a和b互质，f(ab)=f(a)f(b)。所以，分解成质数乘积的思路或许是必然的，这点大家前面都提到了。
这样问题就变成了类似于性能和成本优化的问题了，这时，我们对数值取对数，将求积变成求和或许是方便的。
接下来，我想问题可以转化为这个样子：
某游戏的武器系统设计如下：1、武器配备了无限种属性，分别叫做2、3、5、7、11……，对于每个属性，初始时都是0级，对应的攻击力为0；2、每个属性都可以升级。对于某个属性p，每升一级，都要花费lnp的金钱（ln是自然对数）；3、我们有一个级别和攻击力的对应表：对于n级的p属性，共献的攻击力为ln((p^n-1)/(p^n-p^(n-1)))。
直观起见，我们列出一些攻击力的数据如下，为了加强真实性，数据都乘上了10^5，呵呵：
属性2, 属性3, 属性5, 属性7, 属性11,
0级 0, 0, 0, 0, 0
1级 40547, 28768, 18232, 13353, 8701
2级 55962, 36772, 21511, 15123, 9456
3级 62861, 39304, 22154, 15373, 9524
4级 66140, 40134, 22282, 15409, 9530
5级 67740, 40409, 22308, 15414, 9531
6级 68530, 40501, 22313, 15415, 9531
升级费 69315, 109861, 160944, 194591, 239790
我们的问题有两个：1、如果你玩这个游戏，升级武器的路线是什么？2、如果达到给定的攻击力，就可以战胜某个boss通关，如何以最小的花费去通关？
原问题实际相当于第2个问题，而我们偏要先说第1个问题。在武器升级路线上我们首要考虑的性价比。下面的表里列出了攻击力提升和花费的比例，虽说是比例，我们还是乘上了10^5。
属性2, 属性3, 属性5, 属性7, 属性11,
1级 58496, 26186, 11328, 6862, 3629
2级 22239, 7286, 2037, 910, 315
3级 9954, 2305, 400, 129, 28
4级 4731, 755, 80, 18, 3
5级 2308, 250, 16, 3, 0
6级 1140, 83, 3, 0, 0
因此，我们按照性价比大的方向来升级呗。
按照这个思路可以得到下面的结果：
k 下面的pn表示前n个质数的乘积，例如：p1=2，p2=6，p3=30……
2 p2
3 p1^2 p3
4 p1^2 p2 p5
5 p1^2 p2 p3 p7
6 p1^3 p2 p3 p11
7 p1^3 p2^2 p4 p16
8 p1^3 p2^2 p3 p5 p25
9 p1^4 p2^2 p4 p6 p38
10 p1^4 p2^2 p3 p4 p8 p59
11 p1^5 p2^2 p3 p4 p10 p93
12 p1^5 p2^3 p4 p5 p12 p148
13 p1^4 p2^3 p3 p4 p6 p15 p237
14 p1^5 p2^3 p3 p4 p7 p19 p382
15 p1^5 p2^3 p3 p4^2 p8 p24 p619
16 p1^6 p2^3 p3 p4 p5 p9 p30 p1009
17 p1^6 p2^4 p3 p4 p6 p11 p37 p1653
18 p1^7 p2^3 p3^2 p4 p6 p12 p47 p2719
19 p1^7 p2^4 p3 p4 p5 p8 p15 p61 p4485
20 p1^8 p2^3 p3^2 p4 p5 p8 p16 p74 p7431
21 p1^8 p2^4 p3^2 p4 p6 p9 p20 p95 p12341
上面的结果只算到了21，我想优化后的程序的运行时间会是小于1分钟的，在感觉上，如果思路正确，或许突破23是可能的。
对于第2个问题，实际上可以沿用上面的思路，只要搜索的宽一些就可以了吧。
[解决办法]
先提供一个比较粗糙的版本
不过里面的贪心算法不够好，所以只有对K>=6才能够计算
还有搜索的代码没有对搜索的各个项目先进行排序，这些都可以优化
- C/C++ code
#include <stdio.h>#include <math.h>#define RUN_LENGTH 2000#define MAX (RUN_LENGTH*RUN_LENGTH)int not_prime[MAX];#define the_prime not_primeint count;int best[MAX];int best_n;int b[MAX];double KK;//GREEDY_ITEM gItems[MAX];int greedyInit(int n, double k, double *mr,int* smin,int* smax){ int i,t; double abd=k; double left; double u; double min_m; int cur_index=-1; int estate; for(i=0;i<n;i++){ abd*=(1-1.0/the_prime[i]); if(abd<1.0)break; } *smin=i; //printf("Start from the %dth prime for k=%f\n",i+1, k); for(t=i;t<n-1;t++){ double w=abd; double m=0.0; int isend=0; for(i=t;i>0;i--){ double x=1.0-1.0/((double)the_prime[i]*the_prime[i]); double logi=log(the_prime[i]); m+=logi; w/=x; if(w>=1.0){ isend=1; estate=-1; break; } u=pow(x,i); if(u<=w){ double d; m-=logi; d=1.0-pow(w,1.0/i); d=1.0/d; m+=i*log(d); break; } } if(isend)break; if(2*the_prime[i]*the_prime[i]<the_prime[t]){ estate=0; isend=1; } for(;i>=0;i--){ m-=log(the_prime[i]); } if(cur_index<0||m<min_m){ min_m=m; cur_index=t; } abd*=(1-1.0/the_prime[t+1]); if(isend)break; } *smax=t; *mr=min_m; return cur_index;}int greedyBound(int n, int s, int *b,double k, double *m, int next_index){ int i,t; double abd=k; double left; double u; double min_m; double this_m=0.0; int cur_index=-1; int estate; for(i=n-s;i<n;i++){ abd/=(1.0-1.0/pow(the_prime[i],b[i]+1)); this_m+=b[i]*log(the_prime[i]); } for(i=0;i<n;i++){ abd*=1.0-1.0/the_prime[i]; } if(s==n-1){///Special case, calculate order of 2 and return if(abd>=1.0){ b[0]=0; return -1; }else{ double v=ceil(-log(1.0-abd)/log(the_prime[0])); b[0]=(int)(v-1); this_m+=b[0]*log(the_prime[0]); } *m=this_m; return 0; } //printf("Start from the %dth prime for k=%f\n",i+1, k); for(t=n-s;t>=1;t--){ double w=abd; double m=this_m; int isend=0; for(i=t;i<n-s;i++){ w/=(1.0-1.0/pow(the_prime[i],next_index+1)); m+=next_index*log(the_prime[i]); } if(w>=1.0)continue; u=pow(w, 1.0/t); u=1.0-u; if(u>=pow(1.0/the_prime[t-1],next_index+1)) continue; m+=-log(u)*t; for(i=0;i<t;i++){ m-=log(the_prime[i]); } if(cur_index==-1||min_m>m){ min_m=m; cur_index=t; } } if(cur_index>=0){ b[0]=0; for(i=cur_index;i<n-s;i++){ b[i]=next_index; } *m=min_m; return cur_index; }else{ return -1; }}void find(int n,int s,int step, double *bound){ double bound2; int r,i; r=greedyBound(n,n-s, b, KK, &bound2,step+1); if(r==0){ if(bound2<*bound){ *bound=bound2;///Find a better result best_n=n; memcpy(best,b,sizeof(b[0])*n); } } if(r<0||bound2>=*bound) return; for(i=r;i>=0;i--){ b[i]=step+1; find(n,i,step+1, bound); } for(i=r+1;i<s;i++){ find(n,i,step+1, bound); }}int main(int argc, char *argv[]){ int i; int start,end; int min=1; int step; double v; double bound; int r; sscanf(argv[1],"%lf",&KK); not_prime[0]=not_prime[1]=1; for(i=4;i<MAX;i+=2)not_prime[i]=1; for(i=3;i<RUN_LENGTH;i+=2){ if(!not_prime[i]){ int j; for(j=i*i;j<MAX;j+=2*i){ not_prime[j]=1; } } } for(i=0;i<MAX;i++){ if(!not_prime[i]){ the_prime[count++]=i; } } i = greedyInit(count,KK, &bound,&start,&end); step = 0; min = i-1; if(i<0){ fprintf(stderr,"Out of calculation boundary\n"); return -1; } r=i; while(r>0){ step++; r=greedyBound(i,i-r, b, KK, &bound,step); if(r>=0){ min=r-1; }else{ fprintf(stderr,"Searching failed\n"); return -1; } } best_n=i; memcpy(best,b,sizeof(b[0])*i); for(i=end;i>=start;i--){ find(i,i,0,&bound); } printf("Best boundary %f\n",bound); start=-1; for(i=0;i<best_n;i++){ if(start==-1){ start=i; continue; } if(best[i]!=best[start]){ if(best[start]-best[i]==1){ printf(" p%d",start+1); }else{ printf(" p%d^%d",start+1, best[start]-best[i]); } start=i; } } printf(" p%d\n",best_n); for(i=0;i<best_n;i++){ if(i==0){ start=0; continue; } if(best[i]!=best[start]){ if(start!=i-1) printf(" P[%d-%d]^%d",start+1,i,best[start]); else printf(" P[%d]^%d",start+1,best[start]); start=i; } } printf(" P[%d-%d]\n",start+1,best_n); v=1.0; for(i=0;i<best_n;i++){ v*=1.0-1.0/pow(the_prime[i],best[i]+1); v/=1.0-1.0/the_prime[i]; } printf("%f\n",v);}