package root.sqrt;

/**
 * @author Victor Kuliamin
 */
public class AdvSqrt extends Sqrt {
    private static double eps = 2.25e-16;
    private static long dgmask = 0x7FF0000000000000L;
    private static long mtmask = 0x000FFFFFFFFFFFFFL;
    private static int dgshift = 52;
    private static int sqdgadd = 0x1FF;
    private static double dnrbnd = Double.longBitsToDouble(0x0010000000000000L);
    private static long odddeg = 0x3FF0000000000000L;
    private static long evndeg = 0x3FE0000000000000L;

    public double sqrt(double x) {
        if (Double.isNaN(x) || x < 0) return Double.NaN;
        if (x == -0) return -0.0;
        else if (x == 0 || x == 1 || Double.isInfinite(x)) return x;
        else {
            boolean dnr = false;
            if (x < dnrbnd) {
                x = x * Math.pow(2, dgshift);
                dnr = true;
            }

            long b = Double.doubleToLongBits(x);
            int d = (int) ((b & dgmask) >> dgshift);
            double res, tmp;
            int i = 0;

            if ((d & 1) != 0)
                b = (b & mtmask) | odddeg;
            else
                b = (b & mtmask) | evndeg;

            x = Double.longBitsToDouble(b);
            res = x;

            while (Math.abs(x - res * res) / x > eps && i < 6) {
                i++;
                tmp = res;
                res = (tmp + x / tmp) / 2;
            }
            res = (res + x / res) / 2;
            b = Double.doubleToLongBits(res);
            d = (int) Math.ceil((double) d / 2) + sqdgadd;
            if (dnr) d -= dgshift / 2;

            b = (b & mtmask) | ((long) d << dgshift);
            res = Double.longBitsToDouble(b);
            return res;
        }
    }

}