Studies on Square Roots

Question

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 <thread>
#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);
unsigned int numThreads = std::thread::hardware_concurrency(); // Determine the number of threads based on the processor
std::vector<unsigned int> threadProgress(numThreads, 0);

// 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
        threadProgress[threadIndex] = magicNumber - start + 1;

        // 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 << "Thread " << threadIndex + 1 << ": ";
            showProgressBar(threadProgress[threadIndex], totalIterations);
            std::cout.flush();
        }
    }
    threadProgress[threadIndex] = totalIterations; // Mark thread as done
}

// Function to display the header
void showHeader() {
    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 << "Processor threads available: " << std::thread::hardware_concurrency() << std::endl;
    std::cout << "Threads used: " << numThreads << std::endl;
    std::cout << "-----------------------------------------\n";
}

// Function to display overall progress
void showProgress(unsigned int totalIterations) {
    while (!done) {
        std::this_thread::sleep_for(std::chrono::seconds(SLEEP_DURATION));

        // 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 << ": ";
            showProgressBar(threadProgress[i], totalIterations);
            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;
    std::vector<std::thread> threads;

    unsigned int totalIterations = rangeSize;

    // Display initial header
    showHeader();

    std::thread progressThread(showProgress, totalIterations);

    for (unsigned int i = 0; i < numThreads; ++i) {
        unsigned int start = startMagicNumber + i * rangeSize;
        unsigned int end = (i == numThreads - 1) ? endMagicNumber : start + rangeSize;
        threads.emplace_back(testMagicNumbersRange, start, end, std::ref(testValues), i);
    }

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

    done = true;
    progressThread.join();

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

Answer 1

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.

Answer 2

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.

Answer 3

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
 * it under the terms of the GNU Lesser General Public License as published by
 * 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.
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * 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
```

Answer 4

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.