Some algebraic types for making calculations

alglib.py

from fractions import Fraction, gcd
from itertools import product as cartprod
def memo(f):
    cache = {}
    def g(*argv, **kwargs):
        args = (argv, tuple(sorted(kwargs.items())))
        if args in cache:
            #print("Using {0} |-> {1}".format(args, cache[args]))
            return cache[args]
        else:
            r = f(*argv, **kwargs)
            #print("Saving {0} |-> {1}".format(args, r))
            cache[args] = r
            return r
    return g

class Matrix:
    class SizeMismatchException(BaseException):
        pass
    __slots__ = ["__data", "__w", "__h"]
    def __getitem__(self, ij):
        return self.__data[ij[0]][ij[1]]
    def submatrix(self, rows, cols):
        d = []
        cols = tuple(cols)
        for i in rows:
            d.append([self.__data[i][j] for j in cols])
        return Matrix(d)
    def __len__(self):
        if self.__w == self.__h:
            return self.__w
        raise Matrix.SizeMismatchException()
    @memo
    def det(self):
        if len(self) == 1:
            return self[0, 0]
        d = 0
        k = 1
        cols = tuple(range(1, self.__w))
        rows = tuple(range(self.__h))
        for i in range(self.__h):
            d += self[i, 0] * self.minor(rows[:i] + rows[i + 1:], cols) * k
            k = -k
        return d
    def minor(self, rows, cols):
        return self.submatrix(rows, cols).det()
    def get_size(self):
        return (self.__h, self.__w)
    @memo
    def basis_lines(self):
        rows = [-1]
        cols = [-1]
        flag = True
        while flag:
            flag = False
            for i in range(rows[-1] + 1, self.__h):
                for j in range(cols[-1] + 1, self.__w):
                    if self.minor(rows[1:] + [i], cols[1:] + [j]):
                        rows.append(i)
                        cols.append(j)
                        flag = True
                        break
            if flag:
                break
        return rows[1:], cols[1:]
    def basis_minor(self):
        return self.minor(*self.basis_lines())
    def basis_submatrix(self):
        return self.submatrix(*self.basis_lines())
    def Rg(self):
        return len(self.basis_lines()[0])
    def __init__(self, data):
        if isinstance(data, Matrix):
            data = data.__data
        self.__data = tuple(map(tuple, data))
        self.__h = len(self.__data)
        self.__w = min(map(len, self.__data))
    def __str__(self):
        return '\n'.join(['\t'.join(map(str, i)) for i in self.__data])
    def __repr__(self):
        return '\n'.join(['\t'.join(map(repr, i)) for i in self.__data])
    def __hash__(self):
        return hash(self.__data)
    def __add__(self, other):
        if not isinstance(other, Matrix):
            if other == 0:
                return self
            return NotImplemented
        ss, os = self.get_size(), other.get_size()
        if ss != os:
            raise Matrix.SizeMismatchException()
        return Matrix((self[i, j] + other[i, j] for j in range(ss[1])) for i in range(ss[0]))
    def __radd__(self, other):
        if not isinstance(other, Matrix):
            if other == 0:
                return self
            return NotImplemented
        return other + self
    def __sub__(self, other):
        if not isinstance(other, Matrix):
            if other == 0:
                return self
            return NotImplemented
        ss, os = self.get_size(), other.get_size()
        if ss != os:
            raise Matrix.SizeMismatchException()
        return Matrix((self[i, j] - other[i, j] for j in range(ss[1])) for i in range(ss[0]))
    def transpose(self):
        ss = self.get_size()
        return Matrix((self[j, i] for j in range(ss[0])) for i in range(ss[1]))
    def matrixmul(self, other):
        ss, os = self.get_size(), other.get_size()
        if ss[1] != os[0]:
            raise Matrix.SizeMismatchException()
        return Matrix((sum(self[i, k] * other[k, j] for k in range(ss[1])) for j in range(os[1])) for i in range(ss[0]))
    def numbermul(self, k):
        ss = self.get_size()
        return Matrix((self[i, j] * k for j in range(ss[1])) for i in range(ss[0]))
    def __truediv__(self, k):
        ss = self.get_size()
        return Matrix((self[i, j] / k for j in range(ss[1])) for i in range(ss[0]))
    def __neg__(self):
        ss = self.get_size()
        return Matrix((-self[i, j] for j in range(ss[1])) for i in range(ss[0]))
    def __mul__(self, other):
        if isinstance(other, Matrix):
            if type(self) == Matrix or type(other) == Matrix:
                return self.matrixmul(other)
            return NotImplemented
        return self.numbermul(other)
    def __rmul__(self, other):
        if isinstance(other, Matrix):
            if type(self) == Matrix or type(other) == Matrix:
                return other.matrixmul(self)
            return NotImplemented
        return self.numbermul(other)
    def __eq__(self, other):
        if not isinstance(other, Matrix):
            return NotImplemented
        return self.__data == other.__data
    def cominor(self, i, j):
        rows = list(range(len(self)))
        cols = list(range(len(self)))
        rows.remove(i)
        cols.remove(j)
        return self.minor(rows, cols)
    @memo
    def adjmatrix(self):
        l = len(self)
        return Matrix((((-1) ** (i + j) * self.cominor(j, i) for j in range(l)) for i in range(l)))
    def invmatrix(self):
        return self.adjmatrix() / self.det()
    def __rtruediv__(self, other):
        return other * self.invmatrix()
    @memo
    def __pow__(self, other):
        if type(other) != int:
            return NotImplemented
        if other < 0:
            return self.invmatrix() ** (-other)
        if other == 0:
            return Matrix.identity(len(self))
        if other == 1:
            return self
        if other & 1:
            return self * self ** (other - 1)
        else:
            v = self ** (other >> 1)
            return v * v
    def zeroextend(self, y=0, x=0, yoff=0, xoff=0, getval=None):
        if getval == None:
            getval = lambda i, j: 0
        sz = self.get_size()
        return Matrix((self[i - yoff, j - xoff] if (0 <= i - yoff < sz[0] and 0 <= j - xoff < sz[1]) else getval(i,j) for j in range(x + sz[1])) for i in range(y + sz[0]))
    def wcat(self, other):
        if not isinstance(other, Matrix):
            return NotImplemented
        ss, os = self.get_size(), other.get_size()
        if ss[0] != os[0]:
            raise Matrix.SizeMismatchException()
        return other.zeroextend(x=ss[1], xoff=ss[1], getval=lambda i, j: self[i, j])
    def hcat(self, other):
        if not isinstance(other, Matrix):
            return NotImplemented
        ss, os = self.get_size(), other.get_size()
        if ss[1] != os[1]:
            raise Matrix.SizeMismatchException()
        return other.zeroextend(yoff=ss[0], getval=lambda i, j: self[i, j])
    def substrow(self, i, row):
        if not isinstance(row, Matrix):
            return NotImplemented
        sz = self.get_size()
        return row.zeroextend(yoff=i, y=sz[0] - row.get_size()[0], getval=lambda i, j: self[i, j])
    def substcol(self, j, col):
        if not isinstance(col, Matrix):
            return NotImplemented
        sz = self.get_size()
        return col.zeroextend(xoff=j, x=sz[1] - col.get_size()[1], getval=lambda i, j: self[i, j])
    def rowlc(self, weights):
        sz = self.get_size()
        a = tuple(range(sz[1]))
        return sum(self.submatrix([i], a) * weights[i] for i in range(sz[0]))
    def collc(self, weights):
        sz = self.get_size()
        a = tuple(range(sz[0]))
        return sum(self.submatrix(a, [j]) * weights[j] for j in range(sz[1]))
    def roweltrans(self, d):
        m = self
        for k in d:
            m = m.substrow(k, self.rowlc(tuple(d[k])))
        return m
    def coleltrans(self, d):
        m = self
        for k in d:
            m = m.substcol(k, self.collc(tuple(d[k])))
        return m
    def identity(n):
        if isinstance(n, Matrix):
            raise TypeError('Matrix.identity() is static')
        n = int(n)
        if n < 1:
            raise BaseException(f'Cannot create matrix with size {n}x{n}')
        return Matrix([[1]]).zeroextend(n - 1, n - 1, getval=lambda i, j: int(i == j))
    def scale_rows(self, coefs):
        s = self.get_size()
        return Matrix((self[i, j] * coefs[i] for j in range(s[1])) for i in range(s[0]))
    def scale_cols(self, coefs):
        s = self.get_size()
        return Matrix((self[i, j] * coefs[j] for j in range(s[1])) for i in range(s[0]))
    def diagonal(values):
        return Matrix.identity(len(values)).scale_rows(values)
    def mapvalues(self, f):
        s = self.get_size()
        return Matrix((f(self[i, j]) for j in range(s[1])) for i in range(s[0]))
    def stairs_next_rowlc(self):
        i, j = 0, 0
        s = self.get_size()
        while i < s[0] and j < s[1]:
            if not self[i, j]:
                for k in range(i + 1, s[0]):
                    if self[k, j]:
                        wi = [0] * s[0]
                        wi[i] = 1
                        wi[k] = 1
                        return {i: wi}
                j = j + 1
                continue
            w = {}
            for k in range(i + 1, s[0]):
                if self[k, j]:
                    wk = [0] * s[0]
                    wk[k] = 1
                    wk[i] = -self[k, j] / self[i, j]
                    w[k] = wk
            if len(w) == 0:
                i, j = i + 1, j + 1
                continue
            return w
        return None
    def stairs_row(self, rightskip=0):
        s = self.get_size()
        left = self.submatrix(range(s[0]), range(s[1] - rightskip))
        while True:
            w = left.stairs_next_rowlc()
            if w == None:
                return self
            self = self.roweltrans(w)
            left = left.roweltrans(w)
    def stairs_col(self, bottomskip=0):
        return self.transpose().stairs_row(bottomskip).transpose()
    def canon_row(self, rightskip=0):
        self = self.stairs_row(rightskip)
        h, w = self.get_size()
        if not w * h:
            return self
        coefs = [self[i, 0] for i in range(h)]
        leadidx = [-1] * h
        for i in range(h):
            for j in range(w):
                if self[i, j]:
                    coefs[i] = self[i, j] ** -1
                    leadidx[i] = j
                    break
        self = self.scale_rows(coefs)
        for i in range(h):
            j = leadidx[i]
            if j == -1:
                continue
            d = dict()
            for k in range(i):
                if self[k, j]:
                    coefs = [0] * h
                    coefs[k] = 1
                    coefs[i] = -self[k, j]
                    d[k] = coefs
            if len(d):
                self = self.roweltrans(d)
        return self
    def canon_col(self, bottomskip=0):
        return self.transpose().canon_row(bottomskip).transpose()
    def stairs_row_trace(self, rightskip=0):
        s = self.get_size()
        left = self.submatrix(range(s[0]), range(s[1] - rightskip))
        t = [(None, self)]
        while True:
            w = left.stairs_next_rowlc()
            if w == None:
                return t
            self = self.roweltrans(w)
            left = left.roweltrans(w)
            t.append((w, self))
    def stairs_row_ext(self, right):
            return self.wcat(right).stairs_row(right.get_size()[1])
    def stairs_row_ext_trace(self, right):
            return self.wcat(right).stairs_row_trace(right.get_size()[1])
    def getcol(self, j):
        return self.submatrix(range(self.get_size()[0]), [j])
    def getrow(self, i):
        return self.submatrix([i], range(self.get_size()[1]))
    def orthogonalize_cols(self):
        a = [Vector(self.getcol(j)) for j in range(self.get_size()[1])]
        b = []
        b.append(a[0])
        for ai in a[1:]:
            bi = ai
            for j in range(len(b)):
                bi -= ai.projectto(b[j])
            b.append(bi)
        B = b[0]
        for bi in b[1:]:
            B = B.wcat(bi)
        return B
    def orthogonalize_rows(self):
        return self.transpose().orthogonalize_cols().transpose()
    def orthonorm_cols(self):  # Warning! 0.5 power
        a = [Vector(self.getcol(j)) for j in range(self.get_size()[1])]
        b = []
        b.append(a[0])
        for ai in a[1:]:
            bi = ai
            for j in range(len(b)):
                bi -= ai.projectto(b[j])
            b.append(bi)
        B = b[0] / abs(b[0])
        for bi in b[1:]:
            B = B.wcat(bi / abs(bi))
        return B
    def orthonorm_rows(self):
        return self.transpose().orthonorm_cols().transpose()
    def QR(self):
        Q = self.orthonorm_cols()
        R = Q.transpose() * self
        return Q, R
    def stas(S, A):
        return S.transpose() * A * S
    def fss(self):
        self = self.canon_row()
        s = dict()
        h, w = self.get_size()
        for i in range(h):
            for j in range(w):
                if self[i, j]:
                    for k in range(j + 1, w):
                        if self[i, k]:
                            if k not in s:
                                s[k] = [0] * w
                                s[k][k] = 1
                            s[k][j] = -self[i, k]
                    break
        return tuple(s.values())

