小道具：開平法（二進数）による整数平方根

開平法（二進数）による平方根【C++版】を見つつ、正の整数専用で Python 版を作ってみる。

def lsqrt(x):
    if x < 0:
        raise ArithmeticError
    l, r, s = x.bit_length() >> 1, 0, 0
    for b in range(l, -1, -1):
        r <<= 1
        t = 1 << b
        u = s + t
        v = u << b
        if x < v:
            continue
        x -= v
        r |= 1
        s = u + t
    return (r, x)

テスト

>>> [lsqrt(x) for x in range(1,101)]
[(1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2), (2, 3), (2, 4), (3, 0), (3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6), (4, 7), (4, 8), (5, 0), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (6, 0), (6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6), (6, 7), (6, 8), (6, 9), (6, 10), (6, 11), (6, 12), (7, 0), (7, 1), (7, 2), (7, 3), (7, 4), (7, 5), (7, 6), (7, 7), (7, 8), (7, 9), (7, 10), (7, 11), (7, 12), (7, 13), (7, 14), (8, 0), (8, 1), (8, 2), (8, 3), (8, 4), (8, 5), (8, 6), (8, 7), (8, 8), (8, 9), (8, 10), (8, 11), (8, 12), (8, 13), (8, 14), (8, 15), (8, 16), (9, 0), (9, 1), (9, 2), (9, 3), (9, 4), (9, 5), (9, 6), (9, 7), (9, 8), (9, 9), (9, 10), (9, 11), (9, 12), (9, 13), (9, 14), (9, 15), (9, 16), (9, 17), (9, 18), (10, 0)]
>>> [r * r + s for (r, s) in [lsqrt(x) for x in range(1,101)]]
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]

C++ 版と違って多倍長整数の利用が容易なので、高精度の固定小数点として使ってみる

固定小数点の10進数文字列化

fixdec = (lambda x, b, p=0: (lambda x, b, p, f: f(x >> b, x & ((1 << b) - 1), b, p))
          (x, b, (b, int(b * 0.301), p)[min(2, p + 1)],
           lambda h, l, b, p: (('%d.' % h) + ('%d' % (((l | (1 << b)) * (10**p)) >> b))[1:])))

テスト

>>> fixdec(lsqrt(2 << 256)[0], 128)
'1.41421356237309504880168872420969807856'
>>> fixdec(lsqrt(2 << 256)[0], 128, -1)
'1.41421356237309504880168872420969807856899607931908261748958951509877704601116619082283240582675887253572000190615653991699218750'