寻找快速开方的源码

[貌似掉线]

想找一个快速开方的算法，是对整数开方的，而且只要求精确到整数位即可，不知谁能提供一下，用于统计素数个数的。

[djlxue]

有现成的干嘛不用，最简单的sqrt就可以了

[天下第一好大人]

单纯算素数，用筛数法吧。

[jackyjkchen]

[Quote=引用 2 楼 gogdizzy 的回复:]
单纯算素数，用筛数法吧。
[/Quote]
算法

void InitTable(unsigned int n)
{
    Table = (unsigned char*)malloc( n+14 );
	if(Table == NULL)
	{
		fprintf(stdout,"内存不足！\n");
		exit(-1);
	}
}

void FreeTable()
{
    free(Table);
}

void TestPrime(unsigned int n)
{
	unsigned int i = 0, j = 0;
	unsigned char *p = Table+1;
	unsigned int mpLen = 2*3*5*7*11*13;
	unsigned char *magicPattern = (unsigned char*)malloc(mpLen);
	unsigned int remainder = n%mpLen;
	unsigned char *pstop = p+n-remainder;
	unsigned int step = 0, stop = n/13;
	for (i=0; i<mpLen; i++)
	{
		magicPattern[i++] = 1;
		magicPattern[i++] = 0;
		magicPattern[i++] = 0;
		magicPattern[i++] = 0;
		magicPattern[i++] = 1;
		magicPattern[i] = 0;
	}
	for (i=4; i<=mpLen; i+=5)
		magicPattern[i] = 0;
	for (i=6; i<=mpLen; i+=7)
		magicPattern[i] = 0;
	for (i=10; i<=mpLen; i+=11)
		magicPattern[i] = 0;
	for (i=12; i<=mpLen; i+=13)
		magicPattern[i] = 0;

	while (p < pstop)
	{
		memcpy(p, magicPattern, mpLen);
		p += mpLen;
	}
	if (remainder > 0)
		memcpy(p, magicPattern, remainder);
	Table[2] = 1;
	Table[3] = 1;
	Table[5] = 1;
	Table[7] = 1;
	Table[11] = 1;
	Table[13] = 1;

	for (i=17; i <= stop; i++)
	{
		if (0 == Table[i])
            continue;
		step = i*2;
		for (j=i*17; j <= n; j+=step)
			Table[j] = 0;
	}
	free(magicPattern);
}

调用

int main()
{
	LARGE_INTEGER begintime,endtime,freqtime;
	__int64 resulttime = 0;
	unsigned int top = 0, n = 0, i = 0;
	char filename[MAX_PATH],tmpbuf[MAX_PATH];
	FILE *outfile = NULL;
	do
	{
		fprintf(stdout,"请输入上限：    ");
		if(scanf("%u",&top) <= 0)
		{
			fprintf(stdout,"请输入数字！\n");
			fgets(tmpbuf,MAX_PATH,stdin);
			continue;
		}
		else if(top < 2)
			fprintf(stdout,"请输入>=2的值！\n");
	}while(top < 2);
	fprintf(stdout,"请输入保存位置：");
	fgets(tmpbuf,MAX_PATH,stdin);
	filename[0] = '\0';
	fgets(filename,MAX_PATH,stdin);
	filename[strlen(filename)-1] = '\0';
	outfile = fopen(filename,"w");
	QueryPerformanceCounter(&begintime);

	InitTable(top);
	TestPrime(top);
	for (i=2; i<=top; i++)
	{
	    if(Table[i])
	    {
	        n++;
            if(outfile)
                fprintf(outfile,"%u,",i);
	    }
	}
	FreeTable();

	QueryPerformanceCounter(&endtime);
	QueryPerformanceFrequency(&freqtime);
	fprintf(stdout,"在<=%u时共有%u个素数\n",top,n);
	if(outfile)
	{
		fprintf(outfile,"\n");
		fprintf(outfile,"在<=%u时共有%u个素数\n",top,n);
		fclose(outfile);
		fprintf(stdout,"结果已保存在：  %s\n",filename);
	}
	resulttime = (endtime.QuadPart - begintime.QuadPart) * 1000 / freqtime.QuadPart;
	fprintf(stdout,"处理时间：      %I64d小时%I64d分%I64d.%03I64d秒\n",resulttime / 3600000,(resulttime / 60000) % 60,(resulttime /1000) % 60,resulttime % 1000);
	return 0;
}

[jackyjkchen]

上面的筛选算法比单纯的开放试模快几百倍倍，当然内存开销很大，当然你也可以试试查表试模，快10倍，内存开销小得多

unsigned int *PrimeTable;
unsigned int TableIndex;

void InitTable()
{
	TableIndex = 0;
	PrimeTable = (unsigned int*)malloc(6543*sizeof(unsigned int));
	PrimeTable[0] = 3;
}

void FreeTable()
{
	free(PrimeTable);
}

bool TestPrime(unsigned int n)
{
	unsigned int i = 0;
	if(n == 1)
		return false;
	else if(n == 2)
		return true;
	else if(!(n % 2))
		return false;
	for(i=0; PrimeTable[i] * PrimeTable[i]<=n; i++)
	{
		if(n % PrimeTable[i] == 0)
			return false;
	}
	if(n<=65521)
		PrimeTable[TableIndex++] = n;
	return true;
}