class Vector(Matrix):
    def __init__(self, data):
        if isinstance(data, Matrix):
            super().__init__(data)
            return
        super().__init__((i,) for i in data)
    def __add__(self, other):
        return Vector(super().__add__(other))
    def __sub__(self, other):
        return Vector(super().__sub__(other))
    def __truediv__(self, k):
        return Vector(super().__truediv__(k))
    def __mul__(self, other):
        if isinstance(self, Vector) and isinstance(other, Vector): #Scalar multiply
            return self.transpose().matrixmul(other)[0, 0]
        return Vector(super().__mul__(other))
    def __rmul__(self, other):
        if isinstance(self, Vector) and isinstance(other, Vector): #Scalar multiply
            return self.matrixmul(other.transpose())[0, 0]
        return Vector(super().__rmul__(other))
    def __neg__(self, k):
        return Vector(super().__neg__())
    def __getitem__(self, i):
        if type(i) == tuple:
            return super().__getitem__(i)
        return super().__getitem__((i, 0))
    def __iter__(self):
        for i in range(len(self)):
            yield self[i]
    def __len__(self):
        return self.get_size()[0]
    def __abs__(self):
        return (self * self) ** .5
    def __str__(self):
        return '({0})'.format(', '.join(map(str, self)))
    def __repr__(self):
        return '({0})'.format(', '.join(map(repr, self)))
    def projectto(self, other):
        if not isinstance(other, Vector):
            other = Vector(other)
        return (self * other) * other / (other * other)

