More Square Root

Question

This is a follow up question to my Studies on Square Roots post. I'm posting additional code for solving for square root similar to the other, but with more code to benchmark.

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

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

// 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;
    iterations = 0;

    while (true) {
        nextGuess = 0.5 * (guess + number / guess);
        if (std::fabs(nextGuess - guess) < precision) {
            break;
        }
        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;
        }

        guess = nextGuess;
    }
    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;
}

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

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

    // Refine the estimate using the Newton-Raphson method
    double y = db.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;
    initialSqrtEstimate(number);
    int p = int(number) >> 1;
    double guess = 0;
    double nextGuess = db.value + (number / (1 << p));;
    iterations = 0;
    while (true) {
        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;
        }
        if (std::fabs(nextGuess - guess) < precision) {
            break;
        }
        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) {
    double f = initialSqrtEstimate(number); // Using the initial approximation
    double prev_f = 0;
    iterations = 0;

    while (true) {
        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;
        }

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

        prev_f = f;
        f = 0.5 * (f + number * reciprocal(f)); // Use multiplication for refinement
        
    }
    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) {
        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;
        }
        if (std::fabs(sqrt - prev_sqrt) < precision) {
            break;
        }
        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<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 < 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[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 < NUM_EXECUTIONS; 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 < NUM_EXECUTIONS; 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 < NUM_EXECUTIONS; 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 < NUM_EXECUTIONS; 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());
}

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) {
    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();

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

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;
    std::vector<IterationData> newtonData, anothersData[3], carmackData, carmack2Data;

    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, timesCarmack, timesCarmack2, 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);

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

Sorry for the duplicate code.

Answer 1

use a repo

Sorry for the duplicate code.

Use git for repeated reviews, please. It tracks the edits, revealing what has been changed. A GitHub URL would be welcomed.

Try to keep function definitions in roughly the same order during a refactor, to avoid spurious diffs.

stopping criterion

