# Fastest square root of an arbitrary size

We do seem to have a fastest square root challenge, but it's very restrictive. In this challenge, your program (or function) will be given an arbitrarily sized nonnegative integer, which is the square of an integer.

Input

You may take input in any normally allowed way, including (but not limited to) an arbitrary sized integer type, string in any base, or list of numbers representing digits.

Output

The allowed output methods are similar to the input ones, they just need to be in a way which can also represent a number of arbitrary size.

Requirements

• Your program should theoretically work correctly given any number (which is the square of an integer), assuming the language is capable of this
• All submissions should be in either:
• A language with an online interpreter that can give (or be trivially modified to give) accurate times, and doesn't time out with your solution
• A language with an interpreter/compiler available from the apt package manager, and is less than 1 GiB, in which case it will be run on my Chromebook's Linux VM
• If it doesn't take under 20 minutes, I'll just give a good estimate based on the first test cases and time complexity
• Built in square root functions are allowed (but boring :p)

Scoring

Your program's score will be the amount of time it takes to find the square root of:

• All 4096 perfect squares up to (not including) $$\2^{24}\$$
• A random (but same for all programs) 2048 perfect squares up to $$\2^{48}\$$
• A random (but same for all programs) 1024 perfect squares up to $$\2^{144}\$$
• A random (but same for all programs) 512 perfect squares up to $$\2^{1440}\$$
• A random (but same for all programs) 256 perfect squares up to $$\2^{5760}\$$
• A random (but same for all programs) 128 perfect squares up to $$\2^{14400}\$$
• A random (but same for all programs) 64 perfect squares up to $$\2^{86400}\$$

UPDATE: The numbers have been made bigger. This is a total of 8128 numbers as inputs, specifically these. Remember that fitting the test cases specifically is a standard loophole. The only part that counts toward the time is the math. I/O and setup can be ignored.

