# @file bignum.002.py
# @ingroup experimental
# Experimental big-number representation.
# (sign magnitude version)
# @date 11/22/2024
 
import itertools
import operator
import random
 
class Num:
    base = 2**8
 
    def __init__(self, init):
        if isinstance(init, tuple):
            self.m = init[0]
            self.s = init[1]
        else:
            self.m = iton(abs(init), Num.base)
            self.s = init < 0
 
    def __neg__(self):
        return Num(nneg(self.m, self.s))
 
    def __add__(self, other):
        return Num(nadd(self.m, self.s, other.m, other.s, Num.base))
 
    def __sub__(self, other):
        return Num(nsub(self.m, self.s, other.m, other.s, Num.base))
 
    def __mul__(self, other):
        return Num(nmul(self.m, self.s, other.m, other.s, Num.base))
 
    def __floordiv__(self, other):
        return Num(ndiv(self.m, self.s, other.m, other.s, Num.base)[0])
 
    def __mod__(self, other):
        return Num(ndiv(self.m, self.s, other.m, other.s, Num.base)[1])
 
    def __str__(self):
        return '-' * self.s + ntos(self.m, Num.base)
 
    def __repr__(self):
        return str((self.m, self.s))
 
# Utility.
 
def _iszero(a):
    return len(a) == 1 and a[0] == 0
 
def _fixsign(a, s):
    return a, s and not _iszero(a)
 
def _zeroextend(a, n):
    return a[:] + [0] * (n - len(a))
 
def _zeroreduce(a):
    n = len(a)
    while n != 1 and a[n-1] == 0:
        n -= 1
    return a[:n]
 
def _inplace_mul1(r, x, c, base):
    for i in range(len(r)):
        c, r[i] = divmod(r[i] * x + c, base)
    return c, r
 
# Convert.
 
def iton(x, base):
    r = []
    while True:
        x, y = divmod(x, base)
        r.append(y)
        if x == 0:
            return r
 
def nton(a, base1, base2):
    r = [0]
    for x in reversed(a):
        c = _inplace_mul1(r, base1, x, base2)[0]
        if c != 0:
            r.extend(iton(c, base2))
    return r
 
def ntos(a, base):
    return ''.join(map(str, reversed(nton(a, base, 10))))
 
# Compare.
 
def _cmp(x, y):
    return -1 if x < y else 1
 
def ncmp(a, b):
    n = len(a)
    m = len(b)
    if n != m:
        return _cmp(n, m)
    for x, y in zip(reversed(a), reversed(b)):
        if x != y:
            return _cmp(x, y)
    return 0
 
# Shift.
 
def nshl(a, n):
    return _zeroreduce([0] * n + a[:])
 
def nshr(a, n):
    return _zeroextend(a[n:], 1)
 
# Negate.
 
def nneg(a, s):
    return _fixsign(a.copy(), not s)
 
# Add.
 
def _add(a, b, base, op):
    n = 1 + max(len(a), len(b))
    r = _zeroextend(a, n)
    b = _zeroextend(b, n)
    c = 0
    for i in range(n):
        c, r[i] = divmod(op(r[i], b[i] + abs(c)), base)
    return _zeroreduce(r)
 
def nadd(a, s, b, t, base):
    if s == t:
        op = operator.add
    else:
        op = operator.sub
        if ncmp(a, b) < 0:
            a, b, s = b, a, t
    return _fixsign(_add(a, b, base, op), s)
 
# Subtract.
 
def nsub(a, s, b, t, base):
    return nadd(a, s, b, _fixsign(b, not t)[1], base)
 
# Multiply.
 
def _mul1(a, x, base):
    n = 1 + len(a)
    return _zeroreduce(_inplace_mul1(_zeroextend(a, n), x, 0, base)[1])
 
def _mul(a, b, base):
    if _iszero(b):
        return [0]
    return _add(_mul1(a, b[0], base), _mul(nshl(a, 1), nshr(b, 1), base),
                base, operator.add)
 
def nmul(a, s, b, t, base):
    return _fixsign(_mul(a, b, base), s != t)
 
# Divide.
 
def _div2(a, b):
    if len(a) < len(b):
        q, r = [0], a
    else:
        q, r = _div2(a, nshl(b, 1))
        q = nshl(q, 1)
        if ncmp(r, b) >= 0:
            q = _add(q, [1], 2, operator.add)
            r = _add(r, b, 2, operator.sub)
    return q, r
 
def _div(a, b, base):
    a = nton(a, base, 2)
    b = nton(b, base, 2)
    q, r = _div2(a, b)
    q = nton(q, 2, base)
    r = nton(r, 2, base)
    return q, r
 