It's hard to understand why this line appears in your second submission:

        if (std::fabs(nextGuess - guess) < precision) {

The top-voted review of your first submission, from Chux, made essentially one point, about logarithmic sensitivity. This expression mentions the term precision but it's more about absolute error than sig figs. Minimally it warrants a // comment describing why you rejected the review remark in favor of the current approach.

Rather than comparison to the raw delta, I would expect to see the delta normalized by the input arg. We call this the relative error.

meaning of the error

Caller cares about how far away root * root is from the input argument. But what we're measuring here is step size, speed of convergence. Why should caller be concerned with that? Caller is just interested in an accurate estimate, and shouldn't know or care about the underlying method that was used.

test data

It is wonderful that this codebase offers support for assessing accuracy and speed over diverse inputs in a systematic way. But diversity of inputs, currently, is rather limited. A distribution from \$1\$ to \$10^2\$ covers just two decades. Consider going from tiny numbers, much less than unity, up to very large numbers.

extra step

Consider unconditionally doing one more step of the method before returning an estimate. It won't hurt anything, and it takes just a moment.

Ideally a call of the form newtonSqrt(n * n) with integer argument would return exactly n. Arranging for this with a fixed tolerance test is challenging. But the method offers quadratic convergence, so often when that test shows "we're close enough" a single additional step will clear out all the low order bits and we return a simple integer. In messier cases, like newtonSqrt(1 + n * n), it does no harm.

constants

const int NUM_EXECUTIONS = 1000000;

Prefer to spell a million in this way:

1'000'000

Then it's immediately apparent that we're not talking about e.g. ten million.

meaningful identifier

bool showNewtonIteration = true;
bool showAnothersIteration = true;

The Newton-Raphson method is certainly familiar. But mentioning some other anonymous individual is not helpful. Marcus Müller's insightful answer makes it clear that this technique for approximating the derivative is just Newton's method in another guise. Please give anothersMethod() a sensible name.

I find its output slightly less desirable, as it approaches the root from below. For an irrational root like \$\sqrt{10}\$ this makes little difference, but for a perfect square like \$\sqrt{16}\$ we're likely to see lots of trailing 9's. I find lots of 0's followed by low order noise bits to be a little easier to read.

comments with appropriate level of detail

In initialSqrtEstimate(), this comment is accurate but less helpful than it could be.

    // Bit manipulation to get a rough estimate of the square root
    db.bits = (db.bits & 0x800FFFFFFFFFFFFFULL) | (((db.bits & 0x7FF0000000000000ULL) >> 1) + 0x1FF0000000000000ULL);

We're picking out sign, (biased) exponent, and significand. Naming masks wouldn't necessarily make this code clearer. Do consider naming some of the intermediate results to clarify that we're cutting the exponent in half. Or introduce a helper function, which would let your test suite demonstrate how large and small exponents are correctly handled.

(In contrast, the "// Newton-Raphson iteration for reciprocal" comment was lovely and I thank you for it.)

In reciprocal() we see the very magic number 0x7FDDAC7200000000. What I find troubling is neither the code nor the comments describe the range of anticipated x inputs for which this is a sensible initialization value. John Cormack had a range of inputs that were relevant for Quake III rendering, but here you seem to invite all double values. I'm skeptical that constant is ideal for all use cases.

This is immediately followed by another very magic number,

    for (int i = ...; i < 2; i++) {

After FP gave us a "good" initial guess, we go through two more iterations, each enjoying quadratic convergence. But there's a distressing lack of clarity on what accuracy a caller can expect from this, there's no error bounds.

columnar output

The output of displayAverageTimes() is comprised of elapsed times which are directly comparable to one another. Please line them up so the figures are easy to visually scan and compare.

---

Overall, the biggest aspect I found missing in these "accurate approximation routines" was a description of their promised accuracy. One could use mathematical arguments to make accuracy claims. Or one could use this code's extensive measurement infrastructure to make empirical claims about bounds on the observed error.

Answer 2

Illegal type punning

Using union for type punning is legal in C, but undefined behavior in C++. So this

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

is illegal. The legal way to do type punning in C++ is either using memcpy (in most cases the compiler is smart enough to figure out what you're trying to do) or C++20's std::bit_cast.

Pointless comma operator

Not sure what sortedAnothersData[0, i] and sortherAnothersData[1, i] is supposed to be doing. They are the same thing. In fact, compiling with clang 15.0.0 reports this as a warning.

$ g++ -o main main.cc -std=c++17 -Wall -pedantic 
main.cc:260:25: warning: unused variable 'result' [-Wunused-variable]
        volatile double result = std::sqrt(randomNumbers[i]);
                        ^
main.cc:404:55: warning: left operand of comma operator has no effect [-Wunused-value]
        std::cout << "Number: " << sortedAnothersData[0, i].number
                                                      ^
main.cc:405:55: warning: left operand of comma operator has no effect [-Wunused-value]
            << ", Iterations: " << sortedAnothersData[0, i].iterations
                                                      ^
main.cc:411:55: warning: left operand of comma operator has no effect [-Wunused-value]
        std::cout << "Number: " << sortedAnothersData[1, i].number
                                                      ^
main.cc:412:55: warning: left operand of comma operator has no effect [-Wunused-value]
            << ", Iterations: " << sortedAnothersData[1, i].iterations

You can also use the [[maybe_unused]] attribute on result to remove the compiler's warning.

Avoid new and delete

I'm not sure of the reasoning behind using new and delete to allocate the array for the random numbers (especially since you're using std::vector everywhere else). You can easily create a std::vector containing random numbers using std::generate or std::ranges::generate (if using C++20). C++26 will also come with std::ranges::generate_random.

Avoid duplication

I see a lot of duplicated code that can be refactored. For instance, your displaySomeIterations function is doing the same thing multiple times, the only difference being the object you are using. You can refactor it into something like

static void displaySomeIterationsImpl(const std::vector<IterationData>& data, const char* method) {
   std::sort(data.begin(), data.end(), [](const auto& a, const auto& b) {
        return a.iterations > b.iterations;
        });

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

static void displaySomeIterations(...) {
   displaySomeIterationsImpl(newtonsData, "Newton's data");
   ...
}

Also, the && i < data.size() is unnecessary since the data size is NUM_EXECUTIONS which is 1'000'000 and therefore, the loop will always end at 10.