# Benchmarking Square root functions [closed]

Here's a follow-up to my posts about square roots (Studies on Square Roots and More Square Root). I've tidied up the code, but there's still a lot to do according to the feedback of the other posts.

The main thing I did here was to include the two source codes from IBM and SUN's libc as per the user qwr's suggestion. I adapted the code for C++ and included them in the benchmarks of my code. I took the inroot from the root.tbl file. For EMULV, I had to make an adaptation, and it was found in the dla.h file. The math_force_eval was in the math-barriers.h file. For the SUN code, only one adaptation related to C++ binaries was made. In all the codes, I used a beautifier to improve readability. Here are the files for compilation, considering that VS2022 (New Project - Console Application) was used.

# Code for review

### root.cpp

#include <iomanip>
#include <cmath>
#include <chrono>
#include <limits>
#include <vector>
#include <random>
#include <algorithm>
#include <numeric>
#include <iostream>
#include "ibm_sqrt.h"
#include "sun_sqrtf.h"

// Constants
const int MAX_ITERATIONS = 1000000;
const double precision = 1e-3; // Adjusted precision
bool showNewtonIteration = true;
bool showAnothersIteration = true;

struct IterationData {
double number;
double sqrtResult;
int iterations;
};

// Optimized Newton's method for finding the square root
double newtonSqrt(double number, double precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
double guess = number * 0.5; // Starting guess
double nextGuess = 0;
iterations = 0;

while (true) {
if (std::fabs(guess - nextGuess) < precision) {
break;
}
iterations++;
if (showIteration && showNewtonIteration) {
std::cout << "Newton's method - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << guess << " and more...";
if (std::fabs(guess - nextGuess) < precision) {
std::cout << " reached desired precision";
}
std::cout << std::endl;
}
nextGuess = guess;
guess = 0.5 * (guess + number / guess);
}
return nextGuess;
}

// Iterative another's method with dynamic initial approximation
double anothersMethod(double number, double precision, int& iterations, bool showIteration) {
double f = number * 0.5; // Using the initial approximation
double prev_f = 0;
iterations = 0;

while (true) {
// Check if the value converges within the given precision
if (std::fabs(f - prev_f) < precision) {
break;
}

iterations++;
if (showIteration && showAnothersIteration) {
std::cout << "Another's method - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << f << " and more...";
if (std::fabs(f - prev_f) < precision) {
std::cout << " reached desired precision";
}
std::cout << std::endl;
}
prev_f = f;
// New formula to improve the estimate
f = f - (f * f - number) / (f + number / f);
}
return f;
}

// Union to manipulate double as an integer (bit manipulation)
union {
double value;
uint64_t bits;
} db;

// Initial approximation using bit manipulation
double initialSqrtEstimate(double x) {
db.value = x;

db.bits = (db.bits & 0x800FFFFFFFFFFFFFULL) | (((db.bits & 0x7FF0000000000000ULL) >> 1) + 0x1FF0000000000000ULL);

return db.value;
}

union {
float value;
uint32_t bits;
} db3;

// Function to calculate the reciprocal using bit manipulation and Newton-Raphson method
float reciprocal(float x) {
db3.value = float(x);

// Initialize an estimate of the reciprocal of x using bit manipulation
db3.bits = 0x7eed4f50 - db3.bits;

// Refine the estimate using the Newton-Raphson method
float y = float(db3.value);
for (int i = 0; i < 2; i++) {
y = y * (2.0 - x * y); // Newton-Raphson iteration for reciprocal
}

return y;
}

// Another iterative method with dynamic initial approximation
double anothersMethod2(double number, double precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
int p = int(number) >> 1;
double guess = 0;
double nextGuess = initialSqrtEstimate(number) + (number / (1 << p));
iterations = 0;
while (true) {
if (std::fabs(nextGuess - guess) < precision) {
break;
}
iterations++;
if (showIteration && showAnothersIteration) {
std::cout << "Another's method 2 - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << nextGuess << " and more...";
if (std::fabs(nextGuess - guess) < precision) {
std::cout << " reached desired precision";
std::cout << std::endl;
break;
}
std::cout << std::endl;
}
guess = nextGuess;
nextGuess = 0.5 * (guess + number / guess);
}
return guess;
}

