求smooth数

1.什么是Smooth 数

如果一个整数的所有素因子都不大于B，我们称这个整数为B-Smooth数。用符号Pn表示表示第n个素数，例P1=2,P2=3,P3=5. 我们称素因子不大于Pn的Smooth数叫Pn-Smooth数。关于Smooth数更多的信息，请请参阅 Wiki词条 Smooth number http://en.wikipedia.org/wiki/Smooth_number.  P1 到 P9-Smooth 数可参阅下面的数列，这些数列在On-Line Encycolpedia of Integer Sequences http://oeis.org 给出

    2-smooth numbers: A000079 (2^i, i>=0)

    3-smooth numbers: A003586 (2^i ×3^j,i,j>=0)

    5-smooth numbers: A051037 (2^i ×3^j×5^k, i,j,k>=0)

    7-smooth numbers: A002473 (2^i×3^j×5^k×7^l, i,j,k,l>=0)

    11-smooth numbers: A051038 (etc...)

    13-smooth numbers: A080197

    17-smooth numbers: A080681

    19-smooth numbers: A080682

    23-smooth numbers: A080683

 

2. 如何得到smooth数

  基本上，我们有2种方法来生成一定范围内的smooth数。

1.   穷举法。基本方法是对指定范围内的每个数做分解素因数，如果其所有素因子都小于等于Pn,那个这个数就是Pn-smooth数。这种方法的优点是简单，缺点是速度较慢，只适合检查小范围内的整数。

 

2．构造法。其基本方法是用一个数组表示各个素因子的指数。例arr是一个数组，arr[1]表示素数p1（2）的指数，arr[2]表示素数p2（3）的指数，arr[3]表示素数p3（5）的指数. 欲求min和max之间的所有pn-smooth数, 则需要列举枚举这个数组的各个元素的取值。数组的各个元素的值需满足如下2个条件。

1. p[i]>=0, 且p[i] ^ arr[i] <= max.  （p[i]表示第i个素数，i>=1,i<=n）

2. min <=p[1]^arr[1] * p[2] ^ arr[2] * p[3] ^ arr[3] * …p[n]*arr[n] <=  max

对于一个包含n个元素的数组arr, 每个元素arr[i] （1<=i<=n）各取一个数，就得到一个n维向量。如arr[1]=0, arr[2]=0, arr[3]=0可用3维向量(0,0,0)来表示，而arr[1]=1, arr[2]=0, arr[3]=0 则用3维向量(1,0,0)来表示。为了使得取得的所有n维向量既不重复也不遗漏。我们这里定义了比较2个n维向量的大小的规则

1.  如果2个n维向量a,b的所有分量都是相同的，我们说这两个n维向量是相同的。

2.  对于2个n维向量a,b（a的各个分量为a[1],a[2], …, a[n]）,如果a[n]>b[n]我们说a>b.

3.  如果存在某个i，使得（1）对于所有的j (i<j<=n)都有a[j]=b[j]. (2) a[i]>b[i],那么我们说n维向量a>b.

  我们列举每个n维向量时，总是从小到到的进行，先是(0,0,0), 然后是(1,0,0), 然后是(2,0,0)…，然后是(0,1,0), (1,1,0), (2,1,0).  

未完待续

 

代码：

#include <stdio.h>
#include <limits.h>
#include <string.h>
#include <malloc.h>

#define MAX_LEVEL		32

typedef unsigned long DWORD;

#if defined(_MSC_VER)
	typedef unsigned __int64 UINT64;
	#define UINT64_FMT  "%I64u"

#elif defined(__GNUC__) 
	typedef unsigned long long UINT64;
	#define UINT64_FMT  "%llu"
#else
	#error don't support this compiler
#endif

const DWORD primes[]=
{
	 2,  3,  5,  7, 11, 13, 17, 19,
	23, 29, 31, 37, 41, 43, 47, 53,
	59, 61, 67, 71, 73, 79, 83, 89,
	97,101,103,107,109,113,127,131
};