class Vec3(Vector):
    def __init__(self, x, y, z):
        super().__init__((x, y, z))
    def __getattribute__(self, i):
        if i not in tuple('xyz'):
            return super().__getattribute__(i)
        return self[i]
    def __getitem__(self, i):
        if i not in tuple('xyz'):
            return super().__getitem__(i)
        return self['xyz'.index(i)]
    def triplemul(v1, v2, v3):
        return Matrix((v1, v2, v3)).det()
    def __xor__(self, other):  # Vector product
        return Vec3.triplemul((Vec3(1, 0, 0), Vec3(0, 1, 0), Vec3(0, 0, 1)), self, other)

class Polynom():
    # TODO rewrite working with strings, allow multicharacter variables
    slots = ['names', 'karr', 'field']
    def __init__(self, names, karr=None):
        if isinstance(names, Polynom):
            names, karr = names.names, names.karr
        self.names = tuple(names)
        if karr == None:
            karr = [1] * len(self.names)
        self.karr = tuple(karr)
        self.field = None
    def with_field(self, field):
        p = Polynom(self.names, map(field, self.karr))
        p.field = field
        return p
    def inherit_field(self, *ancestors):
        r = self
        for a in ancestors:
            if a.field != None:
                r = r.with_field(a.field)
        return r
    def from_number(n):
        return Polynom(('',), (n,))
    def to_dict(self):
        try:
            return {self.names[i]: self.karr[i] for i in range(len(self.names))}
        except:
            print(self.names)
            raise
    def _nom_key(nom):
        nom = nom[0].split()
        return str(999999 - len(nom)) + ' '.join(sorted(nom))
    def _nom_fmt(nom):
        m = [str(nom[1])]
        l = ''
        lc = 0
        for v in sorted(nom[0].split()):
            if v == '':
                continue
            if v == l:
                lc += 1
            else:
                if lc == 0:
                    pass
                elif lc == 1:
                    m.append(str(l))
                else:
                    m.append('{0} ** {1}'.format(l, lc))
                l = v
                lc = 1
        if lc == 0:
            pass
        elif lc == 1:
            m.append(str(l))
        else:
            m.append('{0} ** {1}'.format(l, lc))
        return ' * '.join(m)
    def __str__(self):
        if len(self.names) == 0:
            return '(0)'
        d = self.to_dict()
        return ' + '.join(['({0})'.format(Polynom._nom_fmt(nom)) for nom in sorted(d.items(), key=Polynom._nom_key)])
    def __repr__(self):
        return str(self)
    def __neg__(self):
        return Polynom(self.names, map(lambda x: -x, self.karr)).inherit_field(self)
    def remzeros(self):
        karr, names = [], []
        for i in range(len(self.names)):
            if self.karr[i]:
                karr.append(self.karr[i])
                names.append(self.names[i])
        return Polynom(names, karr).inherit_field(self)
    def reduced(self):
        return self // min(self.to_dict().items(), key=Polynom._nom_key)[1]
    def __add__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        names = tuple(set(self.names) | set(other.names))
        sd = self.to_dict()
        od = other.to_dict()
        karr = [sd.get(i, 0) + od.get(i, 0) for i in names]
        return Polynom(names, karr).remzeros().inherit_field(self, other)
    def __radd__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        return other + self
    def __sub__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        return self + (-other)
    def __rsub__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        return other + (-self)
    def __mul__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        d = {}
        sd = self.to_dict()
        od = other.to_dict()
        for i in self.names:
            for j in other.names:
                k = ' '.join(tuple(sorted(i.split() + j.split())))
                a = sd[i] * od[j]
                d[k] = d.get(k, 0) + a
        names, karr = [], []
        for (k, v) in d.items():
            if v != 0:
                names.append(k)
                karr.append(v)
        return Polynom(names, karr).inherit_field(self, other)
    def __rmul__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        d = {}
        sd = self.to_dict()
        od = other.to_dict()
        for j in self.names:
            for i in other.names:
                k = ' '.join(tuple(sorted(i.split() + j.split())))
                a = od[i] * sd[j]
                d[k] = d.get(k, 0) + a
        names, karr = [], []
        for (k, v) in d.items():
            if v != 0:
                names.append(k)
                karr.append(v)
        return Polynom(names, karr).inherit_field(self, other)
    def mul(left, right):
        if type(left) == Polynom:
            return left.__mul__(right)
        if type(right) == Polynom:
            return right.__rmul__(left)
        return Polynom.from_number(left) * Polynom.from_number(right)
    def trypack(self):
        if type(self) != Polynom:
            return self
        self = self.remzeros()
        elems = []
        for i in range(len(self.names)):
            p = Polynom((self.names[i],))
            v = self.karr[i]
            elems.append(v * p) # v * p = v.__mul__(p) ?? p.__rmul__(v)
        if len(elems) == 0:
            return Polynom([])
        return sum(elems)
    def substvalue(self, v, n):
        if type(self) != Polynom:
            return self
        if len(v) != 1:
            raise BaseException('Invalid variable substitution')
        acc = Polynom([])
        for i in range(len(self.names)):
            k = n ** self.names[i].count(v)
            acc += Polynom([self.names[i].replace(v, '')], [k * self.karr[i]])
        return acc
    def constant_as_number(self):
        s = self.remzeros()
        if not s:
            if '' not in self.names:
                return 0
            return self.to_dict()['']
        if '' not in s.names or len(s.names) != 1:
            raise BaseException('Not constant')
        return s.karr[0]
    def __pow__(self, other):
        if type(other) != int or other < 0:
            return NotImplemented()
        if other == 0:
            return Polynom.from_number(1)
        if other % 2 == 0:
            t = self ** (other // 2)
            return t * t
        else:
            return self * self ** (other - 1)
    def __floordiv__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other)
        d = {}
        sd = self.to_dict()
        od = other.to_dict()
        sn = sorted(map(lambda x: x.split(), sd), key=lambda x: (len(x), x), reverse=True)
        on = sorted(map(lambda x: x.split(), od), key=lambda x: (len(x), x), reverse=True)
        if len(on) == 0:
            raise ZeroDivisionError('"{0}" is zero'.format(other))
        if len(sn) == 0:
            return Polynom.from_number(0)
        if len(sn[0]) < len(on[0]):
            return Polynom.from_number(0)
        l = len(on[0])
        p1, p2 = sn[0][:l], sn[0][l:]
        if p1 != on[0]:
            print(p1, '\n', on[0])
            return NotImplemented
        r = Polynom((' '.join(p2),), (sd[' '.join(sn[0])] / od[' '.join(on[0])],))
        if r.karr[0]:
            r += (self - other * r) // other
        return r
    def __mod__(self, other):
        return self - other * (self // other)
    def __hash__(self):
        return hash(tuple(self.to_dict().items()))
    def __eq__(self, other):
        if type(other) != Polynom:
            other = Polynom.from_number(other).inherit_field(self)
        return self.to_dict() == other.to_dict()
    def __bool__(self):
        return bool(len(self.karr))
class F2:
    __slots__ = ('value',)
    def __init__(self, value):
        self.value = int(value) & 1
    def __int__(self):
        return self.value
    def __str__(self):
        return str(self.value)
    def __repr__(self):
        return repr(self.value)
    def __add__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value ^ other.value)
    def __radd__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value ^ other.value)
    def __sub__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value ^ other.value)
    def __rsub__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value ^ other.value)
    def __mul__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value & other.value)
    def __rmul__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value & other.value)
    def __neg__(self):
        return self
    def __bool__(self):
        return bool(self.value)
    def __truediv__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value // other.value)
    def __floordiv__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value // other.value)
    def __mod__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return F2(self.value % other.value)
    def __pow__(self, other):
        return self
    def __hash__(self):
        return hash(self.value)
    def __eq__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return self.value == other.value
    def __lt__(self, other):
        if type(other) != F2:
            if type(other) == int:
                other = F2(other)
            else:
                return NotImplemented
        return self.value < other.value