def ndiv(a, s, b, t, base):
    if _iszero(b):
        raise ZeroDivisionError
    q, r = _div(a, b, base)
    if s != t and not _iszero(r):
        q = _add(q, [1], base, operator.add)
        r = _add(b, r, base, operator.sub)
    return _fixsign(q, s != t), _fixsign(r, t)
 
# Test.
 
def _test(n, z, op, nozeroy=False):
    for x, y in itertools.product(z, repeat=2):
        if not nozeroy or y != 0:
            u = int(str(op(Num(x), Num(y))))
            v = op(x, y)
            assert u == v, f'{op}({x}, {y})\nExpected:{v}, Got:{u}'
    z = Num.base ** 3
    for _ in range(n):
        x = random.randint(-z, z)
        y = random.randint(-z, z)
        if not nozeroy or y != 0:
            u = int(str(op(Num(x), Num(y))))
            v = op(x, y)
            assert u == v, f'{op}({x}, {y})\nExpected:{v}, Got:{u}'
 
n = 1000
w = Num.base
z = (0, 1, -1, w, -w, w-1, 1-w)
_test(n, z, operator.add)
_test(n, z, operator.sub)
_test(n, z, operator.mul)
_test(n, z, operator.floordiv, True)
_test(n, z, operator.mod, True)
 
def show1(a):
    print(repr(a), a, sep='; ')
 
def _show(z, op, nozeroy=False):
    print(op)
    for x, y in itertools.product(z, repeat=2):
        if not nozeroy or y != 0:
            show1(op(Num(x), Num(y)))
 
w = Num.base
z = (0, 1, -1, w, -w)
_show(z, operator.add)
_show(z, operator.sub)
_show(z, operator.mul)
_show(z, operator.floordiv, True)
_show(z, operator.mod, True)
 
def factorial(n):
    r = Num(1)
    for i in range(2, n+1):
        r *= Num(i)
    return r
 
print(factorial)
for i in range(5):
    show1(factorial(10*i))

