Inspired by all answers, decided to implement in pure C++ several best methods from these answers. As everybody knows C++ is always faster than Python.
To glue C++ and Python I used Cython. It allows to make out of C++ a Python module and then call C++ functions directly from Python functions.
Also as complementary I provided not only Python-adopted code, but pure C++ with tests too.
Here are timings from pure C++ tests:
Test 'GMP', bits 64, time 0.000001 sec
Test 'AndersKaseorg', bits 64, time 0.000003 sec
Test 'Babylonian', bits 64, time 0.000006 sec
Test 'ChordTangent', bits 64, time 0.000018 sec
Test 'GMP', bits 50000, time 0.000118 sec
Test 'AndersKaseorg', bits 50000, time 0.002777 sec
Test 'Babylonian', bits 50000, time 0.003062 sec
Test 'ChordTangent', bits 50000, time 0.009120 sec
and same C++ functions but as adopted Python module have timings:
Bits 50000
math.isqrt: 2.819 ms
gmpy2.isqrt: 0.166 ms
ISqrt_GMP: 0.252 ms
ISqrt_AndersKaseorg: 3.338 ms
ISqrt_Babylonian: 3.756 ms
ISqrt_ChordTangent: 10.564 ms
My Cython-C++ is nice in a sence as a framework for those people who want to write and test his own C++ method from Python directly.
As you noticed in above timings as example I used following methods:
math.isqrt, implementation from standard library.
gmpy2.isqrt, GMPY2 library's implementation.
ISqrt_GMP - same as GMPY2, but using my Cython module, there I use C++ GMP library [
] directly.ISqrt_AndersKaseorg, code taken from answer of @AndersKaseorg.
ISqrt_Babylonian, method taken from Wikipedia article, so-called Babylonian method. My own implementation as I understand it.
ISqrt_ChordTangent, it is my own method that I called Chord-Tangent, because it uses chord and tangent line to iteratively shorten interval of search. This method is described in moderate details in my other article. This method is nice because it searches not only square root, but also K-th root for any K. I drew a small picture showing details of this algorithm.
Regarding compiling C++/Cython code, I used GMP library. You
need to install it first, under Linux it is easy through sudo apt install libgmp-dev.
Under Windows easiest is to install really great program VCPKG, this is software Package Manager, similar to APT in Linux. VCPKG compiles all packages from sources using Visual Studio [don't forget to
install Community version of Visual Studio]. After installing VCPKG you can install GMP by vcpkg install gmp. Also you may install MPIR, this is alternative fork of GMP, you can install it through vcpkg install mpir.
After GMP is installed under Windows please edit my Python code and replace path to include directory and library file. VCPKG at the end of installation should show you path to ZIP file with GMP library, there are .lib and .h files.
You may notice in Python code that I also designed special handy cython_compile[] function that I use to compile any C++ code into Python module. This function is really good as it allows for you to easily plug-in any C++ code into Python, this can be reused many times.
If you have any questions or suggestions, or something doesn't work on your PC, please write in comments.
Below first I show code in Python, afterwards in C++. See Try it online! link above C++ code to run code online on GodBolt servers. Both code snippets I fully runnable from scratch as they are, nothing needs to be edited in them.
def cython_compile[srcs]:
import json, hashlib, os, glob, importlib, sys, shutil, tempfile
srch = hashlib.sha256[json.dumps[srcs, sort_keys = True, ensure_ascii = True].encode['utf-8']].hexdigest[].upper[][:12]
pdir = 'cyimp'
if len[glob.glob[f'{pdir}/cy{srch}*']] == 0:
class ChDir:
def __init__[self, newd]:
self.newd = newd
def __enter__[self]:
self.curd = os.getcwd[]
os.chdir[self.newd]
return self
def __exit__[self, ext, exv, tb]:
os.chdir[self.curd]
os.makedirs[pdir, exist_ok = True]
with tempfile.TemporaryDirectory[dir = pdir] as td, ChDir[str[td]] as chd:
os.makedirs[pdir, exist_ok = True]
for k, v in srcs.items[]:
with open[f'cys{srch}_{k}', 'wb'] as f:
f.write[v.replace['{srch}', srch].encode['utf-8']]
import numpy as np
from setuptools import setup, Extension
from Cython.Build import cythonize
sys.argv += ['build_ext', '--inplace']
setup[
ext_modules = cythonize[
Extension[
f'{pdir}.cy{srch}', [f'cys{srch}_{k}' for k in filter[lambda e: e[e.rfind['.'] + 1:] in ['pyx', 'c', 'cpp'], srcs.keys[]]],
depends = [f'cys{srch}_{k}' for k in filter[lambda e: e[e.rfind['.'] + 1:] not in ['pyx', 'c', 'cpp'], srcs.keys[]]],
extra_compile_args = ['/O2', '/std:c++latest',
'/ID:/dev/_3party/vcpkg_bin/gmp/include/',
],
],
compiler_directives = {'language_level': 3, 'embedsignature': True},
annotate = True,
],
include_dirs = [np.get_include[]],
]
del sys.argv[-2:]
for f in glob.glob[f'{pdir}/cy{srch}*']:
shutil.copy[f, f'./../']
print['Cython module:', f'cy{srch}']
return importlib.import_module[f'{pdir}.cy{srch}']
def cython_import[]:
srcs = {
'lib.h': """
#include
#include
#include
#include
#include
#include
#include
#include
#include
#pragma comment[lib, "D:/dev/_3party/vcpkg_bin/gmp/lib/gmp.lib"]
#define ASSERT_MSG[cond, msg] { if [![cond]] throw std::runtime_error["Assertion [" #cond "] failed at line " + std::to_string[__LINE__] + "! Msg '" + std::string[msg] + "'."]; }
#define ASSERT[cond] ASSERT_MSG[cond, ""]
#define LN { std::cout = [1 >= 8;
}
while [n] {
++cnt;
n >>= 1;
}
return cnt;
}
}
template
T ISqrt_Babylonian[T const & y] {
// //en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
if [y >= 1;
}
if [b < a]
std::swap[a, b];
if [b - a > limit]
continue;
++b;
for [size_t i = 0; a y] {
if [i == 0]
break;
else
return a - 1;
}
ASSERT[false];
}
}
template
T ISqrt_AndersKaseorg[T const & n] {
// //stackoverflow.com/a/53983683/941531
if [n > 0] {
T y = 0, x = T[1] > 1];
while [true] {
y = [x + n / x] >> 1;
if [y >= x]
return x;
x = y;
}
} else if [n == 0]
return 0;
else
ASSERT_MSG[false, "square root not defined for negative numbers"];
}
template
T ISqrt_GMP[T const & y] {
// //gmplib.org/manual/Integer-Roots
mpz_class r, n;
bool constexpr is_mpz = std::is_same_v;
if constexpr[is_mpz]
n = y;
else {
static_assert[sizeof[T] > 32];
n 32].get_mpz_t[]]] n] {
x_end = x_mid; y_end = y_mid;
} else {
x_begin = x_mid; y_begin = y_mid;
}
}
// [y_end - y_begin] / [x_end - x_begin] = [n - y_begin] / [x_n - x_begin] ->
x_n = x_begin + [n - y_begin] * [x_end - x_begin] / [y_end - y_begin];
y_n = KthPow[x_n];
tangent_x = x_n + [n - y_n] / KthPowDer[x_n] + 1;
chord_x = x_n + [n - y_n] * [x_end - x_n] / [y_end - y_n];
//ASSERT[chord_x