// Another iterative method using initial approximation and multiplication-only formula
double anothersMethod3(double number, double precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
double f = initialSqrtEstimate(number);
double prev_f = 0;
iterations = 0;

while (iterations < MAX_ITERATIONS) {

if (std::fabs(f - prev_f) < precision) {
break;
}
iterations++;
if (showIteration && showAnothersIteration) {
std::cout << "Another's method 3 - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << f << " and more...";
if (std::fabs(f - prev_f) < precision) {
std::cout << " reached desired precision";
}
std::cout << std::endl;
}

prev_f = f;
f = 0.5 * (f + number * reciprocal(f));
}

if (showIteration && std::fabs(f * f - number) > precision) {
std::cerr << "Warning: Convergence not achieved to desired precision\n";
}

return f;
}

union {
float value;
uint32_t bits;
} db2;

// Function to calculate the inverse square root using Moler-Walsh method
float inverseSqrt(float x) {
db2.value = x;
db2.bits = 0x5f3759df - (db2.bits >> 1);
return db2.value;
}

// Function to calculate the square root using the inverse square root
float molerwalshMethod(float number, float precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
float invSqrt = inverseSqrt(number);
float sqrt = 1 / invSqrt;
iterations = 1;

if (showIteration) {
std::cout << "Moler-Walsh method - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << sqrt << " and more...";
std::cout << " reached desired precision";
std::cout << std::endl;
}
return sqrt;
}

// Moler-Walsh Method with a mix of inverse square root and refinement
float molerwalshMethod2(float number, float precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
float invSqrt = inverseSqrt(number);
float sqrt = 1 / invSqrt;
float prev_sqrt = 0;
iterations = 0;

while (true) {
if (std::fabs(sqrt - prev_sqrt) < precision) {
break;
}
iterations++;
if (showIteration) {
std::cout << "Moler-Walsh method 2 - ";
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) << sqrt << " and more...";
if (std::fabs(sqrt - prev_sqrt) < precision) {
std::cout << " reached desired precision";
}
std::cout << std::endl;
}
prev_sqrt = sqrt;
sqrt = 0.5 * (sqrt + number / sqrt);
}
return sqrt;
}

// Function to calculate the error
double calculateError(double root, double number) {
double square = root * root;
return std::fabs(square - number);
}

