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,uc1,ud1,us1,uc2,ud2,us2> = PolynomialRing(QQ,7,order='invlex')
I = R.ideal([
  mynumerator((us1^2+uc1^2)-(1))
, mynumerator((ua*us1^2+ud1^2)-(1))
, mynumerator((us2^2+uc2^2)-(1))
, mynumerator((ua*us2^2+ud2^2)-(1))
])

J = I + R.ideal([0
, uc1-uc2
, ud1-ud2
, us1-us2
])

ua = fastfrac(ua)
uc1 = fastfrac(uc1)
ud1 = fastfrac(ud1)
us1 = fastfrac(us1)
uc2 = fastfrac(uc2)
ud2 = fastfrac(ud2)
us2 = fastfrac(us2)

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()
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()
a0 = fastfrac((fastfrac(1)))
a1 = fastfrac((fastfrac(0)))
a2 = fastfrac((fastfrac(2)-ua))
a3 = fastfrac((fastfrac(0)))
a4 = fastfrac((fastfrac(1)-ua))
a6 = fastfrac((fastfrac(0)))
wx1 = (((ud1-fastfrac(1))*(fastfrac(1)-ua)/(uc1*ua-ud1+fastfrac(1)-ua))).reduce().sreduce()
wy1 = ((us1*(fastfrac(1)-ua)*ua/(uc1*ua-ud1+fastfrac(1)-ua))).reduce().sreduce()
print(isidentity(a0*(wy1^2)+a1*(wx1*wy1)+a3*wy1-(((wx1+a2)*wx1+a4)*wx1+a6)))
print(isidentity(uc1-(fastfrac(1)-fastfrac(2)/((wy1/wx1)^2+ua)-fastfrac(2)*(fastfrac(1)-ua)/(((wy1/wx1)^2+ua)*wx1))))
print(isidentity(ud1-(fastfrac(1)-fastfrac(2)*ua/((wy1/wx1)^2+ua))))
print(isidentity(us1-(-fastfrac(2)*wy1/(((wy1/wx1)^2+ua)*wx1))))
wx2 = (((ud2-fastfrac(1))*(fastfrac(1)-ua)/(uc2*ua-ud2+fastfrac(1)-ua))).reduce().sreduce()
wy2 = ((us2*(fastfrac(1)-ua)*ua/(uc2*ua-ud2+fastfrac(1)-ua))).reduce().sreduce()
wx3 = (((ud3-fastfrac(1))*(fastfrac(1)-ua)/(uc3*ua-ud3+fastfrac(1)-ua))).reduce().sreduce()
wy3 = ((us3*(fastfrac(1)-ua)*ua/(uc3*ua-ud3+fastfrac(1)-ua))).reduce().sreduce()
wx4 = (((ud4-fastfrac(1))*(fastfrac(1)-ua)/(uc4*ua-ud4+fastfrac(1)-ua))).reduce().sreduce()
wy4 = ((us4*(fastfrac(1)-ua)*ua/(uc4*ua-ud4+fastfrac(1)-ua))).reduce().sreduce()
slope = ((wy2-wy1)/(wx2-wx1)).reduce().sreduce()
print(isidentity(a0*slope^2+a1*slope-a2-wx1-wx2-wx3))
print(isidentity(slope*(wx1-wx3)-wy1-a1*wx3-a3-wy3))
slope = ((fastfrac(3)*wx1^2+fastfrac(2)*a2*wx1+a4-a1*wy1)/(fastfrac(2)*a0*wy1+a1*wx1+a3)).reduce().sreduce()
print(isidentity(a0*slope^2+a1*slope-a2-wx1-wx1-wx4))
print(isidentity(slope*(wx1-wx4)-wy1-a1*wx4-a3-wy4))
