Calculate square roots

Suppose you have an equation like this:

$x^2 = a, \;\;\;\; a \in \mathbb{R}_{\geq 0}$

Your task is to compute the solution $x \in \mathbb{R}_{\geq 0}$.

How do you solve this algorithmically?

Input and output

Write a program that takes two parameters $a, n \in \mathbb{N}_{\geq 1}$:

$a$: The number you should take the square root of. $a$ is not bigger than 65535. To make this a little bit simpler, you can also assume that $a$ is an integer.
$n$: Number of iterations you may use

Your program should output exactly one positive floating point number. The decimal separator is a point. At the end of the number should be a newline character \n.

Reference code

// Thanks to http://stackoverflow.com/a/15363123/562769
#include <stdio.h>
#include <gmp.h>

int main(int argc, char *argv[]) {
    mpf_t res, a;
    mpf_set_default_prec(1000000); // Increase this number.
    mpf_init(res);
    mpf_init(a);
    mpf_set_str(a, "2", 10);
    mpf_sqrt (res, a);
    gmp_printf("%.1000Ff\n\n", res); // Increase this number.
    return 0;
}

You need GMP (libgmp-dev) to compile this.

Compile it like this:

gcc sqrt-reference.c  -lgmp -lm -O0 -g3 -o reference.out

This is the script I use to get the number of correct digits:

#!/usr/bin/env python
# -*- coding: utf-8 -*-


def getScore(program, a, n):
    import os

    os.system("./reference.out " + str(a) + " > reference.txt")
    os.system("./" + program + " " + str(a) + " " + str(n) + " > result.txt")

    f = open("reference.txt", "r")
    reference = f.read()
    f.close()

    f = open("result.txt", "r")
    result = f.read()
    f.close()

    points = 0
    areEqual = True
    while reference[points] != "\n" and result[points] != "\n":
        if reference[points] == result[points]:
            points += 1
        else:
            break
    if points >= 2:
        points -= 1  # decimal point
    return points


if __name__ == "__main__":
    from argparse import ArgumentParser

    parser = ArgumentParser()

    # Add more options if you like
    parser.add_argument(
        "-p",
        "--program",
        dest="program",
        help="your program",
        metavar="FILE",
        required=True,
    )
    parser.add_argument(
        "-a", metavar="A", type=int, required=True, help="calculate squre root of a"
    )
    parser.add_argument(
        "-n", metavar="N", type=int, required=True, help="maximum n iterations"
    )

    args = parser.parse_args()

print(
    "Points for a=%i and n=%i: %i"
    % (args.a, args.n, getScore(args.program, args.a, args.n))
)

Newton's method

How should we choose the initial value? I thought $\frac{a}{2}$ could be ok. In a good implementation you'll probably do this with a lookup table.

With long double:

#include <iostream>
#include <cmath>

using namespace std;

long double newton(int a, int n) {
    long double x = ((long double)a)/2;
    for (int i=0; i<n; i++) {
        x = x - (x*x - a)/(2*x);
    }
    return x;
}

int main(int argc, char *argv[]) {
    if (argc != 3) {
        cout << "Please enter exactly two arguments." << endl;
        return 1;
    }
    int a = atoi(argv[1]);
    int n = atoi(argv[2]);
    printf("%.80Lf\n", newton(a, n));
    return 0;
}

I failed to convert this to a version that uses GMP :-/

But it converges quite fast:

Newton's method for calculating square roots

Exponential identity

According to Wikipedia (Source: Methods of computing square roots) many calculators use the following identity:

$\displaystyle \sqrt{a} = e^{\frac{1}{2}\ln a}$

But to calculate this, you need to be able to calculate $e^x$ for $x \in \mathbb{R}_{\geq 0}$ and $\ln(a)$ for $a \in \mathbb{R}_{\geq 0}$.

I've used the definition of $e^x$ to calculate $e^x$: $\displaystyle e^x := \sum_{i = 0}^\infty \frac{x^i}{i!}$

and:

$\displaystyle \ln(1+x) =- \sum_{k=0}^\infty \frac{(-x)^{k+1}}{k+1}$

So I gave it a try:

#include <iostream>
#include <cmath>

using namespace std;

long double ln(int S, int n) {
    long double tmp = S - 1;
    long double result = tmp;
    long double sign = 1.0;
    for (int i=2; i<n/2; i++) {
        tmp *= (S-1);
        sign *= -1;
        result += sign*tmp/i;
    }
    return result;
}

long double e(long double x, int n) {
    long double numerator = 1;
    long double denominator = 1;
    long double result = 1;
    for (int i=1; i<n/2; i++) {
        numerator *= x;
        denominator *= i;
        result += numerator/denominator;
    }
    return result;
}

long double sqrt(int a, int n) {
    return e(ln(a, n)*0.5, n);
}

int main(int argc, char *argv[]) {
    if (argc != 3) {
        cout << "Please enter exactly two arguments." << endl;
        return 1;
    }
    int a = atoi(argv[1]);
    int n = atoi(argv[2]);
    printf("%.80Lf\n", sqrt(a, n));
    return 0;
}

This converges VERY slowly: For $a = 2$

$n=1$: 1 digit correct
$n=10$: 1 digit correct
$n=100$: 2 digits correct
$n=1000$: 4 digit correct
$n=10,000$: 5 digit correct
$n=100,000$: 5 digit correct
$n=1,000,000$: 6 digit correct

Ok, lets take a better Taylor series for calculating $e^x$: $\displaystyle e^x = \sum_{k=0}^\infty \frac{x^{-1+2 k} (2 \cdot k+x)}{(2\cdot k)!}$

This didn't change anything! I'm very surprised ... as I've calculated $\pi$ for another article, changing the series to something similar improved the speed of convergence drastically.

Input and output

Reference code

Newton's method

Exponential identity

See also