def gcd(a, b):
    while b:
        a %= b
        a, b = b, a
    return a

def extgcd(a, b):
    D = []
    while b:
        D.append(a // b)
        a %= b
        a, b = b, a
        # R = [a, b, a % b,...]
        # R[i] = R[i - 2] % R[i - 1] = R[i - 2] - R[i - 2] // R[i - 1] * R[i - 1] = R[i - 2] - D[i] * R[i - 1]
        # R[i] = R[i - 2] - D[i] * R[i - 3] + D[i] * D[i - 1] * R[i - 2] = ...
    x, y = 1, 1
    for d in reversed(D):
        x, y = y - d * x, x
    return y, x

def modinv(val, ford):
    return extgcd(val, ford)[0] % ford

class __ZnRingBase:
    __slots__ = ('value',)
    _numtype = int
    _one = 1
    def __init__(self, value):
        if type(value) == type(self):
            value = value.value
        self.value = self._numtype(value)
    def __int__(self):
        return int(self.value)
    def __str__(self):
        return str(self.value)
    def __repr__(self):
        return repr(self.value)
    def __add__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value + other.value)
    def __radd__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value + other.value)
    def __sub__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value - other.value)
    def __rsub__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(-self.value + other.value)
    def __mul__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value * other.value)
    def __rmul__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value * other.value)
    def __neg__(self):
        return type(self)(-self.value)
    def __bool__(self):
        return bool(self.value)
    def __natural_0():
        i = 0
        while True:
            yield i
            i += 1
    @memo
    def __truediv__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        if not other.value:
            raise ZeroDivisionError()
        for i in type(self).__natural_0():
            if other * i == self:
                return type(self)(i)
    def __rtruediv__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return other.__truediv__(self)
    def __floordiv__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value // other.value)
    def __mod__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                other = type(self)(other)
            else:
                return NotImplemented
        return type(self)(self.value % other.value)
    @memo
    def __pow__(self, other):
        if type(other) != int:
            return NotImplemented
        if other == 0 and self.value != 0:
            return type(self)(1)
        if other < 0:
            return (self._one / self) ** (-other)
        if other == 1:
            return self
        if other & 1:
            return self * self ** (other - 1)
        else:
            v = self ** (other >> 1)
            return v * v
    def __hash__(self):
        return hash(self.value)
    def __eq__(self, other):
        if type(other) != type(self):
            return False
        return self.value == other.value
    def __lt__(self, other):
        if type(other) != type(self):
            if type(other) == self._numtype:
                oldval = other
                other = type(self)(other)
                if other.value != oldval:
                    return NotImplemented
            else:
                return NotImplemented
        return self.value < other.value

