最多约数问题答案解析与勘误

问题描述： 
正整数x的约数是能整除x的正整数。正整数x 的约数个数记为div(x)。例如，1，2，5，10 都是正整数10 的约数，且div(10)=4。设a 和b 是2 个正整数，a≤b，找出a和b之间约数个数最多的数x。 
编程任务： 
对于给定的2个正整数a≤b，编程计算a 和 b 之间约数个数最多的数。 
数据输入： 
输入数据由文件名为input.txt的文本文件提供。文件的第1 行有2 个正整数 a和 b。 
结果输出: 
程序运行结束时，找到a 和b之间约数个数最多的那个数及最多约数个数。

测试数据：【只给出最多约数个数, time limit: 1s】

问题分析

下面的解析摘自网上一篇博客（来源http://blog.csdn.net/sarah_jieyu/article/details/6512030），标红部分是我自己添加上去的。

本题的要求是，求出一个给定区间内的含约数最多的整数。

    首先明确一下如何求一个数的约数个数：若一个数N满足：N = A1^N1 * A2^N2 * A3^N3 * …… * Am^Nm，则n的约数个数为(N1 + 1) (N2 + 1) (N3 + 1) …… (Nm + 1)。这是可以用乘法原理证明的。

    最浅显的算法是，枚举区间内的每个整数，统计它们的约数个数。这个算法很容易实现，但是时间复杂度却相当高。因为题目约定区间的最大宽度可以达到10^9的数量级，相当庞大。因此，在极限规模时，时间是无法忍受的。所以，我们需要尽量的优化时间。

    分析一下枚举的过程就会发现，如果我们枚举到两个数n和m*n（m为一相对于n较大的质数），那么我们将重复计算n的约数两次。据此，我们发现了枚举效率低的根本所在。为了解决这一重复，我们可以选取另一种搜索方法——以质因子为对象进行深度搜索。

    初始时，令number := 1，然后从最小的质数2开始枚举，枚举包含一个2、两个2……n个2的情况……直至number * 2ⁿ大于区间的上限（max）。对于每种“2^k的情况”，令number := number * 2ⁿ，再枚举3的情况，然后，枚举5的情况、7的情况……方法相同。整个过程是一个深度搜索的过程。当number大于等于区间下限（min）时，我们就找到了一个区间内的数，根据前面介绍的方法，可以得到它的约数个数。所有的区间内的数的约数个数的最大值就是我们要求的目标。

    为什么这种深度搜索可以减少常规枚举过程中的重复问题呢？请看下面的一个例子

     设给定的区间为[6,30],6,18,30为区间内的数，按照常规枚举方法，计算18,30，的时候分别计算了因子6的约数个数，重复计算2次。如果使用上述所说的深度搜索方法，求这3个数的因数个数的路径有一条公共部分，2*3，这一部分只计算了一次，求18只需再乘个3，求30只需再乘个5，相对于常规枚举减少了两次计算2*3的时间。但这种深度搜索也有问题，就是number有可能是无用的，下面的分析便是对这种深搜方法进行无用数据剪枝。

    值得注意的是，我们枚举过程中得到的number可能无用的，即无论用number去乘以多少，都无法得到区间内的数。这样的number如果继续枚举下去，无疑会大大降低效率。那么，能否通过简单的判断，将其剪去呢？答案是可以的。很容易证明，如果(min – 1) div number < max div number，则区间内存在可以被number整除的数。因为，如果区间[min, max]内存在可以被number整除的数，也即是从min到max中至少有一个数能被number整除，那么区间[min – 1, max]内的数被number除得的商肯定不止一种，所以(min – 1) div number必然小于max div number。

反过来，如果(min-1)div number=max div number,则[min,max]内不存在可以被number整除的数。

证明如下：

假设[min,max]内存在可以被number整除的数

如果(min-1) div number=max div number,则min div number = max div number，那么①min等于max，且min和max均可以被number整除，或者②max>min,min可以被number整除，①的情况下可推出，(min-1) div nunber<max div number,与条件矛盾，舍去；②的情况下也可以推出(min-1) div number<max div number,舍去。所以结论成立。因此，我们只需枚举符合要求的number；至于不符合的，可以剪去。

    此外，我们枚举的质数可能会达到很大。因为给出的整数最大可以达到4,000,000,000，它的质因数自然最大也可以到100,000,000的数量级。如果按上面的方法枚举，显然无法承受时间的压力。但是，我们又可以看到，对于区间内的任一整数，它包含的大于SQRT(4,000,000,000) = 2*31623的质因数最多只可能有一个。因此，我们只需枚举小于2*31623的质数。如果对一个符合要求的number（即，可以证明区间内至少存在一个数可以被number整除），无法找到一个小于2*31623的质数p使得number * p * x ∈ [min, max]（x为一正整数）。那么，可知number需要乘以一个大于31623的质数才能得到number’，使得number’ ∈ [min, max]。根据前面介绍的乘法原理，只需将number包含的约数个数乘以2，即得number’包含的约数个数。

    我们还能看到，如果当前搜索状态为(from, number, total)，其中from是指当前枚举到的质因子（按从小到大枚举），total是指number中包含的约数个数。那么剩下的因子数最多为q = [log_from(max / number)]，这些因子组成的约数个数（即上述求约数个数时用到的一串乘积）最大为2^q。当前所能取到的（理想情况）最大约数个数就是total * 2^q，如果这个数仍然无法超过当前最优解，则这一分支可以剪去。

 

    深度搜索的过程，从表面上看是一个指数级的复杂度。其实不然，更准确的说，复杂度应为O(p^logn)。因为，它的指数增长速度是随n成对数级增长，本质上说，还是多项式级的算法。而且我们还进行了一些剪枝，由于枚举中的分支多数为非法的，因此通过剪枝可以去掉绝大多数的分支。这样就大大提高了程序的效率。尽管数据规模很大，但对于极限数据仍然可以做到“一闪即出”。