void showIterations(double number, double& resultNewton, double* resultAnothers, double& resultMolerWalsh, double& resultMolerWalsh2, double& resultIBM, float& resultSun, int& iterationsNewton, int* iterationsAnothers, int& iterationsMolerWalsh, int& iterationsMolerWalsh2) {
resultNewton = newtonSqrt(number, precision, iterationsNewton, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Newton's method final result: " << resultNewton << " reached desired precision, Iterations: " << iterationsNewton << "\n\n";

resultAnothers[0] = anothersMethod(number, precision, iterationsAnothers[0], true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Another's method final result: " << resultAnothers[0] << " reached desired precision, Iterations: " << iterationsAnothers[0] << "\n\n";

resultAnothers[1] = anothersMethod2(number, precision, iterationsAnothers[1], true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Another's method 2 final result: " << resultAnothers[1] << " reached desired precision, Iterations: " << iterationsAnothers[1] << "\n\n";

resultAnothers[2] = anothersMethod3(number, precision, iterationsAnothers[2], true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Another's method 3 final result: " << resultAnothers[2] << " reached desired precision, Iterations: " << iterationsAnothers[2] << "\n\n";

resultMolerWalsh = molerwalshMethod(static_cast<float>(number), static_cast<float>(precision), iterationsMolerWalsh, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Moler-Walsh method final result: " << resultMolerWalsh << " reached desired precision, Iterations: " << iterationsMolerWalsh << "\n\n";

resultMolerWalsh2 = molerwalshMethod2(static_cast<float>(number), static_cast<float>(precision), iterationsMolerWalsh2, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Moler-Walsh method 2 final result: " << resultMolerWalsh2 << " reached desired precision, Iterations: " << iterationsMolerWalsh2 << "\n\n";

resultIBM = __ieee754_sqrt(number);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "IBM method final result: " << resultIBM << " reached desired precision\n\n";

resultSun = __ieee754_sqrtf(static_cast<float>(number));
std::cout << std::setprecision(std::numeric_limits<float>::max_digits10)
<< "Sun method final result: " << resultSun << " reached desired precision\n\n";
}

void measureTime(double* randomNumbers, std::vector<double>& timesSqrt, std::vector<double>& timesNewton, std::vector<double>* timesAnothers, std::vector<double>& timesMolerWalsh, std::vector<double>& timesMolerWalsh2, std::vector<double>& timesIBM, std::vector<double>& timesSun, std::vector<double>& timesDivision, std::vector<double>& timesMultiplication, std::vector<IterationData>& newtonData, std::vector<IterationData>* anothersData, std::vector<IterationData>& molerwalshData, std::vector<IterationData>& molerwalsh2Data, std::vector<IterationData>& ibmData, std::vector<IterationData>& sunData) {
// Measuring time for standard sqrt function
std::cout << "Calculating times...\n";
auto startSqrt = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
volatile double result = std::sqrt(randomNumbers[i]);
}
auto endSqrt = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationSqrt = endSqrt - startSqrt;
timesSqrt.push_back(durationSqrt.count());

// Measuring time for Newton's pure method
auto startNewton = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsNewton = 0;
double result = newtonSqrt(randomNumbers[i], precision, iterationsNewton, false);
newtonData.push_back({ randomNumbers[i], result, iterationsNewton });
}
auto endNewton = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationNewton = endNewton - startNewton;
timesNewton.push_back(durationNewton.count());

// Measuring time for the another's method
auto startAnothers = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsAnothers = 0;
double result = anothersMethod(randomNumbers[i], precision, iterationsAnothers, false);
anothersData[0].push_back({ randomNumbers[i], result, iterationsAnothers });
}
auto endAnothers = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationAnothers = endAnothers - startAnothers;
timesAnothers[0].push_back(durationAnothers.count());

// Measuring time for the another's method 2
auto startAnothers2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsAnothers = 0;
double result = anothersMethod2(randomNumbers[i], precision, iterationsAnothers, false);
anothersData[1].push_back({ randomNumbers[i], result, iterationsAnothers });
}
auto endAnothers2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationAnothers2 = endAnothers2 - startAnothers2;
timesAnothers[1].push_back(durationAnothers2.count());

// Measuring time for the another's method 3
auto startAnothers3 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsAnothers = 0;
double result = anothersMethod3(randomNumbers[i], precision, iterationsAnothers, false);
anothersData[2].push_back({ randomNumbers[i], result, iterationsAnothers });
}
auto endAnothers3 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationAnothers3 = endAnothers3 - startAnothers3;
timesAnothers[2].push_back(durationAnothers3.count());

// Measuring time for Moler-Walsh method
auto startMolerWalsh = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsMolerWalsh = 0;
float result = molerwalshMethod(static_cast<float>(randomNumbers[i]), static_cast<float>(precision), iterationsMolerWalsh, false);
molerwalshData.push_back({ randomNumbers[i], result, iterationsMolerWalsh });
}
auto endMolerWalsh = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationMolerWalsh = endMolerWalsh - startMolerWalsh;
timesMolerWalsh.push_back(durationMolerWalsh.count());

// Measuring time for Moler-Walsh method 2
auto startMolerWalsh2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsMolerWalsh2 = 0;
float result = molerwalshMethod2(static_cast<float>(randomNumbers[i]), static_cast<float>(precision), iterationsMolerWalsh2, false);
molerwalsh2Data.push_back({ randomNumbers[i], result, iterationsMolerWalsh2 });
}
auto endMolerWalsh2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationMolerWalsh2 = endMolerWalsh2 - startMolerWalsh2;
timesMolerWalsh2.push_back(durationMolerWalsh2.count());

// Measuring time for IBM method
auto startIBM = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
double result = __ieee754_sqrt(randomNumbers[i]);
ibmData.push_back({ randomNumbers[i], result, 1 });
}
auto endIBM = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationIBM = endIBM - startIBM;
timesIBM.push_back(durationIBM.count());

// Measuring time for Sun method
auto startSun = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
float result = __ieee754_sqrtf(static_cast<float>(randomNumbers[i]));
sunData.push_back({ randomNumbers[i], result, 1 });
}
auto endSun = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationSun = endSun - startSun;
timesSun.push_back(durationSun.count());

// Measuring time for division by 2
auto startDivision = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
volatile double result = 2.0 / randomNumbers[i];
}
auto endDivision = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationDivision = endDivision - startDivision;
timesDivision.push_back(durationDivision.count());

// Measuring time for multiplication by 2
auto startMultiplication = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
volatile double result = 2.0 * randomNumbers[i];
}
auto endMultiplication = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationMultiplication = endMultiplication - startMultiplication;
timesMultiplication.push_back(durationMultiplication.count());
}