IyBAZmlsZSBiaWdudW0uMDAyLnB5CiMgQGluZ3JvdXAgZXhwZXJpbWVudGFsCiMgRXhwZXJpbWVudGFsIGJpZy1udW1iZXIgcmVwcmVzZW50YXRpb24uCiMgKHNpZ24gbWFnbml0dWRlIHZlcnNpb24pCiMgQGRhdGUgMTEvMjIvMjAyNAoKaW1wb3J0IGl0ZXJ0b29scwppbXBvcnQgb3BlcmF0b3IKaW1wb3J0IHJhbmRvbQoKY2xhc3MgTnVtOgogICAgYmFzZSA9IDIqKjgKCiAgICBkZWYgX19pbml0X18oc2VsZiwgaW5pdCk6CiAgICAgICAgaWYgaXNpbnN0YW5jZShpbml0LCB0dXBsZSk6CiAgICAgICAgICAgIHNlbGYubSA9IGluaXRbMF0KICAgICAgICAgICAgc2VsZi5zID0gaW5pdFsxXQogICAgICAgIGVsc2U6CiAgICAgICAgICAgIHNlbGYubSA9IGl0b24oYWJzKGluaXQpLCBOdW0uYmFzZSkKICAgICAgICAgICAgc2VsZi5zID0gaW5pdCA8IDAKCiAgICBkZWYgX19uZWdfXyhzZWxmKToKICAgICAgICByZXR1cm4gTnVtKG5uZWcoc2VsZi5tLCBzZWxmLnMpKQoKICAgIGRlZiBfX2FkZF9fKHNlbGYsIG90aGVyKToKICAgICAgICByZXR1cm4gTnVtKG5hZGQoc2VsZi5tLCBzZWxmLnMsIG90aGVyLm0sIG90aGVyLnMsIE51bS5iYXNlKSkKCiAgICBkZWYgX19zdWJfXyhzZWxmLCBvdGhlcik6CiAgICAgICAgcmV0dXJuIE51bShuc3ViKHNlbGYubSwgc2VsZi5zLCBvdGhlci5tLCBvdGhlci5zLCBOdW0uYmFzZSkpCgogICAgZGVmIF9fbXVsX18oc2VsZiwgb3RoZXIpOgogICAgICAgIHJldHVybiBOdW0obm11bChzZWxmLm0sIHNlbGYucywgb3RoZXIubSwgb3RoZXIucywgTnVtLmJhc2UpKQoKICAgIGRlZiBfX2Zsb29yZGl2X18oc2VsZiwgb3RoZXIpOgogICAgICAgIHJldHVybiBOdW0obmRpdihzZWxmLm0sIHNlbGYucywgb3RoZXIubSwgb3RoZXIucywgTnVtLmJhc2UpWzBdKQoKICAgIGRlZiBfX21vZF9fKHNlbGYsIG90aGVyKToKICAgICAgICByZXR1cm4gTnVtKG5kaXYoc2VsZi5tLCBzZWxmLnMsIG90aGVyLm0sIG90aGVyLnMsIE51bS5iYXNlKVsxXSkKCiAgICBkZWYgX19zdHJfXyhzZWxmKToKICAgICAgICByZXR1cm4gJy0nICogc2VsZi5zICsgbnRvcyhzZWxmLm0sIE51bS5iYXNlKQoKICAgIGRlZiBfX3JlcHJfXyhzZWxmKToKICAgICAgICByZXR1cm4gc3RyKChzZWxmLm0sIHNlbGYucykpCgojIFV0aWxpdHkuCgpkZWYgX2lzemVybyhhKToKICAgIHJldHVybiBsZW4oYSkgPT0gMSBhbmQgYVswXSA9PSAwCgpkZWYgX2ZpeHNpZ24oYSwgcyk6CiAgICByZXR1cm4gYSwgcyBhbmQgbm90IF9pc3plcm8oYSkKCmRlZiBfemVyb2V4dGVuZChhLCBuKToKICAgIHJldHVybiBhWzpdICsgWzBdICogKG4gLSBsZW4oYSkpCgpkZWYgX3plcm9yZWR1Y2UoYSk6CiAgICBuID0gbGVuKGEpCiAgICB3aGlsZSBuICE9IDEgYW5kIGFbbi0xXSA9PSAwOgogICAgICAgIG4gLT0gMQogICAgcmV0dXJuIGFbOm5dCgpkZWYgX2lucGxhY2VfbXVsMShyLCB4LCBjLCBiYXNlKToKICAgIGZvciBpIGluIHJhbmdlKGxlbihyKSk6CiAgICAgICAgYywgcltpXSA9IGRpdm1vZChyW2ldICogeCArIGMsIGJhc2UpCiAgICByZXR1cm4gYywgcgoKIyBDb252ZXJ0LgoKZGVmIGl0b24oeCwgYmFzZSk6CiAgICByID0gW10KICAgIHdoaWxlIFRydWU6CiAgICAgICAgeCwgeSA9IGRpdm1vZCh4LCBiYXNlKQogICAgICAgIHIuYXBwZW5kKHkpCiAgICAgICAgaWYgeCA9PSAwOgogICAgICAgICAgICByZXR1cm4gcgoKZGVmIG50b24oYSwgYmFzZTEsIGJhc2UyKToKICAgIHIgPSBbMF0KICAgIGZvciB4IGluIHJldmVyc2VkKGEpOgogICAgICAgIGMgPSBfaW5wbGFjZV9tdWwxKHIsIGJhc2UxLCB4LCBiYXNlMilbMF0KICAgICAgICBpZiBjICE9IDA6CiAgICAgICAgICAgIHIuZXh0ZW5kKGl0b24oYywgYmFzZTIpKQogICAgcmV0dXJuIHIKCmRlZiBudG9zKGEsIGJhc2UpOgogICAgcmV0dXJuICcnLmpvaW4obWFwKHN0ciwgcmV2ZXJzZWQobnRvbihhLCBiYXNlLCAxMCkpKSkKCiMgQ29tcGFyZS4KCmRlZiBfY21wKHgsIHkpOgogICAgcmV0dXJuIC0xIGlmIHggPCB5IGVsc2UgMQoKZGVmIG5jbXAoYSwgYik6CiAgICBuID0gbGVuKGEpCiAgICBtID0gbGVuKGIpCiAgICBpZiBuICE9IG06CiAgICAgICAgcmV0dXJuIF9jbXAobiwgbSkKICAgIGZvciB4LCB5IGluIHppcChyZXZlcnNlZChhKSwgcmV2ZXJzZWQoYikpOgogICAgICAgIGlmIHggIT0geToKICAgICAgICAgICAgcmV0dXJuIF9jbXAoeCwgeSkKICAgIHJldHVybiAwCgojIFNoaWZ0LgoKZGVmIG5zaGwoYSwgbik6CiAgICByZXR1cm4gX3plcm9yZWR1Y2UoWzBdICogbiArIGFbOl0pCgpkZWYgbnNocihhLCBuKToKICAgIHJldHVybiBfemVyb2V4dGVuZChhW246XSwgMSkKCiMgTmVnYXRlLgoKZGVmIG5uZWcoYSwgcyk6CiAgICByZXR1cm4gX2ZpeHNpZ24oYS5jb3B5KCksIG5vdCBzKQoKIyBBZGQuCgpkZWYgX2FkZChhLCBiLCBiYXNlLCBvcCk6CiAgICBuID0gMSArIG1heChsZW4oYSksIGxlbihiKSkKICAgIHIgPSBfemVyb2V4dGVuZChhLCBuKQogICAgYiA9IF96ZXJvZXh0ZW5kKGIsIG4pCiAgICBjID0gMAogICAgZm9yIGkgaW4gcmFuZ2Uobik6CiAgICAgICAgYywgcltpXSA9IGRpdm1vZChvcChyW2ldLCBiW2ldICsgYWJzKGMpKSwgYmFzZSkKICAgIHJldHVybiBfemVyb3JlZHVjZShyKQoKZGVmIG5hZGQoYSwgcywgYiwgdCwgYmFzZSk6CiAgICBpZiBzID09IHQ6CiAgICAgICAgb3AgPSBvcGVyYXRvci5hZGQKICAgIGVsc2U6CiAgICAgICAgb3AgPSBvcGVyYXRvci5zdWIKICAgICAgICBpZiBuY21wKGEsIGIpIDwgMDoKICAgICAgICAgICAgYSwgYiwgcyA9IGIsIGEsIHQKICAgIHJldHVybiBfZml4c2lnbihfYWRkKGEsIGIsIGJhc2UsIG9wKSwgcykKCiMgU3VidHJhY3QuCgpkZWYgbnN1YihhLCBzLCBiLCB0LCBiYXNlKToKICAgIHJldHVybiBuYWRkKGEsIHMsIGIsIF9maXhzaWduKGIsIG5vdCB0KVsxXSwgYmFzZSkKCiMgTXVsdGlwbHkuCgpkZWYgX211bDEoYSwgeCwgYmFzZSk6CiAgICBuID0gMSArIGxlbihhKQogICAgcmV0dXJuIF96ZXJvcmVkdWNlKF9pbnBsYWNlX211bDEoX3plcm9leHRlbmQoYSwgbiksIHgsIDAsIGJhc2UpWzFdKQoKZGVmIF9tdWwoYSwgYiwgYmFzZSk6CiAgICBpZiBfaXN6ZXJvKGIpOgogICAgICAgIHJldHVybiBbMF0KICAgIHJldHVybiBfYWRkKF9tdWwxKGEsIGJbMF0sIGJhc2UpLCBfbXVsKG5zaGwoYSwgMSksIG5zaHIoYiwgMSksIGJhc2UpLAogICAgICAgICAgICAgICAgYmFzZSwgb3BlcmF0b3IuYWRkKQoKZGVmIG5tdWwoYSwgcywgYiwgdCwgYmFzZSk6CiAgICByZXR1cm4gX2ZpeHNpZ24oX211bChhLCBiLCBiYXNlKSwgcyAhPSB0KQoKIyBEaXZpZGUuCgpkZWYgX2RpdjIoYSwgYik6CiAgICBpZiBsZW4oYSkgPCBsZW4oYik6CiAgICAgICAgcSwgciA9IFswXSwgYQogICAgZWxzZToKICAgICAgICBxLCByID0gX2RpdjIoYSwgbnNobChiLCAxKSkKICAgICAgICBxID0gbnNobChxLCAxKQogICAgICAgIGlmIG5jbXAociwgYikgPj0gMDoKICAgICAgICAgICAgcSA9IF9hZGQocSwgWzFdLCAyLCBvcGVyYXRvci5hZGQpCiAgICAgICAgICAgIHIgPSBfYWRkKHIsIGIsIDIsIG9wZXJhdG9yLnN1YikKICAgIHJldHVybiBxLCByCgpkZWYgX2RpdihhLCBiLCBiYXNlKToKICAgIGEgPSBudG9uKGEsIGJhc2UsIDIpCiAgICBiID0gbnRvbihiLCBiYXNlLCAyKQogICAgcSwgciA9IF9kaXYyKGEsIGIpCiAgICBxID0gbnRvbihxLCAyLCBiYXNlKQogICAgciA9IG50b24ociwgMiwgYmFzZSkKICAgIHJldHVybiBxLCByCgpkZWYgbmRpdihhLCBzLCBiLCB0LCBiYXNlKToKICAgIGlmIF9pc3plcm8oYik6CiAgICAgICAgcmFpc2UgWmVyb0RpdmlzaW9uRXJyb3IKICAgIHEsIHIgPSBfZGl2KGEsIGIsIGJhc2UpCiAgICBpZiBzICE9IHQgYW5kIG5vdCBfaXN6ZXJvKHIpOgogICAgICAgIHEgPSBfYWRkKHEsIFsxXSwgYmFzZSwgb3BlcmF0b3IuYWRkKQogICAgICAgIHIgPSBfYWRkKGIsIHIsIGJhc2UsIG9wZXJhdG9yLnN1YikKICAgIHJldHVybiBfZml4c2lnbihxLCBzICE9IHQpLCBfZml4c2lnbihyLCB0KQoKIyBUZXN0LgoKZGVmIF90ZXN0KG4sIHosIG9wLCBub3plcm95PUZhbHNlKToKICAgIGZvciB4LCB5IGluIGl0ZXJ0b29scy5wcm9kdWN0KHosIHJlcGVhdD0yKToKICAgICAgICBpZiBub3Qgbm96ZXJveSBvciB5ICE9IDA6CiAgICAgICAgICAgIHUgPSBpbnQoc3RyKG9wKE51bSh4KSwgTnVtKHkpKSkpCiAgICAgICAgICAgIHYgPSBvcCh4LCB5KQogICAgICAgICAgICBhc3NlcnQgdSA9PSB2LCBmJ3tvcH0oe3h9LCB7eX0pXG5FeHBlY3RlZDp7dn0sIEdvdDp7dX0nCiAgICB6ID0gTnVtLmJhc2UgKiogMwogICAgZm9yIF8gaW4gcmFuZ2Uobik6CiAgICAgICAgeCA9IHJhbmRvbS5yYW5kaW50KC16LCB6KQogICAgICAgIHkgPSByYW5kb20ucmFuZGludCgteiwgeikKICAgICAgICBpZiBub3Qgbm96ZXJveSBvciB5ICE9IDA6CiAgICAgICAgICAgIHUgPSBpbnQoc3RyKG9wKE51bSh4KSwgTnVtKHkpKSkpCiAgICAgICAgICAgIHYgPSBvcCh4LCB5KQogICAgICAgICAgICBhc3NlcnQgdSA9PSB2LCBmJ3tvcH0oe3h9LCB7eX0pXG5FeHBlY3RlZDp7dn0sIEdvdDp7dX0nCgpuID0gMTAwMAp3ID0gTnVtLmJhc2UKeiA9ICgwLCAxLCAtMSwgdywgLXcsIHctMSwgMS13KQpfdGVzdChuLCB6LCBvcGVyYXRvci5hZGQpCl90ZXN0KG4sIHosIG9wZXJhdG9yLnN1YikKX3Rlc3Qobiwgeiwgb3BlcmF0b3IubXVsKQpfdGVzdChuLCB6LCBvcGVyYXRvci5mbG9vcmRpdiwgVHJ1ZSkKX3Rlc3Qobiwgeiwgb3BlcmF0b3IubW9kLCBUcnVlKQoKZGVmIHNob3cxKGEpOgogICAgcHJpbnQocmVwcihhKSwgYSwgc2VwPSc7ICcpCgpkZWYgX3Nob3coeiwgb3AsIG5vemVyb3k9RmFsc2UpOgogICAgcHJpbnQob3ApCiAgICBmb3IgeCwgeSBpbiBpdGVydG9vbHMucHJvZHVjdCh6LCByZXBlYXQ9Mik6CiAgICAgICAgaWYgbm90IG5vemVyb3kgb3IgeSAhPSAwOgogICAgICAgICAgICBzaG93MShvcChOdW0oeCksIE51bSh5KSkpCgp3ID0gTnVtLmJhc2UKeiA9ICgwLCAxLCAtMSwgdywgLXcpCl9zaG93KHosIG9wZXJhdG9yLmFkZCkKX3Nob3coeiwgb3BlcmF0b3Iuc3ViKQpfc2hvdyh6LCBvcGVyYXRvci5tdWwpCl9zaG93KHosIG9wZXJhdG9yLmZsb29yZGl2LCBUcnVlKQpfc2hvdyh6LCBvcGVyYXRvci5tb2QsIFRydWUpCgpkZWYgZmFjdG9yaWFsKG4pOgogICAgciA9IE51bSgxKQogICAgZm9yIGkgaW4gcmFuZ2UoMiwgbisxKToKICAgICAgICByICo9IE51bShpKQogICAgcmV0dXJuIHIKCnByaW50KGZhY3RvcmlhbCkKZm9yIGkgaW4gcmFuZ2UoNSk6CiAgICBzaG93MShmYWN0b3JpYWwoMTAqaSkp