# @file bignum.003.py
# @ingroup experimental
# Experimental big-number representation.
# (n's complement version)
# @date 11/24/2024
 
import itertools
import operator
import random
 
class Num:
    base = 2**8
 
    def __init__(self, init):
        if isinstance(init, list):
            self.m = init
        else:
            self.m = iton(init, Num.base)
 
    def __neg__(self):
        return Num(nneg(self.m, Num.base))
 
    def __add__(self, other):
        return Num(nadd(self.m, other.m, Num.base))
 
    def __sub__(self, other):
        return Num(nsub(self.m, other.m, Num.base))
 
    def __mul__(self, other):
        return Num(nmul(self.m, other.m, Num.base))
 
    def __floordiv__(self, other):
        return Num(ndiv(self.m, other.m, Num.base)[0])
 
    def __mod__(self, other):
        return Num(ndiv(self.m, other.m, Num.base)[1])
 
    def __str__(self):
        return '-' * _isneg(self.m) + ntos(self.m, Num.base)
 
    def __repr__(self):
        return str(self.m)
 
# Utility.
 
def _isneg(a):
    return a[-1] != 0
 
def _iszero(a):
    return len(a) == 2 and (a[0] | a[1]) == 0
 
def _samesign(a, b):
    return a[-1] == b[-1]
 
def _signextend(a, n):
    return a[:] + [a[-1]] * (n - len(a))
 
def _signreduce(a):
    n = len(a)
    while n != 2 and a[n-1] == a[n-2]:
        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
 
def _complement(a, base):
    return [base-1-x for x in a]
 
def _abs(a, base):
    return nneg(a, base) if _isneg(a) else a
 
# Convert.
 
def iton(z, base):
    r = []
    x = abs(z)
    while True:
        x, y = divmod(x, base)
        r.append(y)
        if x == 0:
            break
    r.append(0)
    return nneg(r, base) if z < 0 else r
 
def nton(a, base1, base2):
    r = [0]
    for x in reversed(a[:-1]):
        c = _inplace_mul1(r, base1, x, base2)[0]
        while c != 0:
            c, y = divmod(c, base2)
            r.append(y)
    r.append(base2-1 if _isneg(a) else 0)
    return r
 
def ntos(a, base):
    return ''.join(map(str, reversed(nton(_abs(a, base), base, 10)[:-1])))
 
# 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 _signreduce([0] * n + a[:])
 
def nshr(a, n):
    n = min(n, len(a)-1)
    return _signextend(a[n:], 2)
 
# Negate.
 
def nneg(a, base):
    return nadd(_complement(a, base), [1, 0], base)
 
# Add.
 
def nadd(a, b, base):
    n = 1 + max(len(a), len(b))
    r = _signextend(a, n)
    b = _signextend(b, n)
    c = 0
    for i in range(n):
        c, r[i] = divmod(r[i] + b[i] + c, base)
    return _signreduce(r)
 
# Subtract.
 
def nsub(a, b, base):
    return nadd(a, nneg(b, base), base)
 
# Multiply.
 
def _mul1(a, x, base):
    n = 1 + len(a)
    return _signreduce(_inplace_mul1(_signextend(a, n), x, 0, base)[1])
 
def _mul(a, b, base):
    if _iszero(b):
        return [0, 0]
    return nadd(_mul1(a, b[0], base), _mul(nshl(a, 1), nshr(b, 1), base),
                base)
 
def nmul(a, b, base):
    r = _mul(a, _abs(b, base), base)
    if _isneg(b):
        r = nneg(r, base)
    return r
 
# Divide.
 
def _div2(a, b):
    if len(a) < len(b):
        q, r = [0, 0], a
    else:
        q, r = _div2(a, nshl(b, 1))
        q = nshl(q, 1)
        if ncmp(r, b) >= 0:
            q = nadd(q, [1, 0], 2)
            r = nsub(r, b, 2)
    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, b, base):
    if _iszero(b):
        raise ZeroDivisionError
    t = _abs(b, base)
    q, r = _div(_abs(a, base), t, base)
    if not _samesign(a, b) and not _iszero(r):
        q = nadd(q, [1, 0], base)
        r = nsub(t, r, base)
    if not _samesign(a, b):
        q = nneg(q, base)
    if _isneg(b):
        r = nneg(r, base)
    return q, r
 
