# 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. – caird coinheringaahing Nov 11 '20 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 '20 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. – Samathingamajig Nov 11 '20 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. – Sisyphus Nov 12 '20 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. – Anders Kaseorg Nov 12 '20 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 – Redwolf Programs Nov 12 '20 at 14:42
• @RedwolfPrograms Words of wisdom – Alexandre Cassagne Nov 13 '20 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 – Nahuel Fouilleul Nov 11 '20 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 '20 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"? – Samathingamajig Nov 12 '20 at 0:52
• @Neil I fixed that – Samathingamajig Nov 12 '20 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 – Nahuel Fouilleul Nov 12 '20 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/… – Andy Thomason Nov 12 '20 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?. – Peter Cordes Nov 12 '20 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) – Peter Cordes Nov 12 '20 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. – Andy Thomason Nov 12 '20 at 12:37
• I have added larger test cases, might take somewhere in the range of milliseconds now :p – Redwolf Programs Nov 12 '20 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 – Samathingamajig Nov 12 '20 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 '20 at 1:01
• @Neil I/O no longer counts toward the time, sorry for the rule change – Redwolf Programs Nov 12 '20 at 4:41
• I have added larger test cases (sorry for the change halfway :/) – Redwolf Programs Nov 12 '20 at 14:42
• @RedwolfPrograms Well, I tried the largest test case, and it took 25s on TIO... – Neil Nov 14 '20 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? – Samathingamajig Nov 12 '20 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. – Samathingamajig Nov 12 '20 at 5:29
• @Samathingamajig - Okay, done. – vrintle Nov 12 '20 at 5:40
• I have now added larger test cases – Redwolf Programs Nov 12 '20 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. – Nahuel Fouilleul Nov 11 '20 at 21:01
• Wdym by "sys time" @NahuelFouilleul – Samathingamajig Nov 11 '20 at 21:16
• Has anyone tested this to make sure that there's no problems with losing precision at the high numbers? – Samathingamajig Nov 11 '20 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. – Bubbler Nov 12 '20 at 3:38
• Update! OP had specified the I/O as optional in the time count. – vrintle Nov 12 '20 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));
}
$$$$
`