# Studies on Square Roots

This post marks the beginning of my ventures into calculating square roots. However, I have no idea how the processor, more specifically the FPU with the fsqrt instruction, performs the process. I have included a better method for measuring time and also made some improvements to the method I wanted to demonstrate in that post, achieving fewer iterations than the original. However, due to the many divisions, it is still a bit slower than the incredible Newton's method. Therefore, I will post two codes I developed during this time, for you to give your opinion and see what can be improved. Especially regarding the calculation of Carmack's method, I am not sure if there is another way to calculate the magic number without brute force, and I would like to know.

Here is the improved code from the initial post:

#include <iomanip>
#include <cmath>
#include <chrono>
#include <limits>
#include <vector>
#include <random>
#include <algorithm>
#include <numeric>

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

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

// Function to calculate a good initial approximation
double calculateInitial(double number) {
double initial = (number + 1) * 0.5; // Use multiplication instead of division
for (int i = 0; i < 5; i++) {
initial = 0.5 * (initial + number / initial);
}
return initial;
}

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

while (true) {
nextGuess = 0.5 * (guess + number / guess);
iterations++;
if (showIteration && showNewtonIteration) {
std::cout << "Newton's method - ";
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;
}
if (std::fabs(nextGuess - guess) < precision) {
break;
}
guess = nextGuess;
}
return nextGuess;
}

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

while (true) {
f = 0.5 * (f + number / f); // Recalculate initial dynamically
e = f - (f * f - number) / (f + number / f); // New formula

// Check if the value repeats within the precision
if (std::fabs(f - prev_f) < precision) {
break;
}
prev_f = f;
f = e;
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;
}
}
return f;
}

// 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, int &iterationsNewton, int &iterationsAnothers) {
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 = anothersMethod(number, precision, iterationsAnothers, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Another's method final result: " << resultAnothers << " reached desired precision, Iterations: " << iterationsAnothers << "\n\n";
}

void measureTime(double* randomNumbers, std::vector<double> &timesSqrt, std::vector<double> &timesNewton, std::vector<double> &timesAnothers,
std::vector<IterationData> &newtonData, std::vector<IterationData> &anothersData) {
// Measuring time for standard sqrt function
std::cout << "Calculating times...\n";
auto startSqrt = std::chrono::high_resolution_clock::now();
for (int i = 0; i < NUM_EXECUTIONS; 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 < NUM_EXECUTIONS; 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 < NUM_EXECUTIONS; i++) {
int iterationsAnothers = 0;
double result = anothersMethod(randomNumbers[i], precision, iterationsAnothers, false);
anothersData.push_back({randomNumbers[i], result, iterationsAnothers});
}
auto endAnothers = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationAnothers = endAnothers - startAnothers;
timesAnothers.push_back(durationAnothers.count());
}

void displayAverageTimes(const std::vector<double> &timesSqrt, const std::vector<double> &timesNewton, const std::vector<double> &timesAnothers) {
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.begin(), timesAnothers.end(), 0.0) / timesAnothers.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;
}

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::setprecision(std::numeric_limits<double>::max_digits10) << roots[i] << "): " << error << "\n";
}
}

void displaySomeIterations(const std::vector<IterationData> &newtonData, const std::vector<IterationData> &anothersData) {
std::vector<IterationData> sortedNewtonData = newtonData;
std::vector<IterationData> sortedAnothersData = anothersData;

// 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(sortedAnothersData.begin(), sortedAnothersData.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 < sortedAnothersData.size(); ++i) {
std::cout << "Number: " << sortedAnothersData[i].number
<< ", Iterations: " << sortedAnothersData[i].iterations
<< ", Sqrt: " << sortedAnothersData[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;
int iterationsNewton, iterationsAnothers;
showIterations(number, resultNewton, resultAnothers, iterationsNewton, iterationsAnothers);

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

std::vector<double> timesSqrt, timesNewton, timesAnothers;
std::vector<IterationData> newtonData, anothersData;

double* randomNumbers = new double[NUM_EXECUTIONS];
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 < NUM_EXECUTIONS; i++) {
randomNumbers[i] = dist(rng);
}

measureTime(randomNumbers, timesSqrt, timesNewton, timesAnothers, newtonData, anothersData);

delete[] randomNumbers;

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

// Display the error calculations
std::vector<double> roots = {resultNewton, resultAnothers};
calculateAndDisplayErrors(number, roots);

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

// Display the average times
displayAverageTimes(timesSqrt, timesNewton, timesAnothers);

std::cout << "Press any key 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);

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

return 0;
}


