一道有趣的算法问题

[ldj123ldj]

这是一个o(a*R)的做法。
可以证明F(a,b)=F(b,a); 对于X=1e+8,R=1e+9的数据，普通cpu运行，预计需要1500年。
使用CUDA在GPU服务器上运算速度相差在100的量级内，如何充分发挥GPU性能？
求教各位是否有合适算法解决此问题。
这是cuda的代码，可以如何改进程序？

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include <cuda_runtime.h>

#define MAX_Fac_num 1000
#define N 12400
#define N2 100000

__device__ int num[N2];
int result[N2];
__global__ void get_num3(int R,int X,int n,int m){
	//计算每个F(R,X),R<m*n.结果为输出的num数组的元素和，一共有m个元素。(m是总线程个数，n是每个线程计算的r个数）
	int idx,r,a,k,tmp;
	double b,c;
	int i,fac_num;
	int fac1,fac2;
	idx = threadIdx.x + blockIdx.x*blockDim.x;
	//idx = blockIdx.x;
	for(k=0,tmp=0,num[idx]=0;k<n;k++){	
		r=(idx+k*m+1); //r从[1-m],[m+1,2m],....[(n-1)*m+1,n*m]
		if(r<=R){
			for(i=1,fac_num=0;i<=2*r;i++){
				if((4*r*r)%i==0){
					fac_num++;
					fac1=i;
					fac2=4*r*r/i;
					for(a=1;a<=X;a++){
						b = 1.0/4*(fac2 + fac1 + a*1.0/r*(fac2 - fac1));
						c = 1.0/4*(fac2 + fac1 - a*1.0/r*(fac2 - fac1));
						if((fabs(c-round(c))<=1e-6)&&fabs(b-round(b))<=1e-6){
							num[idx]=num[idx]+4-2*(fabs(b-c)<0.01);//b=c加2，b~=c加4；
						}
					} 		
				}	
			}
			//printf("r=%d,fac_num=%d,num=%d\n",r,fac_num,num[r]);
			//对每个r,对所以x的计算结果和
		//	printf("blockIdx=%d,Idx:%d,\tS(%d,%d)=%d\n",blockIdx.x,idx,r,X,num[idx]-tmp);
		//	tmp=num[idx];//存储上一次结果
		//printf("Idx:%d,sum=%d\n",idx,num[idx]);
		}//if
		//else{
		//	printf("idx:%d,blockIdx=%d,threadIdx=%d,r=%d,num[idx]=%d\n",idx,blockIdx.x,threadIdx.x,r,num[idx]);
		//}
	}
}

int main(){
	int count=0;      //传回有计算能力的设备数(≥1)，没有回传回1，device 0是一个仿真装置，不支持CUDA功能
   	cudaGetDeviceCount(&count);
	if(count == 0){   //没有cuda计算能力的设备
        	fprintf(stderr,"There is no device.\n");
        	return 1;
	}
	printf("num of GPU=%d\n",count);
	int R=100,X=100,idx;
	int m,m2,n=5,sum=0;//n is num of r in each threadIdx.
	int thread_num=1;
	clock_t start,end1;
	for(R=10,X=10;R<=1e+5;R=R*10,X=X*10){
		start = clock();
		m2=R/n;//有效线程数
		printf("=======================================R=%d==============================================\n",R);
		m=(int)ceil(1.0*m2/(thread_num*1.));//线程块数目(m*n*thread_num>m2*n=R)
		printf("block_num=%d,thread_num=%d,all_thread_num=%d,n=%d\n",m,thread_num,m*thread_num,n);
		get_num3<<<m, thread_num >>>(R,X,n,m*thread_num);
		
		cudaMemcpyFromSymbol(result, num, N2*sizeof(int));		
		//对所有线程的结果求和
		for(idx=0,sum=0;idx<=m*thread_num; idx++){
			//printf("result[%d]=%d\n",idx,result[idx]);
			sum+=result[idx];
		}
		printf("F(%d,%d)=%d\n",R,X,sum);
		end1 = clock();
		printf("used time=%f\n",(double)(end1-start)/CLOCKS_PER_SEC);
	}
	return 0;	
} 

这是运行结果

==================================demo_1==============================================
gpu006

num of GPU=1
=======================================R=10==============================================
block_num=1,thread_num=32,all_thread_num=32,n=2
F(10,10)=340
used time=0.310000
=======================================R=100==============================================
block_num=2,thread_num=32,all_thread_num=64,n=2
F(100,100)=38412
used time=0.010000
=======================================R=1000==============================================
block_num=16,thread_num=32,all_thread_num=512,n=2
F(1000,1000)=3967260
used time=0.150000
=======================================R=10000==============================================
block_num=157,thread_num=32,all_thread_num=5024,n=2
F(10000,10000)=399499660
used time=2.850000



=====================================demo_2===============================================
gpu006

num of GPU=1
=======================================R=10==============================================
block_num=1,thread_num=512,all_thread_num=512,n=2
F(10,10)=340
used time=0.370000
=======================================R=100==============================================
block_num=1,thread_num=512,all_thread_num=512,n=2
F(100,100)=38412
used time=0.000000
=======================================R=1000==============================================
block_num=1,thread_num=512,all_thread_num=512,n=2
F(1000,1000)=3967260
used time=0.180000
=======================================R=10000==============================================
block_num=10,thread_num=512,all_thread_num=5120,n=2
F(10000,10000)=399499660
used time=3.450000
=======================================R=100000==============================================
block_num=98,thread_num=512,all_thread_num=50176,n=2
F(100000,100000)=800981660
used time=3078.800000


========================================demo_3==============================================
gpu006

num of GPU=1
=======================================R=10==============================================
block_num=2,thread_num=1,all_thread_num=2,n=5
F(10,10)=340
used time=0.400000
=======================================R=100==============================================
block_num=20,thread_num=1,all_thread_num=20,n=5
F(100,100)=38412
used time=0.000000
=======================================R=1000==============================================
block_num=200,thread_num=1,all_thread_num=200,n=5
F(1000,1000)=3967260
used time=0.060000
=======================================R=10000==============================================
block_num=2000,thread_num=1,all_thread_num=2000,n=5
F(10000,10000)=399499660
used time=1.590000


=============================================demo_4============================================
gpu006

num of GPU=1
=======================================R=10==============================================
block_num=1,thread_num=256,all_thread_num=256,n=5
F(10,10)=340
used time=0.410000
=======================================R=100==============================================
block_num=1,thread_num=256,all_thread_num=256,n=5
F(100,100)=38412
used time=0.010000
=======================================R=1000==============================================
block_num=1,thread_num=256,all_thread_num=256,n=5
F(1000,1000)=3967260
used time=0.250000
=======================================R=10000==============================================
block_num=8,thread_num=256,all_thread_num=2048,n=5
F(10000,10000)=399499660
used time=6.430000



============================================demo_5=========================================
gpu006

num of GPU=1
=======================================R=10==============================================
block_num=2,thread_num=1,all_thread_num=2,n=5
F(10,10)=340
used time=0.460000
=======================================R=100==============================================
block_num=20,thread_num=1,all_thread_num=20,n=5
F(100,100)=38412
used time=0.000000
=======================================R=1000==============================================
block_num=200,thread_num=1,all_thread_num=200,n=5
F(1000,1000)=3967260
used time=0.050000
=======================================R=10000==============================================
block_num=2000,thread_num=1,all_thread_num=2000,n=5
F(10000,10000)=399499660
used time=1.610000
=======================================R=100000==============================================
block_num=20000,thread_num=1,all_thread_num=20000,n=5
F(100000,100000)=800981660
used time=3106.660000