Algorithms

void test1(){
	const long MAX_Iter=10000000;//10^7
	int *ar1=new int [MAX_Iter];/// array of integers to find the reciprocal
	double *R=new double [MAX_Iter];//array of double for saving the regular division 1/x
	double *R2=new double [MAX_Iter];//array of double for saving the algorithms division 1/x
	double *ERROR=new double [MAX_Iter];
	double max_error=pow(10,-100.0);;

	/*
	*  generate some number bigger than 2
	*/
	for(int i=0;i<MAX_Iter;i++){
		ar1[i]=rand()+2+clock();
	}
	cout<<"===================================BESME Allah============================="<<endl<<endl;;
	cout<<"==========================regular algorithm 1/integer============================="<<endl;
	clock_t Start = clock();
	/*
	*regualr division 
	*/
	for(int i=0;i<MAX_Iter;i++){
		R[i]=1.0/ar1[i];
	}
	cout<<"				"<<clock()-Start<<endl;
	cout<<"==========================your suggestion algorithm (ASM):1/integer============================="<<endl;
	Start = clock();
	/*
	*suggested division (ASM)
	*/
	float inv;
	for(int i=0;i<MAX_Iter;i++){
		  //float inv=inverse(x);
		float x=ar1[i];
		unsigned int *y = (unsigned int*)&x;
		//*y=0x7EEEEBB3-*y;
		*y=0x7EF39252-*y;
		//inv=ar1[i];
		inv=x;
		  for(int j=0;j<3;j++){
			 inv=inv*(2-inv*ar1[i]);
		  }
		  R2[i]=inv;
	}
	cout<<"				"<<clock()-Start<<endl;
	max_error=pow(10,-100.0);
	for(int i=0;i<MAX_Iter;i++){
		ERROR[i]=abs(R[i]-R2[i]);
		if(max_error<ERROR[i])
			max_error=ERROR[i];
	}
	cout<<"					maximmum   error="<<max_error<<endl;
	cout<<"========================newton raphson division algorithm time=================="<<endl;
	Start = clock();
	/*
	*newton raphson division 
	*/
	int iter=15;
	for(int i=0;i<MAX_Iter;i++){
		//this part is for finding an approximate
		int j=32;
		double x=0.0000000004656612873077392578125;//2^(-32)
		while(!((ar1[i]<<(32-j))& 2147483648)){
			j--;
			x*=2;
		}
		//end of part of finding aproximate
		//newton raphson method
		for(int j=0;j<iter;j++){
			x=x*(2.0-ar1[i]*x);//x=x*(2-d*x);
		}
		R2[i]=x;
		//
	}
	cout<<"					"<<clock()-Start<<endl;
	max_error=pow(10,-100.0);
	for(int i=0;i<MAX_Iter;i++){
		ERROR[i]=abs(R[i]-R2[i]);
		if(max_error<ERROR[i])
			max_error=ERROR[i];
	}
	cout<<"					maximmum   error="<<max_error<<endl;
	cout<<"===========================goldschmidt division algorithm=========================="<<endl;
	/*
	*	goldschmidt division algorithm
	*/
	//
	#define BASE 31

	typedef unsigned int uint32;
	typedef unsigned long long int uint64;
	
	Start = clock();
//s	int iter=20;
	for(int i=0;i<MAX_Iter;i++){
		int shl = 0; // normalization shift
		uint32 nfp = ar1[i];//
		while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead

		uint64 q = 0x100000000ULL; // 2^32
		uint64 e = 0x100000000ULL - (uint64)nfp; // 2^32-NFP
		int j;
		for (j=0;j<5;j++) //for (i=0;i<100;i++) //for (i=0;i<4;i++) // iterate
		{
			// Both multiplications are actually
			// 32x32 bits truncated to the 32 high bits
			q += (q*e)>>(uint64)32;
			e = (e*e)>>(uint64)32;
			//printf("Q=0x%llx E=0x%llx\n",q,e);
		}
		// Here, (Q/2^32) is the inverse of (NFP/2^32).
		// We have 2^31<=NFP<2^32 and 2^32<Q<=2^33
		R2[i]=(uint32)(q>>(64-2*BASE-shl));
	}
	cout<<"					"<<clock()-Start<<endl;
	max_error=pow(10,-100.0);
	for(int i=0;i<MAX_Iter;i++){
		ERROR[i]=abs(R[i]-R2[i]);
		if(max_error<ERROR[i])
			max_error=ERROR[i];
	}
	cout<<"					maximmum    error="<<max_error<<endl;
}

void test3(){
	float x;
	float inv;
	cout<<endl;
	float main_x;
	while(cin>>x){
		main_x=x;
		unsigned int *y = (unsigned int*)&x;
		//*y=0x7EEEEBB3-*y;
		*y=0x7EF39252-*y;
		//inv=ar1[i];
		inv=x;
		  for(int j=0;j<3;j++){
			 inv=inv*(2-inv*main_x);
		  }
		  cout<<"1/x="<<(1.0/main_x)<<"   our inv="<<inv<<"     error="<<(1.0/main_x)-inv<<endl;
	}
}
void test4(){
	float x;
	float inv;
	cout<<endl;
	float main_x;
	while(cin>>x){
		main_x=x;
		//unsigned int *y = (unsigned int*)&x;
		//*y=0x7EF39252-*y;
		x=(0x7EF39252-(unsigned int)x);
		inv=x;
		/*  for(int j=0;j<3;j++){
			 inv=inv*(2-inv*main_x);
		  }*/
		  cout<<"1/x="<<(1.0/main_x)<<"   our inv="<<inv<<"     error="<<(1.0/main_x)-inv<<endl;
	}
}

mid = low +
      ((toFind - sortedArray[low]) * (high - low)) /
      (sortedArray[high] - sortedArray[low]);  //out of range is possible.

Alan Balkany wrote:

So the Wikipedia claim could be correct.