void displayAverageTimes(const std::vector<double>& timesSqrt, const std::vector<double>& timesNewton, const std::vector<double>* timesAnothers, const std::vector<double>& timesMolerWalsh, const std::vector<double>& timesMolerWalsh2, const std::vector<double>& timesIBM, const std::vector<double>& timesSun, const std::vector<double>& timesDivision, const std::vector<double>& timesMultiplication) {
double avgSqrt = std::accumulate(timesSqrt.begin(), timesSqrt.end(), 0.0) / timesSqrt.size();
double avgNewton = std::accumulate(timesNewton.begin(), timesNewton.end(), 0.0) / timesNewton.size();
double avgAnothers = std::accumulate(timesAnothers[0].begin(), timesAnothers[0].end(), 0.0) / timesAnothers[0].size();
double avgAnothers2 = std::accumulate(timesAnothers[1].begin(), timesAnothers[1].end(), 0.0) / timesAnothers[1].size();
double avgAnothers3 = std::accumulate(timesAnothers[2].begin(), timesAnothers[2].end(), 0.0) / timesAnothers[2].size();
double avgMolerWalsh = std::accumulate(timesMolerWalsh.begin(), timesMolerWalsh.end(), 0.0) / timesMolerWalsh.size();
double avgMolerWalsh2 = std::accumulate(timesMolerWalsh2.begin(), timesMolerWalsh2.end(), 0.0) / timesMolerWalsh2.size();
double avgIBM = std::accumulate(timesIBM.begin(), timesIBM.end(), 0.0) / timesIBM.size();
double avgSun = std::accumulate(timesSun.begin(), timesSun.end(), 0.0) / timesSun.size();
double avgDivision = std::accumulate(timesDivision.begin(), timesDivision.end(), 0.0) / timesDivision.size();
double avgMultiplication = std::accumulate(timesMultiplication.begin(), timesMultiplication.end(), 0.0) / timesMultiplication.size();

std::cout << "\nAverage Results:\n";
std::cout << std::scientific << std::setprecision(6); // Ensure values are displayed in scientific notation
std::cout << "Standard sqrt average time (ms): " << avgSqrt << std::endl;
std::cout << "Newton's method average time (ms): " << avgNewton << std::endl;
std::cout << "Another's method average time (ms): " << avgAnothers << std::endl;
std::cout << "Another's method 2 average time (ms): " << avgAnothers2 << std::endl;
std::cout << "Another's method 3 average time (ms): " << avgAnothers3 << std::endl;
std::cout << "Moler-Walsh method average time (ms): " << avgMolerWalsh << std::endl;
std::cout << "Moler-Walsh method 2 average time (ms): " << avgMolerWalsh2 << std::endl;
std::cout << "IBM method average time (ms): " << avgIBM << std::endl;
std::cout << "Sun method average time (ms): " << avgSun << std::endl;
std::cout << "Division by 2 average time (ms): " << avgDivision << std::endl;
std::cout << "Multiplication by 2 average time (ms): " << avgMultiplication << std::endl;
}

void calculateAndDisplayErrors(double number, const std::vector<double>& roots) {
double actualSqrt = std::sqrt(number);
std::cout << "\nError comparison with actual sqrt(" << number << "): " << std::setprecision(std::numeric_limits<double>::max_digits10) << actualSqrt << "\n";
for (size_t i = 0; i < roots.size(); ++i) {
double error = calculateError(roots[i], number);
std::cout << "Error for root " << i + 1 << " (" << std::scientific << std::setprecision(std::numeric_limits<double>::max_digits10) << roots[i] << "): " << error << "\n";
}
}