根据以上分析给出以下代码

#include<iostream>
using namespace std;


void primes(); //用筛选法产生质数存于prim数组中
void search(long from, long total, long num, long low, long up);


const long MAXP = 100000;
long prim[MAXP];
//long max, numb ,PCOUNT; //max存放最多约数个数，numb存放约数个数最多的数
long max_num;
long numb,PCOUNT;


int main()
{
primes();
long low, up;
cin >> low >> up;
if ((low == 1) && (up == 1))
{
max_num = 1;
numb = 1;
}
else
{
max_num = 2;
numb = low;
search(1, 1, 1, low, up);
}
cout << max_num << endl << numb << endl;
return 0;
}


//求小于100000的素数表
void primes() 
{
bool get[MAXP+1];
long i;
for (i = 2; i <= MAXP; i++)
get[i] = true;


for (i = 2; i <= MAXP; i++)
if (get[i])
{
long j = i + i;
while (j <= MAXP) 
{
get[j] = false;
j += i;
}
}
long ii, j;
for (ii = 2, j = 0; ii <= MAXP; ii++)
if (get[ii]) prim[++j] = ii;
PCOUNT = j;
}


// 区间[low,up]上，tot为当前约数最多个数，num为约数个数最多的数，
//from表示现在是第几个质数。
void search(long from, long total, long num, long low, long up)
{
/*cout<<"from is "<<from<<endl;
cout<<"total is "<<total<<endl;
cout<<"num is "<<num<<endl;
cout<<"low is "<<low<<endl;
cout<<"up is "<<up<<endl;
cout<<"numb is "<<numb<<endl<<endl;*/
if (num >= 1){
//下面加注释的是网上广泛流传的代码，其实是错的，因为这种条件判断不能保证最后得到的因数个数最多的数大于区间的下边界
/*
if ( (tot > max) || ((tot == max) && (num < numb)) )
{
max = tot;
numb = num;
}*/
//这是我修改后的代码
if ( (total >=max_num&&num>=numb) )
{
//cout<<"total"<<endl;
max_num = total;
numb = num;
}
}
/**
此处代码较为简练，但包含的内容很多，较难理解.
因为所给区间[low,up]内有可能包含一个数a，a的某个因子大于100000，超出了素数表的搜索范围
有可能搜索不到a,怎么解决这个问题？可以看出来，如果a含有大于100000的因子，则a的因数个数为当前搜索到的因数个数乘以2
如果，一个区间内包含了一个和a含有除了大于100000的因子的数以外的其他因子的数b，那么b的因数个数一定多于a.因此对于
含有a的区间，我们只考虑[a,a]这种可能，其他的区间情况，a不可能是含有因数最多的数。
此处给出的判断条件是low==up,low>num。就包含了[a,a]这种情况，当然此处会有多余的情况，如low=up=b,b可以分解为若干个较小的因子的乘积
但并不影响最终的结果。
*/
if ((low == up) && (low > num)) 
search(from, total*2, num*low, 1, 1);




/*此处就是深度搜索的代码 ，有一个地方需要补充说明一下，在一条搜索路径上，区间边界不断除以一个素数，
会越来越小，因为每次使用的都是除法取整，会有误差。但这里可以证明，累积误差不会导致的除法结果的整数部分不同。
即：a/(x1*x2*..xn)=a/x1*1/x2*1/x3*..*1/xn,证明：设a=c(x1*x2*..*xn)+r(r为余数)，则a/x1*1/x2*..1/xn=c+r/x1*1/x2*..*1/xn,整数部分仍然为c
*/
for (long i = from; i <=PCOUNT; i++)
{
if (prim[i] > up) return;  
else
{
long j = prim[i], x = low - 1, p=low,y = up, n = num, t = total, m = 1;
while (true)
{
m++;
t += total;
x /= j;
y /= j;
p/=j;
/*
此处x==y的判断剪枝根据是如果(min-1)/number=max/number,则[min,max]内不包含可以被number整除的数
*/
if (x == y) 
break;
n *= j;
/*
勘误，此处网上流传的代码为search(i+1, t, n, x+1, y);,这个是不对的，因为(low-1)/j+1和low/j整数部分会产生误差，
例，low=9,j=2,(low-1)/j+1=5,而low/j=4,所以此处进行了更正。下面一些测试数据可以测出错误代码的结果是不正确的
[9999999,10000000],[99999999,100000000],[999999999,1000000000]
*/
//search(i+1, t, n, x+1, y);
search(i+1, t, n, p, y);
}
/*
此处也是剪枝处理
*/

m = 1 << m;
if (total < max_num / m) return;
}//end else
}//end for
}