Round 2 - Final round?
--------------------------
To date, 8 working algorithms for computing square-roots are as follows:
opt.1 (tformed) exp() algorithm
opt.2 (Sturm) looping algorithm
opt.3 (Salem) suggested common log algorithm
opt.4 (firstPerson) natural log algorithm
opt.5 (firstPerson) Taylor's series
opt.6 (firstPerson) Brute force
opt.7 (invisal) Inverse Square Root
opt.8 (Foamgum) Spirit MySQRT
NathanOliver's ambitious nth root of a number algorithm has a bug in line 17, which gave rise to infinite loops. firstPerson's Brute force algorithm is original and simplest and the speed may probably be improved if the first computed estimate of the square root of a number entered can be used in the for loop, instead of starting from 0. Precision may have to be sacrificed to increase speed of computation.
The Sqrt program below is used for comparison of performance of the algorithms. Numbers up to 30, 40, 50 and 60 digits data are entered, and the answers given by various algorithms are used for ranking. e.g. a 30-digit number would look like this, 123456789012345678901234567890.
Slow algorithms are skipped, and are considered out of the race.
This time around, the top two entries can match C++'s square root function to the best of my knowledge. However, I may be wrong. Basing on the results, I would rank -- without any bias -- the algorithms in terms of originality and performance as follows:
Foamgum's Spirit MySQRT algorithm (opt. 8) - First place.
firstPerson's Taylor series algorithm (opt.5) - Second place.
invisal's inverse Square Root algorithm (opt.7) - third place.
Foamgum's algorithm is original and answers for numbers up to 60 digits number are correct.
firstPerson's Taylor's series algorithm's answers for numbers up to 60 digits are correct.
invisal's Inverse Square root algorithm uses a magic hex constant. Its answers for numbers up to 30 digits are correct, and for numbers from 40 digits onwards are not correct. Perhaps invisal can debug the codes, and let us have the debugged version. But, if it can only handle float and not double numeric variable, can greater accuracy be achieved? I am curious to know the answer.
Other algorithms were either too slow, were disqualified or could not work (got bug/s).
Run the program below to confirm the above results yourself.
Please feel free to modify the codes to improve them, and kindly post modified parts for me to update my codes. Thanks.
// sqrt InputCtrl & sqrt YBM 22 Sep 2009
#include <windows.h> // WinApi header
#include <iostream>
#include <string>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <iomanip>
#include <stdio.h>
#define maxloop 5
using namespace std;
string spaz(int&,int,string);
string datainp(string,int);
int escx=0;
int curs=0;
int extended=0;
string concate = "";
double weight;
double x;
char *endptr;
string sweight;
void gotoxy( short x, short y ) {
HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE);
COORD position = { x, y };
SetConsoleCursorPosition( hStdout, position );
}
double sturmsqrt(double num) {
double mod=1;
double c=0;
for(int d=0; d<10000000; c+=mod, d++)
if(c*c>num) {
c-=mod;
mod/=10;
}
return c;
}
float invsqrt(float zx) {
float xhalf = 0.5f*zx;
int i = *(int*)&zx;
i = 0x5f3759df - (i>>1);
zx = *(float*)&i;
for(int k=0; k<5; k++)
zx*=(1.5f - xhalf*zx*zx);
return 1.0f/zx;
}
double taylor(double zx) { // Taylor series
double k = 1, lim2 = 168;
double term;
double x=log(zx)/2;
double y;
double sum = 1;
int i=1;
y=1; sum = 1;
do {
k *= i++; y *= x; term = 1.0*y/k; sum += term;
} while (i <= lim2);
return sum;
}
float SquareRootOf(double num) { // firstPerson's BRUTE FORCE SQRT Check for perfect square root
int i = 0;
for(i = 0;i < 10000 ; i++) {
if(i * i > num) break;
else if(i * i == num) return double(i);
}
i -= 2;
double j = i;
float precision = 0.001f;
while(j * j < num - precision )
j+=precision;
return j;
}
double MySQRT(double number, int *iterations) { // Foamgum's Spirit MySQRT
double b, c, e(1e-12), f;
if(number < 0)
return -1;
while(e > number) // For very small values the criterion for stopping needs to be smaller.
e /= 10;
b = (1.0 + number)/2; // First estimate of square root of number.
c = (b - number/b)/2; // Correction to give a better estimate.
while(c > e && (*iterations) < 1000) {
f = b; // For very small and very large values we will see
b -= c; // if the correction has any effect.
if(f == b) // Applying the correction made no difference to our answer.
return b;
c = (b - number/b)/2; // A revised correction.
++(*iterations); // A counter to measure the effectiveness of our algorithm.
}
return b;
}
int xcz=0;
int decf=0;
int ctr =0;
int cc=0;
int brk=0;
int len=0;
std::string xcon;
char inpg,xconcate;
int getch(void);
int getche(void);
int xc[maxloop]={5,5,8,5,5}; // xy coordinate starting from zero
int yc[maxloop]={11,11,7,17,19};
int lgth[maxloop];
int condi=3;
char efld[maxloop][25] = {" "," ","Number : "," "," "};
double e=2.7182818283;
int main() {
double num;
num=0;
double heyyou,whome;
brk=0;
system ("color 0a");
do {
system ("CLS");
cout << "\n\tEsc:Exit FINDING THE SQUAREROOT OF A NUMBER";
cout << "\n\tOnly numbers and 1 dec. allowed, (codes may not be bug-free) ";
cout << "\n\tNumbers upto 1 trillion seem okay for opt.2 ";
cout << "\n\tNumbers above 10 trillion seem okay for opt.1, 3, 4, 5";
cout << "\n\tUse at your own risk. ";
sweight=datainp(efld[condi],condi);
len=sweight.length();
concate="";
if (brk==1) {break;}
num=weight;
if (num >0) {
whome=num;
heyyou= exp((1.0/2.0)*log(whome));
// double exp(double arg); is the same as e (2.7182818283) raised to the argth power, identical with opt.4.
cout << "\n\t0. (C++ sqrt() Ans.: "<< sqrt(num);
cout << "\n\t1. (tformed)'s exp()sqrt algorithm. Ans.: "<< heyyou;
cout << "\n\t2. (Sturm)'s Sturmsqrt algorithm. Ans.: "<< sturmsqrt(num);
x = pow(10,log10(num)/2);
cout << "\n\t3. (Salem) Common log - pow(10,log10(num)/2). Ans.: " << x;
x = pow(e,log(num)/2);
cout << "\n\t4. (firstPerson) Natural log - pow(e,log(num)/2). Ans.: " << x;
x = taylor(num);
cout << "\n\t5. (firstPerson) Taylor Series Ans.: " << x;
if (num <= 1000000000) {
float(x) = SquareRootOf(num);
cout << "\n\t6. (firstPerson) Brute force Ans.: " << x;
}
else { cout << "\n\t6. (firstPerson) Brute force skip "; }
float fnum;
fnum=float(num);
x = invsqrt(fnum);
cout << "\n\t7. (invisal) Inv.SquareRoot, hex const.0x5f3759df.Ans.: " << x;
double a, b; // Foamgum
int count(0);
a = num;
b = MySQRT(a, &count);
cout << "\n\t8. (Foamgum) Spirit MySQRT Ans.: " << b;
}
else {cout <<"\x7";}
cout << "\n\n\t<Press Any to cont.> ";
getch();
} while ( brk==0);
cout << "\n\n\tBye. Thanks for trying out Homemade SQRT\n";
cout << "\n\tand Homemade InputCtrl (Feedback, pls.still under construction>\n\n";
}
// Data Entry subroutine by YBM
string datainp(string efield,int condi) {
int c;
redo:
extended=0;
curs=0;
decf=0;
//std::string concate = "";
gotoxy(xc[condi-1],yc[condi-1]);
cout << efld[condi-1];
// gotoxy(xc[condi-1]+lgth[condi-1],yc[condi-1]);
do
{
c = getch();
if (c==27) { brk=1;break; }
char inp=c; // assign char from ASCII code Backspace's ASCII code is 8
inpg==inp;
extended==0;
// cout << "ASCII code of " << inpg << "=" << c;
if (c == 0x00 || c == 0XE0 || c==224 || c==0) {
extended = getch();
// cout << " EXTENDED: " << extended << "XXX";
cout << "\x7";
if (extended) {
c= extended;
inp=(char)NULL;
inpg==inp;
}
}
if (condi==3) { // Only Numbers and decimal ASCII 46 allowed
if (extended==0 && (((decf <=1) && ((c >=48 && c <=57) || c==46)) || ((c >=48 && c <=57) || c==46) && (decf==0) && c !=8)) {
if ((decf==0) || (decf >=1 && c != 46)) {
concate += inp; // string concatenation.. The character isn't extended
cout << inp;
}
if (c==46 && decf==0) { ++decf;}
}
else { concate=spaz(condi,c,concate);if (curs==1) {}}
}
} while (c != 13); // ascii 13 is <enter> or <carriage return>
if (c !=27) {
int len = concate.length();
// cout << "\n String length =" << len << " chars" << endl;
if (condi==3) {
// convert string to numeric variable
char xxc[len]; // intialise array
concate.copy(xxc,len); // fills in array with characters from the string
// xxc=concate.c_str();
weight = strtod(xxc,&endptr); // weight declared as global double
}
}
return concate;
} // datainput End
// SUBROUTINE TO MOVE CURSOR POSITION Up/Dn/LF/RT/PgUp/PgDn by YBM
string spaz(int& condi,int xcz,string xconcate) {
if (extended >0) {
cout << "\x7";
if (xconcate.length()==0) {
cout << "\x08";
cout << " ";}
// gotoxy(xc[condi]+lgth[condi],yc[condi]);
if (xcz==72 || xcz==73 || xcz==75 || xcz == 77 || xcz == 80 || xcz==81) {
cout << "\x7";
if (xcz==72 || xcz==77 || xcz==73) {
curs=1;
if (condi >1) {
}
}
if (xcz==80 || xcz==75 || xcz==81) {
curs=1;
if (condi < maxloop) {
}
}
extended=0;
}
extended=0;
}
else
{
if (xcz !=8) { // not backspace
}
else
{
len = xconcate.length();
if (len >0) {
if (condi==3) {
if (decf >=1 && xconcate.substr(len-1,1)=="."){
--decf;
}
}
cout << "\x08";
cout << " ";
cout << "\x08";
xconcate=xconcate.substr(0,len-1);
}
else
{
cout << "";
decf=0;
}
}
}
return xconcate;
} // spaz End