typedef  struct _searchContent_tag
{
	struct search_tag
	{
		UINT64 PI;		// product of some numbers
		UINT64 max;
		int idx;		// the index in array pArr
		int maxIdx;		// the index of last element
		UINT64* pArr;	// the address of primes power table
	}arr[MAX_LEVEL];
	UINT64 min;
	UINT64 max;
	int  maxLevel;		//level
}SRCH_CONTENT_ST;

void initSearchContent(int level, UINT64 min, UINT64 max, SRCH_CONTENT_ST* p)
{
	UINT64 prime,limit,pp;
	int i,c;
	
	memset(p,0,sizeof(SRCH_CONTENT_ST));
	
	p->maxLevel=level-1;
	p->max=max; 	p->min=min;
	p->arr[level-1].PI=1;
	p->arr[level-1].max=max;
	for (i=0;i<level;i++)
	{
		prime=primes[i];
		limit=ULLONG_MAX/prime;
		
		c=1;pp=1;
		while (pp<=limit)	//pp is power of prime
		{ pp*=prime;c++; }
		
		p->arr[i].pArr=(UINT64 *)malloc( sizeof(UINT64) * (c+1));
		
		p->arr[i].pArr[0]=pp=1;	c=1;
		while (pp<=limit)
		{ pp*=prime;p->arr[i].pArr[c++]=pp; }
		
		p->arr[i].maxIdx= c-1;
	}
}

void freeSearchContent( SRCH_CONTENT_ST* p)
{
	int i;
	for (i=0;i<=p->maxLevel;i++)
	{
		if ( p->arr[i].pArr!=NULL)
			free(p->arr[i].pArr);
	}
	memset(p,0,sizeof(SRCH_CONTENT_ST));
}

//recursive form
void searchSmoothNums1( SRCH_CONTENT_ST* p, int order)
{
	int    idx;
	/*
	    recurrence formula
		p->arr[order].PI= 1;
		p->arr[i-1].PI    = arr[i].pArr[ arr[i].idx] * arr[i+1].PI;	
		p->arr[i].max= p->max / p->arr[i].PI;
	*/

	if (order==0)
	{
		UINT64 min;
		UINT64 prod;

		min= (p->min + p->arr[0].PI -1)/p->arr[0].PI;  // min >= p->min/PI( when p->min % PI==0, min= p->min/PI)
		
		if ( min > p->arr[0].max)
			return;

		idx=0;
		while (idx < p->arr[0].maxIdx && p->arr[0].pArr[idx]<min )
			idx++; 
		
		while ( p->arr[0].pArr[idx] <= p->arr[0].max && idx <= p->arr[0].maxIdx )
		{
			prod= p->arr[0].PI * p->arr[0].pArr[idx];
			idx++;
			printf(UINT64_FMT"\n",prod);
		}
	}
	else 
	{
		if ( p->arr[order].max==1 )
		{
			printf(UINT64_FMT"\n",p->arr[order].PI);
			return ;
		}
		
		for (idx=0;  
			idx<= p->arr[order].maxIdx  && p->arr[order].pArr[idx] <= p->arr[order].max;
			idx++ )
		{
			p->arr[order].idx=idx;  
			p->arr[order-1].PI = p->arr[order].PI * p->arr[order].pArr[idx];
			p->arr[order-1].max= p->max / p->arr[order-1].PI;
			searchSmoothNums1(p,order-1);   //call itselft
		}
		
	}
}


#define M_FILL   0    //fill mode
#define M_INC    1    //increase mode

