def mynumerator(x):
  if parent(x) == R:
    return x
  return numerator(x)

class fastfrac:
  def __init__(self,top,bot=1):
    if parent(top) == ZZ or parent(top) == R:
      self.top = R(top)
      self.bot = R(bot)
    elif top.__class__ == fastfrac:
      self.top = top.top
      self.bot = top.bot * bot
    else:
      self.top = R(numerator(top))
      self.bot = R(denominator(top)) * bot
  def reduce(self):
    return fastfrac(self.top / self.bot)
  def sreduce(self):
    return fastfrac(I.reduce(self.top),I.reduce(self.bot))
  def iszero(self):
    return self.top in I and not (self.bot in I)
  def isdoublingzero(self):
    return self.top in J and not (self.bot in J)
  def __add__(self,other):
    if parent(other) == ZZ:
      return fastfrac(self.top + self.bot * other,self.bot)
    if other.__class__ == fastfrac:
      return fastfrac(self.top * other.bot + self.bot * other.top,self.bot * other.bot)
    return NotImplemented
  def __sub__(self,other):
    if parent(other) == ZZ:
      return fastfrac(self.top - self.bot * other,self.bot)
    if other.__class__ == fastfrac:
      return fastfrac(self.top * other.bot - self.bot * other.top,self.bot * other.bot)
    return NotImplemented
  def __neg__(self):
    return fastfrac(-self.top,self.bot)
  def __mul__(self,other):
    if parent(other) == ZZ:
      return fastfrac(self.top * other,self.bot)
    if other.__class__ == fastfrac:
      return fastfrac(self.top * other.top,self.bot * other.bot)
    return NotImplemented
  def __rmul__(self,other):
    return self.__mul__(other)
  def __div__(self,other):
    if parent(other) == ZZ:
      return fastfrac(self.top,self.bot * other)
    if other.__class__ == fastfrac:
      return fastfrac(self.top * other.bot,self.bot * other.top)
    return NotImplemented
  __truediv__ = __div__
  def __pow__(self,other):
    if parent(other) == ZZ:
      return fastfrac(self.top ^ other,self.bot ^ other)
    return NotImplemented

def isidentity(x):
  return x.iszero()

def isdoublingidentity(x):
  return x.isdoublingzero()

R.<ua,uc2,ud2,us2,uc1,ud1,us1,uSC1,uS1,uD1,uC1,uDZ1,uZ1,uSC2,uS2,uD2,uC2,uDZ2,uZ2> = PolynomialRing(QQ,19,order='invlex')
I = R.ideal([
  mynumerator((us1^2+uc1^2)-(1))
, mynumerator((ua*us1^2+ud1^2)-(1))
, mynumerator((us1)-(uS1/uZ1))
, mynumerator((uc1)-(uC1/uZ1))
, mynumerator((ud1)-(uD1/uZ1))
, mynumerator((uSC1)-(uS1*uC1))
, mynumerator((uDZ1)-(uD1*uZ1))
, mynumerator((us2^2+uc2^2)-(1))
, mynumerator((ua*us2^2+ud2^2)-(1))
, mynumerator((us2)-(uS2/uZ2))
, mynumerator((uc2)-(uC2/uZ2))
, mynumerator((ud2)-(uD2/uZ2))
, mynumerator((uSC2)-(uS2*uC2))
, mynumerator((uDZ2)-(uD2*uZ2))
])

J = I + R.ideal([0
, uSC1-uSC2
, uC1-uC2
, uD1-uD2
, uDZ1-uDZ2
, uS1-uS2
, uZ1-uZ2
])

ua = fastfrac(ua)
uc2 = fastfrac(uc2)
ud2 = fastfrac(ud2)
us2 = fastfrac(us2)
uc1 = fastfrac(uc1)
ud1 = fastfrac(ud1)
us1 = fastfrac(us1)
uSC1 = fastfrac(uSC1)
uS1 = fastfrac(uS1)
uD1 = fastfrac(uD1)
uC1 = fastfrac(uC1)
uDZ1 = fastfrac(uDZ1)
uZ1 = fastfrac(uZ1)
uSC2 = fastfrac(uSC2)
uS2 = fastfrac(uS2)
uD2 = fastfrac(uD2)
uC2 = fastfrac(uC2)
uDZ2 = fastfrac(uDZ2)
uZ2 = fastfrac(uZ2)


uS3 = ((uZ1*uC2*uS1*uD2+uD1*uS2*uC1*uZ2))
uC3 = ((uZ1*uC2*uC1*uZ2-uD1*uS2*uS1*uD2))
uD3 = ((uZ1*uD1*uZ2*uD2-ua*uS1*uC1*uS2*uC2))
uZ3 = (((uZ1*uC2)^2+(uD1*uS2)^2))
uSC3 = ((uS3*uC3))
uDZ3 = ((uD3*uZ3))

uc3 = (((uc2*uc1-ud1*us2*us1*ud2)/(uc2^2+(ud1*us2)^2))).reduce()
ud3 = (((ud1*ud2-ua*us1*uc1*us2*uc2)/(uc2^2+(ud1*us2)^2))).reduce()
us3 = (((uc2*us1*ud2+ud1*us2*uc1)/(uc2^2+(ud1*us2)^2))).reduce()

print(isidentity((us3^2+uc3^2)-(fastfrac(1))))
print(isidentity((ua*us3^2+ud3^2)-(fastfrac(1))))
print(isidentity((us3)-(uS3/uZ3)))
print(isidentity((uc3)-(uC3/uZ3)))
print(isidentity((ud3)-(uD3/uZ3)))
print(isidentity((uSC3)-(uS3*uC3)))
print(isidentity((uDZ3)-(uD3*uZ3)))

unified = True
uSC4 = uSC3
uC4 = uC3
uD4 = uD3
uDZ4 = uDZ3
uS4 = uS3
uZ4 = uZ3
uc4 = (((uc1*uc1-ud1*us1*us1*ud1)/(uc1^2+(ud1*us1)^2))).reduce()
ud4 = (((ud1*ud1-ua*us1*uc1*us1*uc1)/(uc1^2+(ud1*us1)^2))).reduce()
us4 = (((uc1*us1*ud1+ud1*us1*uc1)/(uc1^2+(ud1*us1)^2))).reduce()
if unified: unified = isdoublingidentity((us4^2+uc4^2)-(fastfrac(1)))
if unified: unified = isdoublingidentity((ua*us4^2+ud4^2)-(fastfrac(1)))
if unified: unified = isdoublingidentity((us4)-(uS4/uZ4))
if unified: unified = isdoublingidentity((uc4)-(uC4/uZ4))
if unified: unified = isdoublingidentity((ud4)-(uD4/uZ4))
if unified: unified = isdoublingidentity((uSC4)-(uS4*uC4))
if unified: unified = isdoublingidentity((uDZ4)-(uD4*uZ4))
if unified: print("Unified")