Here is the code to calculate the magic number:

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <vector>
#include <mutex>
#include <chrono>
#include <atomic>
#include <sstream>

// Uncomment the following line to compile for Linux
// #define LINUX

#ifdef LINUX
#include <unistd.h>
#define MOVE_CURSOR_TO(x, y) std::cout << "\033[" << (x) << ";" << (y) << "H"
#else
#include <windows.h>
#include <conio.h>
void moveCursorTo(int x, int y) {
COORD coord = {static_cast<SHORT>(y), static_cast<SHORT>(x)};
SetConsoleCursorPosition(GetStdHandle(STD_OUTPUT_HANDLE), coord);
}
#define MOVE_CURSOR_TO(x, y) moveCursorTo((x), (y))
#endif

const int BAR_WIDTH = 70; // Width of the progress bar
const int SLEEP_DURATION = 1; // Duration to sleep between progress updates in seconds
std::mutex mtx;
float globalMinError = std::numeric_limits<float>::max();
unsigned int globalBestMagicNumber = 0;
std::atomic<bool> done(false);

// Function to calculate the inverse square root using Carmack's method with an adjustable magic number
float fastInverseSqrt(float number, unsigned int magicNumber) {
float x2 = number * 0.5F;
float y = number;
union {
float f;
int i;
} u;
u.f = y;
u.i = magicNumber - (u.i >> 1);
y = u.f;
y = y * (1.5F - (x2 * y * y));
return y;
}

// Function to calculate the square root using the inverse square root
float calculateSqrt(float number, unsigned int magicNumber) {
float invSqrt = fastInverseSqrt(number, magicNumber);
return number * invSqrt;
}

// Function to display a single progress bar
void showProgressBar(unsigned int current, unsigned int total) {
float progress = static_cast<float>(current) / total;
int pos = static_cast<int>(BAR_WIDTH * progress);
std::ostringstream bar;
bar << "[";
for (int i = 0; i < BAR_WIDTH; ++i) {
if (i < pos) bar << "=";
else if (i == pos) bar << ">";
else bar << " ";
}
bar << "] " << std::setw(3) << static_cast<int>(progress * 100.0) << " %";
std::cout << bar.str();
}

// Function to test a range of magic numbers and find the best one
void testMagicNumbersRange(unsigned int start, unsigned int end, const std::vector<float>& testValues, unsigned int threadIndex) {
unsigned int totalIterations = end - start;
for (unsigned int magicNumber = start; magicNumber < end; ++magicNumber) {
float totalError = 0.0f;
for (float value : testValues) {
float approxSqrt = calculateSqrt(value, magicNumber);
float exactSqrt = std::sqrt(value);
float error = std::fabs(approxSqrt - exactSqrt);
totalError += error;
}
float averageError = totalError / testValues.size();

{
std::lock_guard<std::mutex> guard(mtx);
if (averageError < globalMinError) {
globalMinError = averageError;
globalBestMagicNumber = magicNumber;
}
}

// Update progress

// Show progress for this thread
if ((magicNumber - start) % 0x10000 == 0) {
std::lock_guard<std::mutex> guard(mtx);
MOVE_CURSOR_TO(threadIndex + 8, 1); // Move cursor to the appropriate line
std::cout.flush();
}
}
}

// Function to display the header
std::cout << "=========================================\n";
std::cout << " Inverse Square Root Calculation Progress\n";
std::cout << "=========================================\n";
std::cout << "Calculating the optimal magic number for\n";
std::cout << "fast inverse square root approximation.\n";
std::cout << "-----------------------------------------\n";
}

// Function to display overall progress
void showProgress(unsigned int totalIterations) {
while (!done) {

// Display progress for all threads
std::lock_guard<std::mutex> guard(mtx);
for (unsigned int i = 0; i < numThreads; ++i) {
MOVE_CURSOR_TO(i + 8, 1); // Move cursor to the appropriate line
std::cout << "Thread " << i + 1 << ": ";
std::cout << std::endl;
}
}
}

int main() {
auto start = std::chrono::high_resolution_clock::now(); // Start timing

std::vector<float> testValues = { 0.5f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 10.0f, 100.0f, 1000.0f };
//unsigned int startMagicNumber = 0x00000000; // Begin with no idea what it can be
//unsigned int endMagicNumber = 0xFFFFFFFF; // End with no idea what it can be
unsigned int startMagicNumber = 0x5F000000; // Approximately close to the real value
unsigned int endMagicNumber = 0x5FFFFFFF; // Approximately close to the real value
unsigned int rangeSize = (endMagicNumber - startMagicNumber) / numThreads;

unsigned int totalIterations = rangeSize;

for (unsigned int i = 0; i < numThreads; ++i) {
unsigned int start = startMagicNumber + i * rangeSize;
unsigned int end = (i == numThreads - 1) ? endMagicNumber : start + rangeSize;
}

for (auto& t : threads) {
t.join();
}

done = true;

auto end = std::chrono::high_resolution_clock::now(); // End timing
std::chrono::duration<double> duration = end - start; // Calculate duration

std::cout << "\nBest magic number: 0x" << std::hex << globalBestMagicNumber << std::dec << std::endl;

// Test the best magic number found
for (float value : testValues) {
float approxSqrt = calculateSqrt(value, globalBestMagicNumber);
float exactSqrt = std::sqrt(value);
std::cout << "Value: " << value
<< ", Approximated sqrt: " << approxSqrt
<< ", Exact sqrt: " << exactSqrt
<< ", Error: " << std::fabs(approxSqrt - exactSqrt) << std::endl;
}

std::cout << "\nTotal execution time: " << duration.count() << " seconds" << std::endl;
std::cout << "Press any key to exit . . .";
std::cin.get();

return 0;
}

• "I am not sure if there is another way to calculate the magic number without brute force, and I would like to know." There are a few different ways to come to a solution for the "magic number", including a relatively recent one by Moroz et al. in 2016 (arxiv.org/abs/1603.04483). The original one was probably found by brute force search over the 32 bit floats. That's what Lomont 2003 suggests, see 0x5f37642f.com for some links and context on the FISR. Commented Jul 27 at 21:37
• Please don't change the code after answers have been given, as this would invalidate those answers. See also "What should I do when someone answers my question?" If you have improved the code and wish for the new version to be reviewed, feel free to post a new question here on Code Review. Commented Jul 28 at 8:15
• @G.Sliepen - Can I add the other code here, without removing the old one? I would like to keep everything on the same page. In fact, I made at least two more versions of this "another's method"... I would like to post the last one that you removed and replaced the old one. In fact, I posted this because I wanted to insist on this formula that I created but which doesn't have much relation to the final result. If you let me I can include a function "anothers2" and comment that it was added later, or I can add it alone in a new code frame... And anyway, you can edit the answers too... Commented Jul 28 at 9:25
• @vanzuita That's not really a good idea. Posting it as a new question allows us to review the new code, without things getting confusing about which answer talks about which version of the code. Creating follow-up questions is totally fine. Once you have created a new question, it is OK to edit this question and to add a link to the new one. Commented Jul 28 at 13:08
• The fast inverse square root is only faster on systems without fast fp squareroot
– qwr
Commented Jul 29 at 15:10

Iteration limit

std::fabs(nextGuess - guess) < precision takes the floating out of floating point objects.

This iteration test make sense if floating point values were linearly distributed. Overall, floating point values are distributed logarithmically over its type range. There are about as many values in the [0.0001 0.0002) range as [1000.0 2000.0).

Try finding the square roots of 2.0e-100, 2.0, 2.0e+100. I'd expect answers like 1.41421356237e-50, 1.41421356237 , 1.41421356237e+50.

Instead a test, something like fabs(nextGuess - guess) < precision*nextGuess makes more sense.

• Advanced: IMO, the quality of a floating point function should assess the number of floating point values between it and the mathematical best answer.

Test values

testValues[] deserves to include 0.0, the smallest possible double and the largest finite value.

# Square root algorithms and hardware implementations

Almost any software or hardware implementation of (inverse) square roots, but also for trigonometric functions, can be split into several steps:

1. Range reduction. This maps the input to a small number which is easier to work with. For square roots, you can multiply the input by a power of four, and compensate by dividing the result by a corresponding power of two.
2. A rough estimate, which can be done in various ways. For floating point numbers, which are comprised of a mantissa and an exponent, a simple way is to just divide the exponent by two. Or if it's in a small range, for example by using range reduction first, you could do a table lookup.
3. Iterate using Newton-Raphson or some other method until the desired precision is reached.

Some CPU architectures have separate instructions for the initial guess and for the iteration step, allowing you to choose a trade-off between performance and accuracy.

In any case, you should be able to use a fixed number of steps to get the desired precision, without having to check the difference between two iterations.

# Use constexpr for constants

Use constexpr or even static constexpr when declaring constants.

# Printing

Some improvements can be made in the way you print results:

• Use '\n' instead of std::endl.
• Be consistent and end every line with \n. It should almost never be necessary to start a print statement with a \n.
• Use std::format() to format the output, or even better, use std::print() if possible.
• Don't use system() to clear the screen; it's not portable as you've already noticed, it is very inefficient, and depending on your environment variabels and shell configuration, "clear" or "cls" might not even do what you think. I recommend you just don't clear the screen; it's not essential for this code.
• Printing will take some amount of time, so it's not a good idea to do that in code which you are trying to benchmark.

# Avoid code duplication

While you created some helper functions, there is still some code duplication that could be avoided by creating a few more. For example, displaySomeIterations() basically contains the same code twice, once for Newton's methond, once for Another's method. In displayAverageTimes(), you are calculating the average for three different vectors. Even if you can do this in a one-liner, I would just create a function to return the average of a vector.

You can also create some code to do timing. Either a function that takes another function as input, or you could create a class that keeps track of elapsed time.

# Use std::chrono::steady_clock

Unfortunately, there is no guarantee that std::chrono::high_resolution_clock is a steady clock that does not suffer from time jumps (for example because of NTP updates, daylight saving time changes and other events). Instead, prefer to use std::chrono::steady_clock; it doesn't suffer from these problems, and it often has the same resolution as std::chrono::high_resolution_clock anyway.

For reference, this is how glibc e_sqrt.c implements correctly-rounded sqrt when there's no hardware builtin (a modern 64-bit x86 CPU will use something like a single SSE2 SQRTSD instruction). From SO question on how the function works, it's really not documented well - apparently Newton's method with range optimizations and some bit trickery (but no bit casting magic of the fast invsqrt).

/*
* 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: dla.h endian.h mydefs.h                             */
/*               uroot.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 "endian.h"
#include "mydefs.h"
#include <dla.h>
#include "root.tbl"
#include <math-barriers.h>
#include <math_private.h>
#include <fenv_private.h>
#include <libm-alias-finite.h>
#include <math-use-builtins.h>

/*********************************************************************/
/* An ultimate sqrt routine. Given an IEEE double machine number x   */
/* it computes the correctly rounded (to nearest) value of square    */
/* root of x.                                                        */
/*********************************************************************/
double
__ieee754_sqrt (double x)
{
#if USE_SQRT_BUILTIN
return __builtin_sqrt (x);
#else
/* Use generic implementation.  */
static const double
rt0 = 9.99999999859990725855365213134618E-01,
rt1 = 4.99999999495955425917856814202739E-01,
rt2 = 3.75017500867345182581453026130850E-01,
rt3 = 3.12523626554518656309172508769531E-01;
static const double big = 134217728.0;
double y, t, del, res, res1, hy, z, zz, s;
mynumber a, c = { { 0, 0 } };
int4 k;

a.x = x;
k = a.i[HIGH_HALF];
a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000;
t = inroot[(k & 0x001fffff) >> 14];
s = a.x;
/*----------------- 2^-1022  <= | x |< 2^1024  -----------------*/
if (k > 0x000fffff && k < 0x7ff00000)
{
int rm = __fegetround ();
fenv_t env;
libc_feholdexcept_setround (&env, FE_TONEAREST);
double ret;
y = 1.0 - t * (t * s);
t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
c.i[HIGH_HALF] = 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); /* (z+zz)=res*res1 */
res = ((((z - s) + zz) < 0) ? max (res, res1) :
min (res, res1));
ret = res * c.x;
}
math_force_eval (ret);
libc_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);
}
#endif /* ! USE_SQRT_BUILTIN  */
}
#ifndef __ieee754_sqrt
libm_alias_finite (__ieee754_sqrt, __sqrt)
#endif


Interestingly, the flt-32 version is written by a completely different author and is simpler and better-documented.

/* 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 <math.h>
#include <math_private.h>
#include <libm-alias-finite.h>
#include <math-use-builtins.h>

float
__ieee754_sqrtf(float x)
{
#if USE_SQRTF_BUILTIN
return __builtin_sqrtf (x);
#else
/* Use generic implementation.  */
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  */
}
#ifndef __ieee754_sqrtf
libm_alias_finite (__ieee754_sqrtf, __sqrtf)
#endif
$$$$

• One could benchmark this. Commented Jul 29 at 15:35
• @vanzuita you could, but only for academic purposes as a modern x86 CPU will use a single SQRTSD instruction.
– qwr
Commented Jul 29 at 16:02
• The 32 bit version is from fdlibm (netlib.org/fdlibm), which is interestingly also where the original idea for the bit fiddling in the FISR comes from (see the comment block in netlib.org/fdlibm/e_sqrt.c). Some of the code in fdlibm, including the sqrt code, was derived from the libm released for BSD 4.3 in 1986. Commented Jul 29 at 23:02
• Ok, you are saying it's not Carmack's method, but at least he made it very knowledgeble by releasing his source code. Commented Jul 30 at 15:52
• @vanzuita It never was Carmack's, though people attributed it to him because it leaked from the Q3A source code and sometimes people forget that games are developed by more than one person. ;) the source code for Q3A was formally released and open sourced, but that was a few years (2005) after the FISR code got super popular online. Commented Jul 30 at 16:14

Oh, this site could be more fluid. I think this thing slows down on multiplications and speeds up on divisions. This why it is nothing today. It speed very much on fsqrt, so I don't know this is very well documented, but if the new ones (my is old) do all things in all speed, may you will see the real difference between they. Including bit shift.

#include <iomanip>
#include <cmath>
#include <chrono>
#include <limits>
#include <vector>
#include <random>
#include <algorithm>
#include <numeric>
#include <iostream>

// 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;

// Bit manipulation to get a rough estimate of the square root
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); // Using the initial approximation
double prev_f = 0;
iterations = 0;

while (iterations < MAX_ITERATIONS) {
// 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 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)); // Use multiplication for refinement
}

// Ensure the final value meets the precision requirement
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 Carmack's method
float inverseSqrt(float x) {
db2.value = x;
db2.bits = 0x5f3759df - (db2.bits >> 1); // Initial approximation using magic number
return db2.value;
}

// Function to calculate the square root using the inverse square root
float carmacksMethod(float number, float precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
float invSqrt = inverseSqrt(number); // Get the inverse square root
float sqrt = 1 / invSqrt; // Calculate the square root by taking the reciprocal of the inverse
iterations = 1;

if (showIteration) {
std::cout << "Carmack's 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;
}

// Carmack's Method with a mix of inverse square root and refinement
float carmacksMethod2(float number, float precision, int& iterations, bool showIteration) {
if (number == 0) return 0;
float invSqrt = inverseSqrt(number); // Get the inverse square root
float sqrt = 1 / invSqrt; // Calculate the square root by taking the reciprocal of the inverse
float prev_sqrt = 0;
iterations = 0;

while (true) {
if (std::fabs(sqrt - prev_sqrt) < precision) {
break;
}
iterations++;
if (showIteration) {
std::cout << "Carmack's 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& resultCarmack, double& resultCarmack2, int& iterationsNewton, int* iterationsAnothers, int& iterationsCarmack, int& iterationsCarmack2) {
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";

resultCarmack = carmacksMethod(number, precision, iterationsCarmack, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Carmack's method final result: " << resultCarmack << " reached desired precision, Iterations: " << iterationsCarmack << "\n\n";

resultCarmack2 = carmacksMethod2(number, precision, iterationsCarmack2, true);
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
<< "Carmack's method 2 final result: " << resultCarmack2 << " reached desired precision, Iterations: " << iterationsCarmack2 << "\n\n";
}

void measureTime(double* randomNumbers, std::vector<double>& timesSqrt, std::vector<double>& timesNewton, std::vector<double>* timesAnothers, std::vector<double>& timesCarmack, std::vector<double>& timesCarmack2, std::vector<double>& timesDivision, std::vector<double>& timesMultiplication, std::vector<IterationData>& newtonData, std::vector<IterationData>* anothersData, std::vector<IterationData>& carmackData, std::vector<IterationData>& carmack2Data) {
// 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 Carmack's method
auto startCarmack = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsCarmack = 0;
double result = carmacksMethod(randomNumbers[i], precision, iterationsCarmack, false);
carmackData.push_back({ randomNumbers[i], result, iterationsCarmack });
}
auto endCarmack = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationCarmack = endCarmack - startCarmack;
timesCarmack.push_back(durationCarmack.count());

// Measuring time for Carmack's method 2
auto startCarmack2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsCarmack2 = 0;
double result = carmacksMethod2(randomNumbers[i], precision, iterationsCarmack2, false);
carmack2Data.push_back({ randomNumbers[i], result, iterationsCarmack2 });
}
auto endCarmack2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationCarmack2 = endCarmack2 - startCarmack2;
timesCarmack2.push_back(durationCarmack2.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>& timesCarmack, const std::vector<double>& timesCarmack2, 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 avgCarmack = std::accumulate(timesCarmack.begin(), timesCarmack.end(), 0.0) / timesCarmack.size();
double avgCarmack2 = std::accumulate(timesCarmack2.begin(), timesCarmack2.end(), 0.0) / timesCarmack2.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 << "Carmack's method average time (ms): " << avgCarmack << std::endl;
std::cout << "Carmack's method 2 average time (ms): " << avgCarmack2 << 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>* anothersData, const std::vector<IterationData>& carmackData, const std::vector<IterationData>& carmack2Data) {
std::vector<IterationData> sortedNewtonData = newtonData;
std::vector<IterationData> sortedAnothersData = anothersData[0];
std::vector<IterationData> sortedAnothersData1 = anothersData[1];
std::vector<IterationData> sortedAnothersData2 = anothersData[2];
std::vector<IterationData> sortedCarmackData = carmackData;
std::vector<IterationData> sortedCarmack2Data = carmack2Data;

// 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(sortedAnothersData.begin(), sortedAnothersData.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(sortedCarmackData.begin(), sortedCarmackData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedCarmack2Data.begin(), sortedCarmack2Data.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 < sortedAnothersData.size(); ++i) {
std::cout << "Number: " << sortedAnothersData[0, i].number
<< ", Iterations: " << sortedAnothersData[0, i].iterations
<< ", Sqrt: " << sortedAnothersData[0, 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 < sortedAnothersData.size(); ++i) {
std::cout << "Number: " << sortedAnothersData[1, i].number
<< ", Iterations: " << sortedAnothersData[1, i].iterations
<< ", Sqrt: " << sortedAnothersData[1, 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 < sortedAnothersData.size(); ++i) {
std::cout << "Number: " << sortedAnothersData[2, i].number
<< ", Iterations: " << sortedAnothersData[2, i].iterations
<< ", Sqrt: " << sortedAnothersData[2, i].sqrtResult << "\n";
}

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

std::cout << "\nSome random numbers and their iterations for Carmack's method 2:\n";
for (size_t i = 0; i < 10 && i < sortedCarmack2Data.size(); ++i) {
std::cout << "Number: " << sortedCarmack2Data[i].number
<< ", Iterations: " << sortedCarmack2Data[i].iterations
<< ", Sqrt: " << sortedCarmack2Data[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], resultCarmack, resultCarmack2;
int iterationsNewton, iterationsAnothers[3], iterationsCarmack, iterationsCarmack2;
showIterations(number, resultNewton, resultAnothers, resultCarmack, resultCarmack2, iterationsNewton, iterationsAnothers, iterationsCarmack, iterationsCarmack2);

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

std::vector<double> timesSqrt, timesNewton, timesAnothers[3], timesCarmack, timesCarmack2, timesDivision, timesMultiplication;
std::vector<IterationData> newtonData, anothersData[3], carmackData, carmack2Data;

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, timesCarmack, timesCarmack2, timesDivision, timesMultiplication, newtonData, anothersData, carmackData, carmack2Data);

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], resultCarmack, resultCarmack2 };
calculateAndDisplayErrors(number, roots);

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

// Display the average times
displayAverageTimes(timesSqrt, timesNewton, timesAnothers, timesCarmack, timesCarmack2, 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, carmackData, carmack2Data);

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

return 0;
}


I came from Pascal, and following the recommendations of user Rish, this wasn't working; I don't know what I was thinking. I think in Pascal it could work this way, but I don't remember. In other words, that double indexing in the another's variables... Another thing I did was to name the most likely discoverers, as advised by Adam Hyland.

#include <iomanip>
#include <cmath>
#include <chrono>
#include <limits>
#include <vector>
#include <random>
#include <algorithm>
#include <numeric>
#include <iostream>

// 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;

// Bit manipulation to get a rough estimate of the square root
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); // Using the initial approximation
double prev_f = 0;
iterations = 0;

while (iterations < MAX_ITERATIONS) {
// 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 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)); // Use multiplication for refinement
}

// Ensure the final value meets the precision requirement
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); // Initial approximation using magic number
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); // Get the inverse square root
float sqrt = 1 / invSqrt; // Calculate the square root by taking the reciprocal of the inverse
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); // Get the inverse square root
float sqrt = 1 / invSqrt; // Calculate the square root by taking the reciprocal of the inverse
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& resultCarmack, double& resultCarmack2, int& iterationsNewton, int* iterationsAnothers, int& iterationsCarmack, int& iterationsCarmack2) {
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";

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

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

void measureTime(double* randomNumbers, std::vector<double>& timesSqrt, std::vector<double>& timesNewton, std::vector<double>* timesAnothers, std::vector<double>& timesCarmack, std::vector<double>& timesCarmack2, std::vector<double>& timesDivision, std::vector<double>& timesMultiplication, std::vector<IterationData>& newtonData, std::vector<IterationData>* anothersData, std::vector<IterationData>& carmackData, std::vector<IterationData>& carmack2Data) {
// 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 startCarmack = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsCarmack = 0;
double result = molerwalshMethod(randomNumbers[i], precision, iterationsCarmack, false);
carmackData.push_back({ randomNumbers[i], result, iterationsCarmack });
}
auto endCarmack = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationCarmack = endCarmack - startCarmack;
timesCarmack.push_back(durationCarmack.count());

// Measuring time for Moler-Walsh method 2
auto startCarmack2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < MAX_ITERATIONS; i++) {
int iterationsCarmack2 = 0;
double result = molerwalshMethod2(randomNumbers[i], precision, iterationsCarmack2, false);
carmack2Data.push_back({ randomNumbers[i], result, iterationsCarmack2 });
}
auto endCarmack2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> durationCarmack2 = endCarmack2 - startCarmack2;
timesCarmack2.push_back(durationCarmack2.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>& timesCarmack, const std::vector<double>& timesCarmack2, 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 avgCarmack = std::accumulate(timesCarmack.begin(), timesCarmack.end(), 0.0) / timesCarmack.size();
double avgCarmack2 = std::accumulate(timesCarmack2.begin(), timesCarmack2.end(), 0.0) / timesCarmack2.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): " << avgCarmack << std::endl;
std::cout << "Moler-Walsh method 2 average time (ms): " << avgCarmack2 << 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>& carmackData, const std::vector<IterationData>& carmack2Data) {
std::vector<IterationData> sortedNewtonData = newtonData;
std::vector<IterationData> sortedAnothersData1 = anothersData1;
std::vector<IterationData> sortedAnothersData2 = anothersData2;
std::vector<IterationData> sortedAnothersData3 = anothersData3;
std::vector<IterationData> sortedCarmackData = carmackData;
std::vector<IterationData> sortedCarmack2Data = carmack2Data;

// 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(sortedCarmackData.begin(), sortedCarmackData.end(), [](const IterationData& a, const IterationData& b) {
return a.iterations > b.iterations;
});

std::sort(sortedCarmack2Data.begin(), sortedCarmack2Data.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 < sortedCarmackData.size(); ++i) {
std::cout << "Number: " << sortedCarmackData[i].number
<< ", Iterations: " << sortedCarmackData[i].iterations
<< ", Sqrt: " << sortedCarmackData[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 < sortedCarmack2Data.size(); ++i) {
std::cout << "Number: " << sortedCarmack2Data[i].number
<< ", Iterations: " << sortedCarmack2Data[i].iterations
<< ", Sqrt: " << sortedCarmack2Data[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], resultCarmack, resultCarmack2;
int iterationsNewton, iterationsAnothers[3], iterationsCarmack, iterationsCarmack2;
showIterations(number, resultNewton, resultAnothers, resultCarmack, resultCarmack2, iterationsNewton, iterationsAnothers, iterationsCarmack, iterationsCarmack2);

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

std::vector<double> timesSqrt, timesNewton, timesAnothers[3], timesCarmack, timesCarmack2, timesDivision, timesMultiplication;
std::vector<IterationData> newtonData, anothersData[3], carmackData, carmack2Data;

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, timesCarmack, timesCarmack2, timesDivision, timesMultiplication, newtonData, { anothersData }, carmackData, carmack2Data);

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], resultCarmack, resultCarmack2 };
calculateAndDisplayErrors(number, roots);

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