# 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))

IyBAZmlsZSBiaWdudW0uMDAzLnB5CiMgQGluZ3JvdXAgZXhwZXJpbWVudGFsCiMgRXhwZXJpbWVudGFsIGJpZy1udW1iZXIgcmVwcmVzZW50YXRpb24uCiMgKG4ncyBjb21wbGVtZW50IHZlcnNpb24pCiMgQGRhdGUgMTEvMjQvMjAyNAoKaW1wb3J0IGl0ZXJ0b29scwppbXBvcnQgb3BlcmF0b3IKaW1wb3J0IHJhbmRvbQoKY2xhc3MgTnVtOgogICAgYmFzZSA9IDIqKjgKCiAgICBkZWYgX19pbml0X18oc2VsZiwgaW5pdCk6CiAgICAgICAgaWYgaXNpbnN0YW5jZShpbml0LCBsaXN0KToKICAgICAgICAgICAgc2VsZi5tID0gaW5pdAogICAgICAgIGVsc2U6CiAgICAgICAgICAgIHNlbGYubSA9IGl0b24oaW5pdCwgTnVtLmJhc2UpCgogICAgZGVmIF9fbmVnX18oc2VsZik6CiAgICAgICAgcmV0dXJuIE51bShubmVnKHNlbGYubSwgTnVtLmJhc2UpKQoKICAgIGRlZiBfX2FkZF9fKHNlbGYsIG90aGVyKToKICAgICAgICByZXR1cm4gTnVtKG5hZGQoc2VsZi5tLCBvdGhlci5tLCBOdW0uYmFzZSkpCgogICAgZGVmIF9fc3ViX18oc2VsZiwgb3RoZXIpOgogICAgICAgIHJldHVybiBOdW0obnN1YihzZWxmLm0sIG90aGVyLm0sIE51bS5iYXNlKSkKCiAgICBkZWYgX19tdWxfXyhzZWxmLCBvdGhlcik6CiAgICAgICAgcmV0dXJuIE51bShubXVsKHNlbGYubSwgb3RoZXIubSwgTnVtLmJhc2UpKQoKICAgIGRlZiBfX2Zsb29yZGl2X18oc2VsZiwgb3RoZXIpOgogICAgICAgIHJldHVybiBOdW0obmRpdihzZWxmLm0sIG90aGVyLm0sIE51bS5iYXNlKVswXSkKCiAgICBkZWYgX19tb2RfXyhzZWxmLCBvdGhlcik6CiAgICAgICAgcmV0dXJuIE51bShuZGl2KHNlbGYubSwgb3RoZXIubSwgTnVtLmJhc2UpWzFdKQoKICAgIGRlZiBfX3N0cl9fKHNlbGYpOgogICAgICAgIHJldHVybiAnLScgKiBfaXNuZWcoc2VsZi5tKSArIG50b3Moc2VsZi5tLCBOdW0uYmFzZSkKCiAgICBkZWYgX19yZXByX18oc2VsZik6CiAgICAgICAgcmV0dXJuIHN0cihzZWxmLm0pCgojIFV0aWxpdHkuCgpkZWYgX2lzbmVnKGEpOgogICAgcmV0dXJuIGFbLTFdICE9IDAKCmRlZiBfaXN6ZXJvKGEpOgogICAgcmV0dXJuIGxlbihhKSA9PSAyIGFuZCAoYVswXSB8IGFbMV0pID09IDAKCmRlZiBfc2FtZXNpZ24oYSwgYik6CiAgICByZXR1cm4gYVstMV0gPT0gYlstMV0KCmRlZiBfc2lnbmV4dGVuZChhLCBuKToKICAgIHJldHVybiBhWzpdICsgW2FbLTFdXSAqIChuIC0gbGVuKGEpKQoKZGVmIF9zaWducmVkdWNlKGEpOgogICAgbiA9IGxlbihhKQogICAgd2hpbGUgbiAhPSAyIGFuZCBhW24tMV0gPT0gYVtuLTJdOgogICAgICAgIG4gLT0gMQogICAgcmV0dXJuIGFbOm5dCgpkZWYgX2lucGxhY2VfbXVsMShyLCB4LCBjLCBiYXNlKToKICAgIGZvciBpIGluIHJhbmdlKGxlbihyKSk6CiAgICAgICAgYywgcltpXSA9IGRpdm1vZChyW2ldICogeCArIGMsIGJhc2UpCiAgICByZXR1cm4gYywgcgoKZGVmIF9jb21wbGVtZW50KGEsIGJhc2UpOgogICAgcmV0dXJuIFtiYXNlLTEteCBmb3IgeCBpbiBhXQoKZGVmIF9hYnMoYSwgYmFzZSk6CiAgICByZXR1cm4gbm5lZyhhLCBiYXNlKSBpZiBfaXNuZWcoYSkgZWxzZSBhCgojIENvbnZlcnQuCgpkZWYgaXRvbih6LCBiYXNlKToKICAgIHIgPSBbXQogICAgeCA9IGFicyh6KQogICAgd2hpbGUgVHJ1ZToKICAgICAgICB4LCB5ID0gZGl2bW9kKHgsIGJhc2UpCiAgICAgICAgci5hcHBlbmQoeSkKICAgICAgICBpZiB4ID09IDA6CiAgICAgICAgICAgIGJyZWFrCiAgICByLmFwcGVuZCgwKQogICAgcmV0dXJuIG5uZWcociwgYmFzZSkgaWYgeiA8IDAgZWxzZSByCgpkZWYgbnRvbihhLCBiYXNlMSwgYmFzZTIpOgogICAgciA9IFswXQogICAgZm9yIHggaW4gcmV2ZXJzZWQoYVs6LTFdKToKICAgICAgICBjID0gX2lucGxhY2VfbXVsMShyLCBiYXNlMSwgeCwgYmFzZTIpWzBdCiAgICAgICAgd2hpbGUgYyAhPSAwOgogICAgICAgICAgICBjLCB5ID0gZGl2bW9kKGMsIGJhc2UyKQogICAgICAgICAgICByLmFwcGVuZCh5KQogICAgci5hcHBlbmQoYmFzZTItMSBpZiBfaXNuZWcoYSkgZWxzZSAwKQogICAgcmV0dXJuIHIKCmRlZiBudG9zKGEsIGJhc2UpOgogICAgcmV0dXJuICcnLmpvaW4obWFwKHN0ciwgcmV2ZXJzZWQobnRvbihfYWJzKGEsIGJhc2UpLCBiYXNlLCAxMClbOi0xXSkpKQoKIyBDb21wYXJlLgoKZGVmIF9jbXAoeCwgeSk6CiAgICByZXR1cm4gLTEgaWYgeCA8IHkgZWxzZSAxCgpkZWYgbmNtcChhLCBiKToKICAgIG4gPSBsZW4oYSkKICAgIG0gPSBsZW4oYikKICAgIGlmIG4gIT0gbToKICAgICAgICByZXR1cm4gX2NtcChuLCBtKQogICAgZm9yIHgsIHkgaW4gemlwKHJldmVyc2VkKGEpLCByZXZlcnNlZChiKSk6CiAgICAgICAgaWYgeCAhPSB5OgogICAgICAgICAgICByZXR1cm4gX2NtcCh4LCB5KQogICAgcmV0dXJuIDAKCiMgU2hpZnQuCgpkZWYgbnNobChhLCBuKToKICAgIHJldHVybiBfc2lnbnJlZHVjZShbMF0gKiBuICsgYVs6XSkKCmRlZiBuc2hyKGEsIG4pOgogICAgbiA9IG1pbihuLCBsZW4oYSktMSkKICAgIHJldHVybiBfc2lnbmV4dGVuZChhW246XSwgMikKCiMgTmVnYXRlLgoKZGVmIG5uZWcoYSwgYmFzZSk6CiAgICByZXR1cm4gbmFkZChfY29tcGxlbWVudChhLCBiYXNlKSwgWzEsIDBdLCBiYXNlKQoKIyBBZGQuCgpkZWYgbmFkZChhLCBiLCBiYXNlKToKICAgIG4gPSAxICsgbWF4KGxlbihhKSwgbGVuKGIpKQogICAgciA9IF9zaWduZXh0ZW5kKGEsIG4pCiAgICBiID0gX3NpZ25leHRlbmQoYiwgbikKICAgIGMgPSAwCiAgICBmb3IgaSBpbiByYW5nZShuKToKICAgICAgICBjLCByW2ldID0gZGl2bW9kKHJbaV0gKyBiW2ldICsgYywgYmFzZSkKICAgIHJldHVybiBfc2lnbnJlZHVjZShyKQoKIyBTdWJ0cmFjdC4KCmRlZiBuc3ViKGEsIGIsIGJhc2UpOgogICAgcmV0dXJuIG5hZGQoYSwgbm5lZyhiLCBiYXNlKSwgYmFzZSkKCiMgTXVsdGlwbHkuCgpkZWYgX211bDEoYSwgeCwgYmFzZSk6CiAgICBuID0gMSArIGxlbihhKQogICAgcmV0dXJuIF9zaWducmVkdWNlKF9pbnBsYWNlX211bDEoX3NpZ25leHRlbmQoYSwgbiksIHgsIDAsIGJhc2UpWzFdKQoKZGVmIF9tdWwoYSwgYiwgYmFzZSk6CiAgICBpZiBfaXN6ZXJvKGIpOgogICAgICAgIHJldHVybiBbMCwgMF0KICAgIHJldHVybiBuYWRkKF9tdWwxKGEsIGJbMF0sIGJhc2UpLCBfbXVsKG5zaGwoYSwgMSksIG5zaHIoYiwgMSksIGJhc2UpLAogICAgICAgICAgICAgICAgYmFzZSkKCmRlZiBubXVsKGEsIGIsIGJhc2UpOgogICAgciA9IF9tdWwoYSwgX2FicyhiLCBiYXNlKSwgYmFzZSkKICAgIGlmIF9pc25lZyhiKToKICAgICAgICByID0gbm5lZyhyLCBiYXNlKQogICAgcmV0dXJuIHIKCiMgRGl2aWRlLgoKZGVmIF9kaXYyKGEsIGIpOgogICAgaWYgbGVuKGEpIDwgbGVuKGIpOgogICAgICAgIHEsIHIgPSBbMCwgMF0sIGEKICAgIGVsc2U6CiAgICAgICAgcSwgciA9IF9kaXYyKGEsIG5zaGwoYiwgMSkpCiAgICAgICAgcSA9IG5zaGwocSwgMSkKICAgICAgICBpZiBuY21wKHIsIGIpID49IDA6CiAgICAgICAgICAgIHEgPSBuYWRkKHEsIFsxLCAwXSwgMikKICAgICAgICAgICAgciA9IG5zdWIociwgYiwgMikKICAgIHJldHVybiBxLCByCgpkZWYgX2RpdihhLCBiLCBiYXNlKToKICAgIGEgPSBudG9uKGEsIGJhc2UsIDIpCiAgICBiID0gbnRvbihiLCBiYXNlLCAyKQogICAgcSwgciA9IF9kaXYyKGEsIGIpCiAgICBxID0gbnRvbihxLCAyLCBiYXNlKQogICAgciA9IG50b24ociwgMiwgYmFzZSkKICAgIHJldHVybiBxLCByCgpkZWYgbmRpdihhLCBiLCBiYXNlKToKICAgIGlmIF9pc3plcm8oYik6CiAgICAgICAgcmFpc2UgWmVyb0RpdmlzaW9uRXJyb3IKICAgIHQgPSBfYWJzKGIsIGJhc2UpCiAgICBxLCByID0gX2RpdihfYWJzKGEsIGJhc2UpLCB0LCBiYXNlKQogICAgaWYgbm90IF9zYW1lc2lnbihhLCBiKSBhbmQgbm90IF9pc3plcm8ocik6CiAgICAgICAgcSA9IG5hZGQocSwgWzEsIDBdLCBiYXNlKQogICAgICAgIHIgPSBuc3ViKHQsIHIsIGJhc2UpCiAgICBpZiBub3QgX3NhbWVzaWduKGEsIGIpOgogICAgICAgIHEgPSBubmVnKHEsIGJhc2UpCiAgICBpZiBfaXNuZWcoYik6CiAgICAgICAgciA9IG5uZWcociwgYmFzZSkKICAgIHJldHVybiBxLCByCgojIFRlc3QuCgpkZWYgX3Rlc3Qobiwgeiwgb3AsIG5vemVyb3k9RmFsc2UpOgogICAgZm9yIHgsIHkgaW4gaXRlcnRvb2xzLnByb2R1Y3QoeiwgcmVwZWF0PTIpOgogICAgICAgIGlmIG5vdCBub3plcm95IG9yIHkgIT0gMDoKICAgICAgICAgICAgdSA9IGludChzdHIob3AoTnVtKHgpLCBOdW0oeSkpKSkKICAgICAgICAgICAgdiA9IG9wKHgsIHkpCiAgICAgICAgICAgIGFzc2VydCB1ID09IHYsIGYne29wfSh7eH0sIHt5fSlcbkV4cGVjdGVkOnt2fSwgR290Ont1fScKICAgIHogPSBOdW0uYmFzZSAqKiAzCiAgICBmb3IgXyBpbiByYW5nZShuKToKICAgICAgICB4ID0gcmFuZG9tLnJhbmRpbnQoLXosIHopCiAgICAgICAgeSA9IHJhbmRvbS5yYW5kaW50KC16LCB6KQogICAgICAgIGlmIG5vdCBub3plcm95IG9yIHkgIT0gMDoKICAgICAgICAgICAgdSA9IGludChzdHIob3AoTnVtKHgpLCBOdW0oeSkpKSkKICAgICAgICAgICAgdiA9IG9wKHgsIHkpCiAgICAgICAgICAgIGFzc2VydCB1ID09IHYsIGYne29wfSh7eH0sIHt5fSlcbkV4cGVjdGVkOnt2fSwgR290Ont1fScKCm4gPSAxMDAwCncgPSBOdW0uYmFzZQp6ID0gKDAsIDEsIC0xLCB3LCAtdywgdy0xLCAxLXcpCl90ZXN0KG4sIHosIG9wZXJhdG9yLmFkZCkKX3Rlc3Qobiwgeiwgb3BlcmF0b3Iuc3ViKQpfdGVzdChuLCB6LCBvcGVyYXRvci5tdWwpCl90ZXN0KG4sIHosIG9wZXJhdG9yLmZsb29yZGl2LCBUcnVlKQpfdGVzdChuLCB6LCBvcGVyYXRvci5tb2QsIFRydWUpCgpkZWYgc2hvdzEoYSk6CiAgICBwcmludChyZXByKGEpLCBhLCBzZXA9JzsgJykKCmRlZiBfc2hvdyh6LCBvcCwgbm96ZXJveT1GYWxzZSk6CiAgICBwcmludChvcCkKICAgIGZvciB4LCB5IGluIGl0ZXJ0b29scy5wcm9kdWN0KHosIHJlcGVhdD0yKToKICAgICAgICBpZiBub3Qgbm96ZXJveSBvciB5ICE9IDA6CiAgICAgICAgICAgIHNob3cxKG9wKE51bSh4KSwgTnVtKHkpKSkKCncgPSBOdW0uYmFzZQp6ID0gKDAsIDEsIC0xLCB3LCAtdykKX3Nob3coeiwgb3BlcmF0b3IuYWRkKQpfc2hvdyh6LCBvcGVyYXRvci5zdWIpCl9zaG93KHosIG9wZXJhdG9yLm11bCkKX3Nob3coeiwgb3BlcmF0b3IuZmxvb3JkaXYsIFRydWUpCl9zaG93KHosIG9wZXJhdG9yLm1vZCwgVHJ1ZSkKCmRlZiBmYWN0b3JpYWwobik6CiAgICByID0gTnVtKDEpCiAgICBmb3IgaSBpbiByYW5nZSgyLCBuKzEpOgogICAgICAgIHIgKj0gTnVtKGkpCiAgICByZXR1cm4gcgoKcHJpbnQoZmFjdG9yaWFsKQpmb3IgaSBpbiByYW5nZSg1KToKICAgIHNob3cxKGZhY3RvcmlhbCgxMCppKSk=