@memo
def ZnRing(n, numtype=int, one=1):
    class Zn(__ZnRingBase):
        def __init__(self, value):
            if isinstance(value, Zn):
                value = value.value
            super().__init__(self._numtype(value) % n)
        def __truediv__(self, value):
            try:
                if isinstance(value, Zn):
                    value = value.value
                i = modinv(value, n)
                if (value * i) % n == self._one:
                    return self * Zn(i)
                raise ZeroDivisionError(f'Value {value} is not invertible in {numtype} / <{n}>')
            except BaseException:
                raise
            return super().__truediv__(value)
    Zn._numtype = numtype
    Zn._one = one
    return Zn

def product(it):
    acc = 1
    flag = False
    for i in it:
        if flag:
            acc *= i
        else:
            acc = i
            flag = True
    return acc

def nAk(n, k):
    return product(range(n - k + 1, n + 1))

def nCk(n, k):
    return nAk(n, k) // product(range(1, k + 1))

def chinese(*rapairs):
    r_arr, a_arr = zip(*rapairs)
    M = product(a_arr)
    x = 0
    for r, a in rapairs:
        Mi = M // a
        x += r * Mi * modinv(Mi, a)
    return x % M

class PrimesGenerator:
    def __init__(self, bound=2):
        self._primes = [2]
        self._true_bound = 2
        self.update_bound(bound)
    def update_bound(self, bound):
        old_bound = self._true_bound
        self._bound = bound
        if bound <= old_bound:
            return
        self._true_bound = bound
        slen = bound - old_bound
        sieve = [True] * slen
        off = -(old_bound + 1)
        for p in self._primes:
            for i in range(off % p, slen, p):
                sieve[i] = False
        for k in range(slen):
            if sieve[k]:
                p = k - off
                self._primes.append(p)
                for i in range(off % p, slen, p):
                    sieve[i] = False
    def __iter__(self):
        return iter(self.generator())
    def generator(self):
        for i in self._primes:
            if i > self._bound:
                break
            yield i
    def get_bound(self):
        return self._bound