int main()
{
	LARGE_INTEGER begintime,endtime,freqtime;
	__int64 resulttime = 0;
	unsigned int num = 2, top = 0, n = 1;
	char filename[MAX_PATH],tmpbuf[MAX_PATH];
	FILE *outfile = NULL;
	do
	{
		fprintf(stdout,"请输入上限：    ");
		if(scanf("%u",&top) <= 0)
		{
			fprintf(stdout,"请输入数字！\n");
			fgets(tmpbuf,MAX_PATH,stdin);
			continue;
		}
		else if(top < 2)
			fprintf(stdout,"请输入>=2的值！\n");
	}while(top < 2);
	fprintf(stdout,"请输入保存位置：");
	fgets(tmpbuf,MAX_PATH,stdin);
	filename[0] = '\0';
	fgets(filename,MAX_PATH,stdin);
	filename[strlen(filename)-1] = '\0';
	outfile = fopen(filename,"w");
	QueryPerformanceCounter(&begintime);
	InitTable();
	if(outfile)
		fprintf(outfile,"%u,",num);
	for(num=3; num<=top; num+=2)
	{
		if(TestPrime(num))
		{
			if(outfile)
				fprintf(outfile,"%u,",num);
			n++;
		}
	}
	FreeTable();
	QueryPerformanceCounter(&endtime);
	QueryPerformanceFrequency(&freqtime);
	fprintf(stdout,"在<=%u时共有%u个素数\n",top,n);
	if(outfile)
	{
		fprintf(outfile,"\n");
		fprintf(outfile,"在<=%u时共有%u个素数\n",top,n);
		fclose(outfile);
		fprintf(stdout,"结果已保存在：  %s\n",filename);
	}
	resulttime = (endtime.QuadPart - begintime.QuadPart) * 1000 / freqtime.QuadPart;
	fprintf(stdout,"处理时间：      %I64d小时%I64d分%I64d.%03I64d秒\n",resulttime / 3600000,(resulttime / 60000) % 60,(resulttime /1000) % 60,resulttime % 1000);
	return 0;
}

[jackyjkchen]

即便是开方试模，相信如果代码细节处理不好速度也是会相差几倍的

给个参考实现，和你自己的比比

bool TestPrime(unsigned int n)
{
	unsigned int i = 0, sq = 0;
	if(n == 1)
		return false;
	else if(n == 2)
		return true;
	else if(!(n % 2))
		return false;
	sq = (unsigned int)sqrt((double)n);
	for(i=3; i<=sq; i+=2)
	{
		if(n % i == 0)
			return false;
	}
	return true;
}

[fallening]

unsigned short sqrt(unsigned long a){ 
  unsigned long rem = 0; 
  unsigned long root = 0; 
  unsigned long divisor = 0; 
  for(int i=0; i<16; i++){ 
    root <<= 1; 
    rem = ((rem << 2) + (a >> 30)); 
    a <<= 2; 
    divisor = (root<<1) + 1; 
    if(divisor <= rem){ 
      rem -= divisor; 
      root++; 
    } 
  } 
  return (unsigned short)(root); 
}   

[貌似掉线]

谢谢楼上各位，不过我只想看有没有算开方的函数，等我学习了更多的知识后我会用筛选法，现在只是想用我目前的知识看用穷举法能算到什么程度。sqrt听说是一个慢的函数，所以想看下有没有更快的。或者给个思路也行。

[liyaobinRyan]

帮顶

[BillLeecn]

提醒LZ最好不要怀疑STL里的东西
STL是高高高高高手弄出来的，经过严格测试。
STL里有的不要自己弄（如果你自己能弄出比STL里的更好更快的，那你也就没必要学C++了）

[panghuhu250]

引用 7 楼 maosidiaoxian 的回复:

谢谢楼上各位，不过我只想看有没有算开方的函数，等我学习了更多的知识后我会用筛选法，现在只是想用我目前的知识看用穷举法能算到什么程度。sqrt听说是一个慢的函数，所以想看下有没有更快的。或者给个思路也行。

1. 在你的程序里,性能瓶颈应该不是sqrt函数,所以优化sqrt函数不会有太大作用.应该先进行性能测试,确认sqrt函数确实是性能瓶颈,再考虑怎么优化它.

2. 如9楼所说,STL已经很成熟了,要比STL还快可不大容易.

[taodm]

“sqrt听说是一个慢的函数”，这都是怎么道听途说来的啊？
慢和快是互相比较的，你不能吃了那和谁比的那半句啊。

[Mingci_why]

膜拜jackyjkchen大牛..

[gdstcwl]

找本数值计算的书看看,牛顿迭代什么的

[saramand9]

LZ 可以尝试下二分答案，对N进行开放的话LOG（N）的复杂度，不知道能不能满足需求

[tian545882735]

现成的为什么不用呢？
头文件#include<cmath>中的sqrt函数

[downmooner]

先搜索下 手动整数开平方 方法！

很简单。按照那个写程序就可以了