// iterative form
// search all p[order+1] smooth number which between p->min and p->max
// For exmaple
//   order=0, search 2-smooth numbers
//   order=1, search 3-smooth numbers
//   order=2, search 5-smooth number
void searchSmoothNums2( SRCH_CONTENT_ST* p, int order)
{
	UINT64 prod;
	UINT64 PI;
	UINT64 min;
	int idx;

	int i=order;
	int mode=M_FILL;
	/*
	    recurrence formula
		p->arr[order].PI= 1;
		p->arr[i-1].PI   = arr[i].pArr[ arr[i].idx] * arr[i+1].PI;	
		p->arr[i].max    = p->max / arr[i].PI;	
	*/
	
	p->arr[order].PI=1;
	p->arr[order].max= p->max;

	while (i<=order)
	{
		
		if (mode==M_FILL)  //fill mode
		{
			PI = p->arr[i].PI;
			if (i>0)
			{
				p->arr[i].idx=0;
				p->arr[i-1].PI = p->arr[i].PI;	//p->arr[i-1].PI= p->arr[i].PI;
				p->arr[i-1].max= p->arr[i].max; 
				i--;						//forward to next level
			}
			else  // i==0
			{
				min= (p->min + PI-1)/PI;	// min >= p->min/PI( when p->min % PI==0, min= p->min/PI)
			
				if ( min > p->arr[0].max)
				{
					i++;					// backtrace
					mode=M_INC;
				}
				else
				{
					idx=0;
					while (idx < p->arr[0].maxIdx && p->arr[0].pArr[idx]<min )
						idx++; 
					if (  p->arr[0].pArr[idx] > p->arr[0].max || idx > p->arr[0].maxIdx)  
						i++;					// backtrace
					else   //in this case, varible i don't change
					{
						p->arr[0].idx=idx;
						printf(UINT64_FMT"\n",PI * p->arr[0].pArr[idx]);

					}
					mode=M_INC;		//change to increase mode
				}
			}
		}
		else  //is increase mode
		{
			idx= p->arr[i].idx + 1;
			if (p->arr[i].pArr[idx] > p->arr[i].max || idx> p->arr[i].maxIdx)
				i++;		// backtrace, mode don't change, still is increase mode
			else
			{
				p->arr[i].idx=idx;
				prod= p->arr[i].PI * p->arr[i].pArr[idx];

				if (i==0)
					printf(UINT64_FMT"\n",prod);
				else
				{
					p->arr[i-1].max= p->max / prod;
					if ( p->arr[i-1].max ==1)
					{
						printf(UINT64_FMT"\n",prod);
						i++;			// backtrace, mode don't change, still is increase mode
					}
					else
					{
						p->arr[i-1].PI=prod;
						mode=M_FILL;	// change mode to fill mode
						i--;			// forward to next level
					}
				}
			}
		}
	}
}

void ref_searchSmoothNums(int order,UINT64 min, UINT64 max)
{
	UINT64 i,m;
	int j,isSmooth;

	for (i=min;i<=max;i++)
	{
		for (m=i,isSmooth=0,j=0;j<order;j++)
		{
			while (m % primes[j]==0) m/=primes[j];
			if (m==1)
			{	isSmooth=1;	break;	}
		}
		if (isSmooth)
			printf(UINT64_FMT"\n",i);
	}
}

void test()
{
	SRCH_CONTENT_ST searchContent;
	UINT64 nMin,nMax;
	int cs;//case

	for (cs=0;cs<=1;cs++)
	{
		if (cs==0)
		{ nMin=1;  nMax=100; }
		else if (cs==1)
		{ nMin=100;	nMax=200;	}

		initSearchContent(4,nMin,nMax,&searchContent);			
		
		printf("\n\nSearch smooth number ref_searchSmoothNums-------\n"); 
		ref_searchSmoothNums(4,nMin,nMax);
		
		printf("\n\nSearch smooth number searchSmoothNums1-------\n"); 
		searchSmoothNums1( &searchContent,4-1);
		
		printf("\n\nSearch smooth number searchSmoothNums2-------\n"); 
		searchSmoothNums2( &searchContent,4-1);

		freeSearchContent( &searchContent);
	}
}

int main(int argc, char* argv[])
{
	test();
	return 0;
}