class __EllipticGroupBase:
    __slots__ = ['x', 'y']
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __iter__(self):
        yield self.x
        yield self.y
    def is_zero(self):
        return self.x is None and self.y is None
    def get_order(self, limit=10000):
        acc = self
        for i in range(limit):
            if acc.is_zero():
                return i + 1
            acc += self
        return 0
    def __eq__(self, other):
        if self.is_zero():
            return other.is_zero()
        return self.x == other.x and self.y == other.y
    def __neg__(self):
        if self.is_zero():
            return self
        return type(self)(self.x, -self.y)
    def __sub__(self, other):
        return self + (-other)
    def __mul__(self, n):
        if n < 0:
            return (-self) * (-n)
        if n == 0:
            return type(self).ZERO
        if n == 1:
            return self
        if n & 1:
            return self + self * (n - 1)
        else:
            v = self * (n >> 1)
            return v + v
    def __rmul__(self, n):
        return self.__mul__(n)
    def __bool__(self):
        return not self.is_zero()
    def __str__(self):
        if self.is_zero():
            return "O"
        return "(({0}), ({1}))".format(*self)
    def __repr__(self):
        if self.is_zero():
            return "O"
        return "(({0}), ({1}))".format(*map(repr, self))

@memo
def EllipticGroup(a, b, field=None):
    if field is None:
        field = type(a)
    class EG(__EllipticGroupBase):
        ZERO = None
        def get_field(*args):
            return field
        def get_a(*args):
            return a
        def get_b(*args):
            return b
        def __init__(self, *args):
            if len(args) == 1:
                args = args[0]
            notzero = True
            try:
                x, y = args
            except TypeError:
                x, y = None, None
            if x is not None and y is not None:
                x, y = field(x), field(y)
            super().__init__(x, y)
        def __add__(self, other):
            if self.is_zero():
                return other
            if other.is_zero():
                return self
            x1, y1 = self
            x2, y2 = other
            if self == other:
                if y1 == 0:
                    return type(self).ZERO
                k = (3 * x1 ** 2 + a) / (2 * y1)
            else:
                if x1 == x2:
                    return type(self).ZERO
                k = (y1 - y2) / (x1 - x2)
            x3 = k ** 2 - x1 - x2
            y3 = -y1 + k * (x1 - x3)
            return EG(x3, y3)
        def validate(self):
            if self.is_zero():
                return True
            maybe_zero = self.y ** 2 - self.x ** 3 - a * self.x - b
            if not maybe_zero or maybe_zero == 0:
                return True
            if isinstance(maybe_zero, float):
                x = abs(maybe_zero)
                return x < 0.000001 or x < (abs(self.x) + abs(self.y)) * 0.0000000001
            return False
    EG.ZERO = EG(None, None)
    return EG