void displaySomeIterations(const std::vector<IterationData>& newtonData, const std::vector<IterationData>& anothersData1, const std::vector<IterationData>& anothersData2, const std::vector<IterationData>& anothersData3, const std::vector<IterationData>& molerwalshData, const std::vector<IterationData>& molerwalsh2Data, const std::vector<IterationData>& ibmData, const std::vector<IterationData>& sunData) {
std::vector<IterationData> sortedNewtonData = newtonData;
std::vector<IterationData> sortedAnothersData1 = anothersData1;
std::vector<IterationData> sortedAnothersData2 = anothersData2;
std::vector<IterationData> sortedAnothersData3 = anothersData3;
std::vector<IterationData> sortedMolerWalshData = molerwalshData;
std::vector<IterationData> sortedMolerWalsh2Data = molerwalsh2Data;
std::vector<IterationData> sortedIbmData = ibmData;
std::vector<IterationData> sortedSunData = sunData;

// Sort the data based on the number of iterations
std::sort(sortedNewtonData.begin(), sortedNewtonData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedAnothersData1.begin(), sortedAnothersData1.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedAnothersData2.begin(), sortedAnothersData2.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedAnothersData3.begin(), sortedAnothersData3.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedMolerWalshData.begin(), sortedMolerWalshData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedMolerWalsh2Data.begin(), sortedMolerWalsh2Data.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedIbmData.begin(), sortedIbmData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedSunData.begin(), sortedSunData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::cout << "\nSome random numbers and their iterations for Newton's method:\n";
for (size_t i = 0; i < 10 && i < sortedNewtonData.size(); ++i) {
std::cout << "Number: " << sortedNewtonData[i].number
<< ", Iterations: " << sortedNewtonData[i].iterations
<< ", Sqrt: " << sortedNewtonData[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Another's method:\n";
for (size_t i = 0; i < 10 && i < sortedAnothersData1.size(); ++i) {
std::cout << "Number: " << sortedAnothersData1[i].number
<< ", Iterations: " << sortedAnothersData1[i].iterations
<< ", Sqrt: " << sortedAnothersData1[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Another's method 2:\n";
for (size_t i = 0; i < 10 && i < sortedAnothersData2.size(); ++i) {
std::cout << "Number: " << sortedAnothersData2[i].number
<< ", Iterations: " << sortedAnothersData2[i].iterations
<< ", Sqrt: " << sortedAnothersData2[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Another's method 3:\n";
for (size_t i = 0; i < 10 && i < sortedAnothersData3.size(); ++i) {
std::cout << "Number: " << sortedAnothersData3[i].number
<< ", Iterations: " << sortedAnothersData3[i].iterations
<< ", Sqrt: " << sortedAnothersData3[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Moler-Walsh method:\n";
for (size_t i = 0; i < 10 && i < sortedMolerWalshData.size(); ++i) {
std::cout << "Number: " << sortedMolerWalshData[i].number
<< ", Iterations: " << sortedMolerWalshData[i].iterations
<< ", Sqrt: " << sortedMolerWalshData[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Moler-Walsh method 2:\n";
for (size_t i = 0; i < 10 && i < sortedMolerWalsh2Data.size(); ++i) {
std::cout << "Number: " << sortedMolerWalsh2Data[i].number
<< ", Iterations: " << sortedMolerWalsh2Data[i].iterations
<< ", Sqrt: " << sortedMolerWalsh2Data[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for IBM method:\n";
for (size_t i = 0; i < 10 && i < sortedIbmData.size(); ++i) {
std::cout << "Number: " << sortedIbmData[i].number
<< ", Sqrt: " << sortedIbmData[i].sqrtResult << "\n";
}

std::cout << "\nSome random numbers and their iterations for Sun method:\n";
for (size_t i = 0; i < 10 && i < sortedSunData.size(); ++i) {
std::cout << "Number: " << sortedSunData[i].number
<< ", Sqrt: " << sortedSunData[i].sqrtResult << "\n";
}
}

void clearScreen() {
#ifdef _WIN32
system("cls");
#else
system("clear");
#endif
}

int main() {
double number;

std::cout << "Enter the number to find the square root: ";
std::cin >> number;

// Show iterations for each method with user input
double resultNewton, resultAnothers[3], resultMolerWalsh, resultMolerWalsh2, resultIBM;
float resultSun;
int iterationsNewton, iterationsAnothers[3], iterationsMolerWalsh, iterationsMolerWalsh2;
showIterations(number, resultNewton, resultAnothers, resultMolerWalsh, resultMolerWalsh2, resultIBM, resultSun, iterationsNewton, iterationsAnothers, iterationsMolerWalsh, iterationsMolerWalsh2);

std::cout << "Press enter to continue . . .";
std::cin.get();
std::cin.get();

std::vector<double> timesSqrt, timesNewton, timesAnothers[3], timesMolerWalsh, timesMolerWalsh2, timesIBM, timesSun, timesDivision, timesMultiplication;
std::vector<IterationData> newtonData, anothersData[3], molerwalshData, molerwalsh2Data, ibmData, sunData;

double* randomNumbers = new double[MAX_ITERATIONS];
std::mt19937_64 rng;
std::uniform_real_distribution<double> dist(1.0, 100.0); // Random numbers between 1 and 100

for (int i = 0; i < MAX_ITERATIONS; i++) {
randomNumbers[i] = dist(rng);
}

measureTime(randomNumbers, timesSqrt, timesNewton, timesAnothers, timesMolerWalsh, timesMolerWalsh2, timesIBM, timesSun, timesDivision, timesMultiplication, newtonData, anothersData, molerwalshData, molerwalsh2Data, ibmData, sunData);

delete[] randomNumbers;

std::cout << "Press enter to continue to error calculations . . .";
std::cin.get();

// Display the error calculations
std::vector<double> roots = { resultNewton, resultAnothers[0], resultAnothers[1], resultAnothers[2], resultMolerWalsh, resultMolerWalsh2, resultIBM, static_cast<double>(resultSun) };
calculateAndDisplayErrors(number, roots);

std::cout << "Press enter to continue to time measurements . . .";
std::cin.get();

// Display the average times
displayAverageTimes(timesSqrt, timesNewton, timesAnothers, timesMolerWalsh, timesMolerWalsh2, timesIBM, timesSun, timesDivision, timesMultiplication);

std::cout << "Press enter to continue to top iterations display . . .";
std::cin.get();

// Display the top iterations
std::cout << "Some random numbers and their iterations:\n";
displaySomeIterations(newtonData, anothersData[0], anothersData[1], anothersData[2], molerwalshData, molerwalsh2Data, ibmData, sunData);

std::cout << "Press enter to exit . . .";
std::cin.get();

return 0;
}


## Reference implementations for benchmarks

(This is not my code, so not subject to review)

### ibm_sqrt.h

#ifndef IBM_SQRT_H
#define IBM_SQRT_H

#include <cmath>
#include <cfenv>
#include <cstdint>
#include <algorithm>
#include <fenv.h>
#include <vector>
#include <mutex>

double __ieee754_sqrt(double x);

#ifndef GENERATE
static const double inroot[128] = {
1.40872145012100,  1.39792649065766,  1.38737595123859,  1.37706074531819,
1.36697225234682,  1.35710228748795,  1.34744307370643,  1.33798721601135,
1.32872767765984,  1.31965775814772,  1.31077107283046,  1.30206153403386,
1.29352333352711,  1.28515092624400,  1.27693901514820,  1.26888253714903,
1.26097664998256,  1.25321671998073,  1.24559831065844,  1.23811717205462,
1.23076923076923,  1.22355058064300,  1.21645747403153,  1.20948631362953,
1.20263364480453,  1.19589614840310,  1.18927063399547,  1.18275403352732,
1.17634339535009,  1.17003587860341,  1.16382874792529,  1.15771936846787,
1.15170520119791,  1.14578379846309,  1.13995279980655,  1.13420992801334,
1.12855298537376,  1.12297985014975,  1.11748847323133,  1.11207687497107,
1.10674314218572,  1.10148542531442,  1.09630193572405,  1.09119094315276,
1.08615077328341,  1.08117980543918,  1.07627647039410,  1.07143924829188,
1.06666666666667,  1.06195729855996,  1.05730976072814,  1.05272271193563,
1.04819485132867,  1.04372491688551,  1.03931168393861,  1.03495396376504,
1.03065060224133,  1.02640047855933,  1.02220250399990,  1.01805562076124,
1.01395880083916,  1.00991104495649,  1.00591138153909,  1.00195886573624,
0.99611649018350,  0.98848330114434,  0.98102294317595,  0.97372899112030,
0.96659534932828,  0.95961623024651,  0.95278613468066,  0.94609983358253,
0.93955235122353,  0.93313894963169,  0.92685511418159,  0.92069654023750,
0.91465912076005,  0.90873893479530,  0.90293223677296,  0.89723544654727,
0.89164514012056,  0.88615804099474,  0.88077101210109,  0.87548104826333,
0.87028526915267,  0.86518091269740,  0.86016532891275,  0.85523597411976,
0.85039040552437,  0.84562627613070,  0.84094132996422,  0.83633339758291,
0.83180039185606,  0.82734030399203,  0.82295119979782,  0.81863121615464,
0.81437855769486,  0.81019149366693,  0.80606835497581,  0.80200753138734,
0.79800746888611,  0.79406666717674,  0.79018367731967,  0.78635709949278,
0.78258558087123,  0.77886781361798,  0.77520253297841,  0.77158851547266,
0.76802457717971,  0.76450957210799,  0.76104239064719,  0.75762195809661,
0.75424723326565,  0.75091720714229,  0.74763090162560,  0.74438736831878,
0.74118568737933,  0.73802496642311,  0.73490433947940,  0.73182296599416,
0.72878002987884,  0.72577473860242,  0.72280632232420,  0.71987403306536,
0.71697714391715,  0.71411494828392,  0.71128675915902,  0.70849190843208 };
#endif
#endif // IBM_SQRT_H


### sun_sqrtf.h

#ifndef SUN_SQRTF_H
#define SUN_SQRTF_H

float __ieee754_sqrtf(float x);

#endif // SUN_SQRTF_H


### ibm_sqrt.cpp

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001-2024 Free Software Foundation, Inc.
*
* This program is free software; you can redistribute it and/or modify
* the Free Software Foundation; either version 2.1 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, see <https://www.gnu.org/licenses/>.
*/
/*********************************************************************/
/* MODULE_NAME: uroot.c                                              */
/*                                                                   */
/* FUNCTION:    usqrt                                                */
/*                                                                   */
/* FILES NEEDED: endian.h mydefs.h dla.h                             */
/*               root.tbl                                            */
/*                                                                   */
/* An ultimate sqrt routine. Given an IEEE double machine number x   */
/* it computes the correctly rounded (to nearest) value of square    */
/* root of x.                                                        */
/* Assumption: Machine arithmetic operations are performed in        */
/* round to nearest mode of IEEE 754 standard.                       */
/*                                                                   */
/*********************************************************************/

#include "ibm_sqrt.h"

#define USE_SQRT_BUILTIN 0

/*#ifdef DLA_FMS
# define  EMULV(x, y, z, zz)          \
z = x * y; zz = DLA_FMS (x, y, z);
#else
# define  EMULV(x, y, z, zz)          \
({  __typeof__ (x) __p, hx, tx, hy, ty;          \
__p = CN * (x);  hx = ((x) - __p) + __p;  tx = (x) - hx; \
__p = CN * (y);  hy = ((y) - __p) + __p;  ty = (y) - hy; \
z = (x) * (y); zz = (((hx * hy - z) + hx * ty) + tx * hy) + tx * ty; \
})
#endif*/

#define EMULV(a, b, c, d) \
{ \
c = a * b; \
d = std::fma(a, b, -c); \
}

inline void math_force_eval(double x) {
volatile double y = x;
(void)y;
}

union mynumber {
double x;
uint32_t i[2];
};

#ifdef GENERATE
std::vector<double> generate_inroot_table() {
std::vector<double> table(128);
for (int i = 0; i < 128; ++i) {
double value = 1.0 + static_cast<double>(i) / 128.0;
table[i] = 1.0 / std::sqrt(value);
}
return table;
}

static std::vector<double> inroot = generate_inroot_table();
#endif
static const double rt0 = 9.99999999859990725855365213134618E-01;
static const double rt1 = 4.99999999495955425917856814202739E-01;
static const double rt2 = 3.75017500867345182581453026130850E-01;
static const double rt3 = 3.12523626554518656309172508769531E-01;
static const double big = 134217728.0;

double __ieee754_sqrt(double x) {

double y, t, del, res, res1, hy, z, zz, s;
mynumber a, c = { {0} };
int32_t k;

a.x = x;
k = a.i[1];
a.i[1] = (k & 0x001fffff) | 0x3fe00000;
t = inroot[(k & 0x001fffff) >> 14];
s = a.x;
/*----------------- 2^-1022  <= | x |< 2^1024  -----------------*/
if (k > 0x000fffff && k < 0x7ff00000) {
int rm = std::fegetround();
std::fenv_t env;
feholdexcept(&env);
fesetround(FE_TONEAREST);
double ret;

y = 1.0 - t * (t * s);
t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
c.i[1] = 0x20000000 + ((k & 0x7fe00000) >> 1);
y = t * s;
hy = (y + big) - big;
del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy));
res = y + del;
if (res == (res + 1.002 * ((y - res) + del)))
ret = res * c.x;
else {
res1 = res + 1.5 * ((y - res) + del);
EMULV(res, res1, z, zz);
res = ((((z - s) + zz) < 0) ? std::max(res, res1) : std::min(res, res1));
ret = res * c.x;
}
math_force_eval(ret);
fesetenv(&env);
double dret = x / ret;
if (dret != ret) {
double force_inexact = 1.0 / 3.0;
math_force_eval(force_inexact);
/* The square root is inexact, ret is the round-to-nearest
value which may need adjusting for other rounding
modes.  */
switch (rm) {
#ifdef FE_UPWARD
case FE_UPWARD:
if (dret > ret)
ret = (res + 0x1p-1022) * c.x;
break;
#endif
#ifdef FE_DOWNWARD
case FE_DOWNWARD:
#endif
#ifdef FE_TOWARDZERO
case FE_TOWARDZERO:
#endif
#if defined FE_DOWNWARD || defined FE_TOWARDZERO
if (dret < ret)
ret = (res - 0x1p-1022) * c.x;
break;
#endif
default:
break;
}
}
/* Otherwise (x / ret == ret), either the square root was exact or
the division was inexact.  */
return ret;
}
else {
if ((k & 0x7ff00000) == 0x7ff00000)
return x * x + x; // sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN
if (x == 0)
return x;       // sqrt(+0)=+0, sqrt(-0)=-0
if (k < 0)
return (x - x) / (x - x); // sqrt(-ve)=sNaN
return 0x1p-256 * __ieee754_sqrt(x * 0x1p512);
}
}


### sun_sqrtf.cpp

/* e_sqrtf.c -- float version of e_sqrt.c.
*/

/*
* ====================================================
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/

#include "sun_sqrtf.h"
#include <cmath>
#include <cstdint>

// Macros to extract and set float word
#define GET_FLOAT_WORD(i, d)  \
do {                      \
union {               \
float f;          \
uint32_t i;       \
} _u;                 \
_u.f = (d);           \
(i) = _u.i;           \
} while (0)

#define SET_FLOAT_WORD(d, i)  \
do {                      \
union {               \
float f;          \
uint32_t i;       \
} _u;                 \
_u.i = (i);           \
(d) = _u.f;           \
} while (0)

float __ieee754_sqrtf(float x) {
#if USE_SQRTF_BUILTIN
return __builtin_sqrtf(x);
#else
float z;
int32_t sign = (int)0x80000000;
int32_t ix, s, q, m, t, i;
uint32_t r;

GET_FLOAT_WORD(ix, x);

/* take care of Inf and NaN */
if ((ix & 0x7f800000) == 0x7f800000) {
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */
}
/* take care of zero */
if (ix <= 0) {
if ((ix & (~sign)) == 0) return x; /* sqrt(+-0) = +-0 */
else if (ix < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix >> 23);
if (m == 0) { /* subnormal x */
for (i = 0; (ix & 0x00800000) == 0; i++) ix <<= 1;
m -= i - 1;
}
m -= 127; /* unbias exponent */
ix = (ix & 0x007fffff) | 0x00800000;
if (m & 1) /* odd m, double x to make it even */
ix += ix;
m >>= 1; /* m = [m/2] */

/* generate sqrt(x) bit by bit */
ix += ix;
q = s = 0; /* q = sqrt(x) */
r = 0x01000000; /* r = moving bit from right to left */

while (r != 0) {
t = s + r;
if (t <= ix) {
s = t + r;
ix -= t;
q += r;
}
ix += ix;
r >>= 1;
}

/* use floating add to find out rounding direction */
if (ix != 0) {
z = 0x1p0 - 0x1.4484cp-100; /* trigger inexact flag. */
if (z >= 0x1p0) {
z = 0x1p0 + 0x1.4484cp-100;
if (z > 0x1p0)
q += 2;
else
q += (q & 1);
}
}
ix = (q >> 1) + 0x3f000000;
ix += (m << 23);
SET_FLOAT_WORD(z, ix);
return z;
#endif /* ! USE_SQRTF_BUILTIN */
}

• Both codes are not suitable for x86-64 architecture. I think this architecture is being optimized for gaming. Commented Jul 30 at 23:58
• (not my code, so not subject to review using a block quote would make this very apparent. (Tedious without a small tool.)) It may be better to put all 3rd party code ("another's method") in "block-quoted code blocks". Commented Jul 31 at 8:22
• Actually, this is for review too. Anyone who wants to explain the mannerisms and code frameworks of these two companies, and what hardware they were designed for, feel free. Commented Jul 31 at 8:43
• Commented Jul 31 at 8:47
• Well, I want to say that I learned a lot about programming business here in this site, there types of things that make feel here is not the best (not even close) to exchange... :) Commented Jul 31 at 9:07