• Just to note, Dennis has said in the past that TIO is not a good timing method for fastest-code challenges, which presumably would be applicable to other online interpreters. It might be best for you to test all submissions on your machine, requiring people to include testing instructions. Nov 11, 2020 at 16:26
• I think you need more or bigger test cases. With Python 3.8's math.isqrt, the actual computation seems to take little time relative to that for input/output.
– xnor
Nov 11, 2020 at 20:34
• This is a very important question: What parts are required to be in the time: the actual maths, the outputting (print/console.log/System.out.println/prinf/etc.), the setup (reading from a file/initializing an array? I would assume only the maths AND the outputting should be included in the timing, but we need to know for sure to stay consistent for everyone. Nov 11, 2020 at 23:11
• I also think this needs larger testcases. Currently I/O takes more time than the actual square root. Test cases on the order of $2^{10^4}$ to $2^{10^5}$ should do. Nov 12, 2020 at 0:45
• Can I suggest a more scalable benchmarking strategy? Describe an infinite sequence of squares of exponentially increasing length, using some simple deterministic algorithm with a seed provided as input. Your score is the number of these square roots you can take in 10 seconds. Nov 12, 2020 at 20:44

# C++, 140 µs

This is a trivial answer that uses GMP, but I’m submitting it anyway to drive home the point that the test cases are way too small. Scores in the range of literally nanoseconds per input cannot meaningfully be compared. We need larger test cases for this to be an interesting challenge.

#include <chrono>
#include <gmpxx.h>
#include <iostream>
#include <vector>

int main() {
std::vector<mpz_class> numbers;
mpz_class number;
while (std::cin >> number)
numbers.push_back(number);
auto t0 = std::chrono::high_resolution_clock::now();
for (auto &number : numbers)
number = sqrt(number);
auto t1 = std::chrono::high_resolution_clock::now();
for (const auto &number : numbers)
std::cout << number << '\n';
std::cout
<< std::chrono::duration_cast<std::chrono::nanoseconds>(t1 - t0).count()
<< "ns\n";
}


Try it online!

• I have added larger test cases. Sorry about that, I guess I shouldn't have underestimated C/C++ :p Nov 12, 2020 at 14:42
• @RedwolfPrograms Words of wisdom Nov 13, 2020 at 16:19

# JavaScript (Node), Unknown time due to TIO not liking 2.6 million characters as the input.

If you want to see the previous tests, then Try it online!. Currently, we can't use TIO with all of the inputs.

function iroot(base) {
let s = base + 1n;
let u = base;
if (!base) { console.log(0n); return; }
while (u < s) {
s = u;
u = (u + base / u) / 2n;
}
return s;
}

const inputs = [/* the gist goes here*/]

console.time("sqrt");

inputs.map(iroot);

console.timeEnd("sqrt");

for (let i = 0; i < iroot.length; i++) {
console.log(inputs[i]);
}


This is based off of a StackOverflow answer I did, which was based off of Newton's formula. This function is better than the built in Math.sqrt() or Math.pow() or ** since this gives precision even to the last digits instead of 1.457746e57 or something like that. Also with BigInt, all numbers are integers, so you can't use ** a decimal/fraction.

This code expects a long array to be prefilled with the numbers as bigints: https://gist.github.com/Samathingamajig/23b8107f299649d6bd3399b04f65c257 , which can be compressed into one line through this gist: https://gist.github.com/Samathingamajig/88f64e58ba21bd8d392be23bbb85dd9d

A Gist with all the data as raw JSON can be found here: https://gist.github.com/Samathingamajig/c15fd1f15a7bca6181a2e792c469503d

• why the second argument as the challenge is for square root? just using 2n instead and removing k1 gives 111ms, also it's seems timer gives user time not sys time, which is less Nov 11, 2020 at 18:50
• Your TIO link suffers from rounding errors in its inputs array; you need to suffix n on all of those numbers otherwise they'll get converted to float and back to BigInt.
– Neil
Nov 11, 2020 at 22:37
• @NahuelFouilleul fixed, no improvements in time were given, it is still randomly between 100-200ms. What do you mean by "also it's seems timer gives user time not sys time"? Nov 12, 2020 at 0:52
• @Neil I fixed that Nov 12, 2020 at 0:52
• @Samathingamajig about the different times i'm refering the debug section given by tio. for more details stackoverflow.com/questions/556405/… finally it's user+sys time which should be taken Nov 12, 2020 at 13:55

# Rust 450 μs

use std::time::SystemTime;
use num_bigint::BigUint;

fn main() {
// Parse the numbers
let mut numbers: Vec<_> = INPUT.iter().map(|i| {
i.parse::<BigUint>().unwrap()
}).collect();

let start = SystemTime::now();

// Do the calculation
numbers.iter_mut().for_each(|s| *s = s.sqrt());

// Print the time.
println!("{}us", start.elapsed().unwrap().as_micros());
}


I agree that the input is too small. If performance was an issue, we would not use a library, but create a custom integer algorithm for monotonic series of integers, do them on SIMD hardware and multithreaded on 64 cores.

• Rust playground: play.rust-lang.org/… Nov 12, 2020 at 11:42
• The OP's Chromebook won't have 64 cores, but yes good point about threading the computation since the test is guaranteed to be multiple separate sqrt calculations. No need to find parallelism within one single sqrt. And yes maybe even interleave limbs of 2 or 4 integers to allow vertical SIMD, perhaps using double-precision FMA instructions to do 52 integer bits per 64-bit chunk, like @Mysticial's code does for prime95 / y-cruncher, see his answer on Can I use the AVX FMA units to do bit-exact 52 bit integer multiplications?. Nov 12, 2020 at 12:11
• Although Goldmont Atom still only has SSE4.x, not AVX/FMA, so SIMD is much harder to use for extended-precision, beyond how GMP already uses it (e.g. for right shifts, not multiplies or divides IIRC) Nov 12, 2020 at 12:13
• The impl of BigUint::sqrt starts well with a f64 approximation, but then uses a refinement algorithm with a divide which significantly increases the cost. github.com/rust-num/num-bigint/blob/master/src/biguint.rs#L1735 I'd need to check to see what gmp does. Nov 12, 2020 at 12:37
• I have added larger test cases, might take somewhere in the range of milliseconds now :p Nov 12, 2020 at 14:44

# JavaScript (Node.js), < 60 ms on TIO

function iroot(input) {
let bit = 1n << BigInt(input.toString(2).length - 1 >> 1);
let root = 0n;
while (bit) {
let diff = (root + root + bit) * bit;
if (input >= diff) {
input -= diff;
root += bit;
}
bit /= 2n;
}
return root;
}

console.time("sqrt");
let roots = inputs.map(iroot);
console.timeEnd("sqrt");

for (let root of roots) {
console.log(root);
}


Try it online! Can't remember the name of this method sorry, but it's very handy for square rooting floating-point numbers by shifting as the bit is always 1 thus reducing the complexity even further.

• If you move the logging into the timing section, like here (copy the link and paste to get the tio) then you get the same timing range as mine. We are awaiting the OP's ruling on if timing should count for I/O, and I will change mine if the decision is to disregard I/O Nov 12, 2020 at 0:51
• @Samathingamajig And I will change mine if the decision is to include I/O, but you'll need to ping me since I'm not watching the question.
– Neil
Nov 12, 2020 at 1:01
• @Neil I/O no longer counts toward the time, sorry for the rule change Nov 12, 2020 at 4:41
• I have added larger test cases (sorry for the change halfway :/) Nov 12, 2020 at 14:42
• @RedwolfPrograms Well, I tried the largest test case, and it took 25s on TIO...
– Neil
Nov 14, 2020 at 2:36

# Ruby, ~16-28 ms on TIO

def newtonRoot base
s = base + 1
u = base
if base == 0
return 0
end
while u < s do
s = u
u = (u + base / u) / 2
end
return s
end

t = Time.now
outputs = inputs.map{ |i| newtonRoot i }
puts "> time elasped: #{ (Time.now - t) * 1000 } ms.\n\n"

puts outputs.map{ |v| v * v } == inputs

# un-comment below line to see outputs!
# outputs.map{ |i| print "#{i}, " }


The former code is a modification of latter code, which is specific to the square root of a "perfect square", and a bit faster!

Try it online!

# Ruby, ~22-36 ms on TIO

def newtonRoot (base, root)
s = base + 1
k1 = root - 1
u = base
if (base == 0)
print "0, "
return
end
while (u < s) do
s = u
u = ((u * k1) + base / (u ** k1)) / root
end
print "#{s}, "
end

t = Time.now

inputs.each{|i| newtonRoot(i, 2)}

puts "\n\n> time elasped: #{ (Time.now - t) * 1000 } ms."


A remix of Newton's formula in Ruby. I posted this as I think Ruby is very fast.

• Can you please provide the TIO link for your second solution? Nov 12, 2020 at 3:42
• Maybe you could move your new solution to the top, remove the old TIO link and use the characters to paste the new TIO link? Usually when people improve their solution, they unshift their new solution to the top. Nov 12, 2020 at 5:29
• @Samathingamajig - Okay, done. Nov 12, 2020 at 5:40
• I have now added larger test cases Nov 12, 2020 at 14:41

# Python 3, unknown time

Finishes the first 7979 cases on my machine after 19.2 minutes. (Using head -n 7979 input.txt | python3 sqrt.py)

from time import time

numbers = []
try:
while True:
numbers += [int(input())]
except:
pass

t1 = time()

for i in range(len(numbers)):
base = numbers[i]
u, s = base, base + 1
if u:
while u < s:
u, s = (u + base // u) // 2, u
numbers[i] = u

t2 = time()
print(t2 - t1, "seconds")


Uses the same method as Samathingamajig's answer.

• @Samathingamajig, no, it's faster because it doesn't multiply by 1 and doesn't raise to the power 1, as in my other comment, also the most significant time is the sys time. Nov 11, 2020 at 21:01
• Wdym by "sys time" @NahuelFouilleul Nov 11, 2020 at 21:16
• Has anyone tested this to make sure that there's no problems with losing precision at the high numbers? Nov 11, 2020 at 23:06
• @Samathingamajig I don't think it would lose precision anywhere in the code; it's all arbitrary-precision integers from start to finish. Nov 12, 2020 at 3:38
• Update! OP had specified the I/O as optional in the time count. Nov 12, 2020 at 5:58

# Rust Babylonian + New data 33 ms

This uses the new, larger data set and a novel recursive version of the Babylonian algorithm.

The algorithm refines the guess for the babylonian by calculating the square root of the top half of the number, which is much less expensive.

use num_bigint::BigUint;
use std::time::SystemTime;
use num_traits::{FromPrimitive, ToPrimitive};
use rayon::prelude::*;

fn recursive_babylonian_sqrt(a: &BigUint) -> BigUint {
if a.bits() < 53 {
//println!("a={} exact={}", a, BigUint::from_f64(a.to_f64().unwrap().sqrt()).unwrap());
return BigUint::from_f64(a.to_f64().unwrap().sqrt()).unwrap();
}

let mut guess = recursive_babylonian_sqrt(&(a >> a.bits()/4*2)) << (a.bits()/4);
//println!("a={} sqrt={} guess={}", a, a.sqrt(), guess);

//let mut nits = 0;
loop {
let mut new = a / &guess;
new += &guess;
new >>= 1;
//nits += 1;
if new == guess {
break;
}
guess = new;
}
//println!("{} nits", nits);
guess
}

pub fn main() {
.unwrap();

let numbers : Vec<_> = text[0..text.len()-1]
.split(|c| *c == b'\n')
.map(|s| std::str::from_utf8(s).unwrap().parse::<BigUint>().unwrap() )
.collect();

let start = SystemTime::now();

let result : Vec<_> = numbers
.par_iter()
.map(|s| recursive_babylonian_sqrt(s))
.collect();

// Print the time.
println!("{}us", start.elapsed().unwrap().as_micros());

numbers
.iter()
.zip(result.iter())
.for_each(|(n, r)| assert_eq!(&n.sqrt(), r));
}
$$$$
`