# Newton-raphson to find the square root of an arbitrarily sized number with arbitrary precision

I read here that "Furthermore, for a zero of multiplicity 1, the convergence is at least quadratic (see rate of convergence) in a neighbourhood of the zero, which intuitively means that the number of correct digits roughly doubles in every step.". The code which I've written to find the square root using Newton-Raphson has it's time increase quadratically (a.k.a. if I double the desired accuracy, the time taken approximately quadruples): effectively no better than the pen-and-paper square root algorithm.

Here are some graphs I've made.
(The code is immediately below these graphs, if you'd like to look at that first)

# Graph 1

Surprisingly, the number of iterations if your initial seed value is the number itself to calculate the square root of a number with n digits seems to be directly proportional to n.

(The number's always been '1', '10', '100', '1000', ...)

# Graph 2

Here's a graph of the time it takes to find the square root of a number with n digits; it seems to be polynomial in nature.

(Again, the number's always been '1', '10', '100', '1000', ...)

# Graph 3

And finally, a graph of time vs accuracy (I deem this graph to be quadratic in nature).

(I increase the accuracy by 1 every time, and the number I've found the square root of each time is '1000')

I find it interesting how many artifacts this graph has, but then again, these could just be artifacts.

# Code

I've used my own arbitrary precision library for this. To locally reproduce this code, just go to the 'releases' section and download the library's '.h', '.hpp', and '.a' files.

#include "basic_math_operations.hpp"
#include <chrono>
#include <iostream>

std::string f(std::string c, std::string x, size_t _accuracy) {
_accuracy /= 2;
std::string prev;
size_t accuracy = 1;

size_t i = 0;
while (true) {
auto start = std::chrono::high_resolution_clock::now();
i++;
x = basic_math_operations::multiply(
x, basic_math_operations::divide(c, x, accuracy * 2)),
"0.5"); // x = 0.5(x + c/x)
auto end = std::chrono::high_resolution_clock::now();
if (x.find('.') != std::string::npos)
x = x.substr(0, x.find('.') + 2 * accuracy);
if (x.substr(x.find('.') + 1, accuracy) == prev)
break;
prev = x.substr(x.find('.') + 1, accuracy);

if (accuracy * 2 <= _accuracy)
accuracy *= 2;
else
accuracy = _accuracy;
}
return x;
}

int main() {
std::string n;
size_t x;
std::cout << "Enter a number: ";
std::cin >> n;
std::cout << "Enter the number of decimal places: ";
std::cin >> x;
auto start = std::chrono::high_resolution_clock::now();
std::string answer = f(n, n, x);
auto end = std::chrono::high_resolution_clock::now();

std::cout << "sqrt(" << n << ") = " << answer << " (completed in "
<< std::chrono::duration_cast<std::chrono::microseconds>(end -
start)
.count() *
1e-3
<< "ms)\n";
}


How do I improve my implementation of the algorithm?

# Examples

Since @greybeard asked for this, here's an example of the program finding the square root of 10 with an accuracy of 4.

Enter a number: 10
Enter the number of decimal places: 4
x_0 = 10
x_1 = 5.5 (completed in 15μs)
x_2 = 3.659 (completed in 13μs)
x_3 = 3.195 (completed in 22μs)
x_4 = 3.162 (completed in 18μs)
x_5 = 3.162 (completed in 18μs)
sqrt(10) = 3.162 (completed in 0.253ms)


Note that it says completed in 0.253 ms even though each of the individual calculations only took microseconds due to the overhead of having to print each of those statements. x_n is the value of x after the nth loop.

• How does your library compare to, say, the GNU Multiple Precision arithmetic library? What is your purpose in creating another biginteger library, in part in an architecture and toolset dependent language? Commented Oct 29, 2022 at 5:35
• Base 10 big endian terminated ascii. I realise this isn’t the most efficient method, but it’s what I’ve done, and it seems far too tedious to change it now. Commented Oct 29, 2022 at 7:12
• My library does not compare at all to GMP, it’s several times better. This is more of a personal project than an actual competitor to GMP. Commented Oct 29, 2022 at 7:12
• @greybeard just added that example you requested Commented Oct 29, 2022 at 7:29
• Alright yeah, that's a tiny bit faster. Commented Oct 29, 2022 at 8:53

# Avoid mixing high and low-level ways of manipulating numbers

On one hand you have basic_math_operations to do operations on the arbitrary precision numbers, which result in some relatively clean code, but then to check the accuracy, you are using low-level string manipulation. It would be much nicer if you added high-level functions to basic_math_operations to do what you want. Consider being able to write:

x = basic_math_operations::limit_precision(x, 2 * accuracy);
if (basic_math_operations::fraction(x) == prev)
break;


Where limit_precision() would return the value with the precision limited to the given number of digits, and fraction() would return the fractional part.

# Avoid repeating yourself

Avoid repeating the same code or expressions where possible. For example, when checking whether you reached the required precision, you wrote x.substr(x.find('.') + 1, accuracy) twice. In this case, just create a temporary variable:

auto fractional_part = x.substr(x.find('.') + 1, accuracy);

if (fractional_part == prev)
break;

prev = fractional_part;


Reducing duplication will often reduce the possibility of bugs, and by naming the things you are reusing, the code becomes more self-documenting. It is also easier to change the code later; if you want to use a different part of the number to compare to the previous iteration, you now only have to change one line.

# Checking for convergence

The way you check for convergence works for the case of the square root. However, if you have a more complicated function, consider that the Newton-Rhapson method might converge, but oscillate around the limit. If that limit happens to be an integer value, for example 1, then you'll have a series that might look like: 0.9 1.09 0.999 1.0009 1.99999... If you just look at the digits after the decimal point, then they will same to flip between zeroes and nines with each step, seemingly never converging.

The proper way would be to subtract the current value from the previous one, and check whether the absolute difference is less than $$\10^{-\mathtt{accuracy}}\$$.

# Consider using namespace basic_math_operations in your code

Writing basic_math_operations:: every time you need to do an operation on a big number requires a lot of typing, and clutters the code. You can avoid that by using namespace basic_math_operations. Note that you don't have to put that at the top of the source file, you can do this inside f(), so the effect of it is limited to that function only:

std::string f(std::string c, std::string x, size_t _accuracy) {
using namespace basic_math_operations;
...
x = multiply(add(x, divide(c, x, accuracy * 2)), "0.5");
...
}

• So my actual implementation of the algorithm (speed wise for only square roots) is perfect? (Again, only speed wise) Commented Oct 29, 2022 at 15:11
• Perfect in what way? The algorithm is just Newton-Rhapson, nothing to change there. You can probably calculate rather exactly how many steps you need for a given number (as you already have seen, it's proportional to the number of digits), so you can avoid checking at each step if you have reached the desired accuracy. Using strings with decimal numbers is not optimal though, and you can make a better initial seed value as well (consider that the number of digits in the square root is approximately half of that in the original number). Commented Oct 29, 2022 at 18:41