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 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. = 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)