// Display the average times
displayAverageTimes(timesSqrt, timesNewton, timesAnothers, timesCarmack, timesCarmack2, 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], carmackData, carmack2Data);

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

return 0;
}


Using bit manipulation and basic algebra, I developed a faster method to find the magic number, including for direct square root calculation:

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdint>

const int NUM_TESTS = 9;

union FloatInt {
float value;
uint32_t bits;
};

union DoubleInt {
double value;
uint64_t bits;
};

// Function to calculate the inverse square root
float inverseSqrt(float x) {
FloatInt num;
num.value = x;
num.bits = 0x5f3759df - (num.bits >> 1);

float y = num.value;
//y = y * (1.5f - (x * 0.5f * y * y)); // Newton-Raphson method
return y;
}

// Function to calculate the magic number
uint32_t calculateMagicNumber(const std::vector<float>& numbers, bool inverse) {
std::vector<uint32_t> bitManipulations(NUM_TESTS);
double sum = 0;

// Calculate bitManipulations from known (inverse) square roots
for (size_t i = 0; i < NUM_TESTS; ++i) {
FloatInt fi;
float rootValue;
if (inverse) {
rootValue = 1.0f / sqrt(numbers[i]);
}
else {
rootValue = sqrt(numbers[i]);
}
fi.value = rootValue;
bitManipulations[i] = fi.bits;
}

DoubleInt doubleSum;
doubleSum.value = 0.0;

// Calculate the magic number from bitManipulations
for (size_t i = 0; i < NUM_TESTS; ++i) {
FloatInt fi;
fi.value = numbers[i];
if (inverse) {
fi.value = static_cast<float>(bitManipulations[i] + (fi.bits >> 1));
}
else {
fi.value = static_cast<float>(bitManipulations[i] - (fi.bits >> 1));
}
doubleSum.value += fi.value;
}

// Calculate the average
uint32_t averageMagic = static_cast<uint32_t>(doubleSum.value / NUM_TESTS);
return averageMagic;
}

int main() {
float number;

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

// Vector of known numbers
std::vector<float> knownNumbers = { 0.5f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 10.0f, 100.0f, 1000.0f };

// Calculate the approximate magic number for inverse square root
uint32_t magic = calculateMagicNumber(knownNumbers, true);

// Display the calculated magic number for inverse square root
std::cout << "Calculated magic number for inverse square root: 0x" << std::hex << magic << std::endl;

// Test the inverse square root function with the calculated magic number
FloatInt testNum;
testNum.value = number;
testNum.bits = magic - (testNum.bits >> 1);

// Calculate the approximate magic number for direct square root
uint32_t magic2 = calculateMagicNumber(knownNumbers, false);

// Display the calculated magic number for direct square root
std::cout << "Calculated magic number for direct square root: 0x" << std::hex << magic2 << std::endl;

// Test the square root function with the calculated magic number
FloatInt testNum2;
testNum2.value = number;
testNum2.bits = magic2 + (testNum2.bits >> 1);

// Test the original inverse square root function
float testNum3 = inverseSqrt(number);

std::cout << "Square root of " << number << " with calculated inverse magic number is approximately: " << std::dec << 1 / testNum.value << std::endl;
std::cout << "Square root of " << number << " with calculated magic number is approximately: " << std::dec << testNum2.value << std::endl;
std::cout << "Square root of " << number << " with original (0x5f3759df) magic number is approximately: " << std::dec << 1 / testNum3 << std::endl;
std::cout << "Square root of " << number << " with standard sqrt() is approximately: " << std::dec << sqrt(number) << std::endl;

return 0;
}
`

I believe it was discovered using a similar method, where, given many different inputs, the result was a similar magic number.