Skip to content

Instantly share code, notes, and snippets.

@Vftdan
Created Sep 3, 2021
Embed
What would you like to do?
Some algebraic types for making calculations
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment