Fastest square root of an arbitrary size

Question

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.

Answer 1

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!

Answer 2

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

Answer 3

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.

Answer 4

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.

Answer 5

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.

Answer 6

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.

Answer 7

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() {
    let text = std::fs::read("numbers-updated.txt")
        .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));
}
```