'Return the square root of self.
Uses a converging algorithm (Xn+1 = 0.5*(Xn + self / Xn))
Should quadratically approach the right answer.'

def sqrt(self, context=None):

    if self._is_special:
        ans = self._check_nans(context=context)
        if ans:
            return ans
        if (self._isinfinity() and (self._sign == 0)):
            return Decimal(self)
    if (not self):
        exp = (self._exp // 2)
        if (self._sign == 1):
            return Decimal((1, (0,), exp))
        else:
            return Decimal((0, (0,), exp))
    if (context is None):
        context = getcontext()
    if (self._sign == 1):
        return context._raise_error(InvalidOperation, 'sqrt(-x),    x    >    0')
    tmp = Decimal(self)
    expadd = (tmp._exp // 2)
    if (tmp._exp & 1):
        tmp._int += (0,)
        tmp._exp = 0
    else:
        tmp._exp = 0
    context = context._shallow_copy()
    flags = context._ignore_all_flags()
    firstprec = context.prec
    context.prec = 3
    if ((tmp.adjusted() & 1) == 0):
        ans = Decimal((0, (8, 1, 9), (tmp.adjusted() - 2)))
        ans = ans.__add__(tmp.__mul__(Decimal((0, (2, 5, 9), (-2))), context=context), context=context)
        ans._exp -= (1 + (tmp.adjusted() // 2))
    else:
        ans = Decimal((0, (2, 5, 9), ((tmp._exp + len(tmp._int)) - 3)))
        ans = ans.__add__(tmp.__mul__(Decimal((0, (8, 1, 9), (-3))), context=context), context=context)
        ans._exp -= (1 + (tmp.adjusted() // 2))
    (Emax, Emin) = (context.Emax, context.Emin)
    (context.Emax, context.Emin) = (DefaultContext.Emax, DefaultContext.Emin)
    half = Decimal('0.5')
    maxp = (firstprec + 2)
    rounding = context._set_rounding(ROUND_HALF_EVEN)
    while 1:
        context.prec = min(((2 * context.prec) - 2), maxp)
        ans = half.__mul__(ans.__add__(tmp.__div__(ans, context=context), context=context), context=context)
        if (context.prec == maxp):
            break
    context.prec = firstprec
    prevexp = ans.adjusted()
    ans = ans._round(context=context)
    context.prec = (firstprec + 1)
    if (prevexp != ans.adjusted()):
        ans._int += (0,)
        ans._exp -= 1
    lower = ans.__sub__(Decimal((0, (5,), (ans._exp - 1))), context=context)
    context._set_rounding(ROUND_UP)
    if (lower.__mul__(lower, context=context) > tmp):
        ans = ans.__sub__(Decimal((0, (1,), ans._exp)), context=context)
    else:
        upper = ans.__add__(Decimal((0, (5,), (ans._exp - 1))), context=context)
        context._set_rounding(ROUND_DOWN)
        if (upper.__mul__(upper, context=context) < tmp):
            ans = ans.__add__(Decimal((0, (1,), ans._exp)), context=context)
    ans._exp += expadd
    context.prec = firstprec
    context.rounding = rounding
    ans = ans._fix(context)
    rounding = context._set_rounding_decision(NEVER_ROUND)
    if (not (ans.__mul__(ans, context=context) == self)):
        context._regard_flags(flags)
        context._raise_error(Rounded)
        context._raise_error(Inexact)
    else:
        exp = (self._exp // 2)
        context.prec += (ans._exp - exp)
        ans = ans._rescale(exp, context=context)
        context.prec = firstprec
        context._regard_flags(flags)
    (context.Emax, context.Emin) = (Emax, Emin)
    return ans._fix(context)
