from math import gcd, pi
from functools import reduce
from re import sub as re_sub

def newChar(a,b):
    d = gcd(a,b)
    return d and a*b//d

reduc = lambda a,b: a and ((not b and a) or a%b)

def coeff(n,s):
    if s == '':
        return str(n)
    elif n == 1:
        return ''+s
    elif n == -1:
        return '-'+s
    else:
        return str(n)+'*'+s

ID0 = lambda *args : args[0]

TRUE = lambda *args : True

class AbGrpElt:
    def __init__(self,values):
        self.values = {key:reduc(value,key[1]) for key, value in values.items() if reduc(value,key[1])}
        self.STR = tuple(sorted(self.values.items()))
    
    def debug_str(self): #Just in case
        return "AbGrpElt("+str(self.values)+")"
    
    def __lt__(self,other):
        return self.STR < other.STR
    
    def __add__(self,other,preSetDual=None):
        if type(other) != type(self):
            other = AbGrpElt({(AbGrpElt({}),0):other})
            other.dual = other
        intersect = {key:(self.values[key]+other.values[key]) for key in self.values.keys() & other.values.keys()}
        #return AbGrpElt(self.values | other.values | intersect)
        ret = AbGrpElt(self.values | other.values | intersect)
        if preSetDual != None:
            ret.dual = preSetDual
        elif hasattr(self,"dual") and hasattr(other,"dual"):
            ret.dual = other.dual.__add__(self.dual,ret)
        return ret
        
    def __radd__(self,other):
        return self+other
        
    def __mul__(self,other,preSetDual=None):
        if type(other) != type(self):
            #return AbGrpElt({key:value*other for key, value in self.values.items()})
            ret = AbGrpElt({key:value*other for key, value in self.values.items()})
            if type(ret) != type(self):
                return ret ###No explicit dual for scalars
            if preSetDual != None:
                ret.dual = preSetDual
            elif hasattr(self,"dual"):
                ret.dual = self.dual.__mul__(other,ret)
            return ret
        
        #return sum(AbGrpElt({((x[0]+y[0]),newChar(x[1],y[1])):(self.values[x]*other.values[y])}) for x in self.values for y in other.values)
        ret = sum(AbGrpElt({((x[0]+y[0]),newChar(x[1],y[1])):(self.values[x]*other.values[y])}) for x in self.values for y in other.values)
        if type(ret) != type(self):
            return ret ###No explicit dual for scalars
        if preSetDual != None:
            ret.dual = preSetDual
        elif hasattr(self,"dual") and hasattr(other,"dual"):
            ret.dual = other.dual.__mul__(self.dual,ret)
        return ret

    
    def __rmul__(self,other):
        return self*other
    
    def __neg__(self):
        return self*(-1)
    
    def __sub__(self,other):
        return self + -other
    
    def __rsub__(self,other):
        return -self+other
    
    def __pow__(self,n):
        if n==0:
            return AbGrpElt({(AbGrpElt({}),0):1})
        elif n == 1:
            return self
        elif n & 1:
            return self*(self**(n-1))
        else:
            return (self*self)**(n//2)
    
    
    def __eq__(self,other):
        if type(other) != type(self):
            return False
        return self.STR == other.STR
    
    def __hash__(self):
        return hash(self.STR)
    
    def monomial_str(self):
        return '*'.join(key[0]+("**"+str(value) if value != 1 else '') for key,value in self.values.items())
    
    def __repr__(self):
        try:
            return ' + '.join(coeff(value,key[0].monomial_str()) for key, value in self.values.items())
        except AttributeError:
            return self.monomial_str()
        except:
            return self.debug_str()
    
    def __str__(self):
        R = self.__repr__()
        if '+' in R:
            return '(' + R + ')'
        else:
            return R
    
    def modify(self,g=ID0,f=ID0,condition=TRUE):
        return AbGrpElt({f(key,value):g(value,key) for key, value in self.values.items() if condition(key,value)})
    
    def d(self,var,order=1,char=0):
        """DERIVATIVE w.r.t. a given variable"""
        if order <= 0:
            return self
        if type(var) == type(self):
            key0 = tuple(tuple(var.values.keys())[0][0].values.keys())[0]
        elif type(var) == type(''):
            key0 = (var, char)
        elif type(var) == type((1,)):
            key0 = var
        
        f1_g0 = lambda value, key: value-1 if key == key0 else value
        #f1_cond0 = lambda key, value: key == key0
        f1 = lambda key, value : (key[0].modify(f1_g0,ID0,TRUE),key[1])
        g1 = lambda value, key : value*key[0].values[key0]
        cond1 = lambda key, value : key0 in key[0].values.keys()
        return self.modify(g1,f1,cond1).d(var,order-1,char)
    
    def __getitem__(self,other):
        """Differential Operator: corresponds to other.d(self); note: self must be `simple` e.g. x"""
        return other.d(self)
    
    def __call__(self,**formulas):
        """To evaluate - VERY BARBARIC (and lazy) implementation"""
        pattern = "|".join(key for key in formulas.keys())
        res = re_sub(pattern, lambda match: '('+formulas[match.group(0)].__repr__()+')', str(self))
        return eval(res)
    
    def __invert__(self):
        """Returns ~self if `dual` attribute defined; otherwise, returns self"""
        if hasattr(self,"dual"):
            return self.dual
        else:
            return self
    

        
        
    
    


def Ring(char=0,**variables):
    for key, value in variables.items():
        globals()[key] = AbGrpElt({(AbGrpElt({(key,value):1}),char):1})

def createVars(container,char=0,values=0):
    """Creates algebraic variables from an iterable object such as a string"""
    Ring(char,**{x:values for x in container})

### Quaternions ###

class Vec3:
    """Can have length superior to 3, but no less than 3 (zeros are automatically added)"""
    def __init__(self,*xyz):
        self.xyz = xyz + (0,)*(3-len(xyz))
    
    def __repr__(self):
        return "Vec3" + str(self.xyz)
    
    def __iter__(self):
        return self.xyz.__iter__()
    
    def __eq__(self,other):
        return all(s == t for (s,t) in zip(self,other))
    
    def __getitem__(self,n):
        return self.xyz[n]
    
    def __mul__(self,other):
        """Scalar Multiplication"""
        return Vec3(*(s*other for s in self))
    
    def __neg__(self):
        return self*(-1)
    
    def __rmul__(self,other):
        return Vec3(*(other*s for s in self))
    
    def __or__(self,other):
        """Dot Product"""
        return sum(s*t for (s,t) in zip(self,other)) #self.x*other.x + self.y*other.y + self.z*other.z
    
    def __add__(self,other):
        return Vec3(*(s+t for (s,t) in zip(self,other)))
    
    def __sub__(self,other):
        return self + -other
    
    def __and__(self,other):
        return Vec3(*(s*t for (s,t) in zip(self,other)))
    
    def __lshift__(self,n):
        return Vec3(*(self[n:]+self[:n]))
    
    def __rshift__(self,n):
        return self << -n
    
    def __matmul__(self,other):
        """Cross Product (defined in a general way, although only makes sense for length 3)"""
        return ((self << 1)&(other >> 1)) - ((other << 1)&(self >> 1))
    
    def squaredNorm(self):
        return (self|self)
    
    def __abs__(self):
        return (self|self)**0.5
    
    def unit(self):
        """Return a unit vector; `self` must obviously be nonzero"""
        return self*(1/abs(self))
    
    def __call__(self,cplx):
        """Turns a Complex into a Quaternion"""
        return Quaternion(cplx.real*abs(self),self*cplx.imag)
        
        

class Quaternion:
    def __init__(self,time=0,space=Vec3()):
        self.time = time
        self.space = space
    
    def __str__(self):
        return " + ".join(str(n)+L for (n,L) in zip((self.time,*self.space.xyz),('','*I','*J','*K')) if n != 0) or '0'
    
    def __repr__(self):
        S = str(self)
        return '(' + S + ')' if '+' in S else S
    
    def __invert__(self):
        """Conjugate"""
        return Quaternion(self.time,-self.space)
    
    def __neg__(self):
        """Conjugate"""
        return Quaternion(-self.time,-self.space)
    
    def squaredNorm(self):
        return self.time**2 + self.space.squaredNorm()
    
    def __abs__(self):
        return self.squaredNorm()**0.5
    
    def __add__(self,other):
        return Quaternion(self.time+other.time,self.space+other.space) if type(other) == type(self) else Quaternion(self.time+other,self.space)
    
    def __radd__(self,other):
        return self+other
    
    def __sub__(self,other):
        return self+ -other
    
    def __rsub__(self,other):
        return -self + other
    
    def __mul__(self,other):
        if type(other) == type(self):
            return Quaternion(self.time*other.time - (self.space|other.space),other.space.__rmul__(self.time) + self.space.__mul__(other.time) + (self.space @ other.space))
        else:
            return Quaternion(self.time*other,self.space*other)
    
    def __rmul__(self,other):
        if type(other) == type(self):
            print("??? Not supposed to happen???")
            return other*self
        else:
            return self*other
    
    def inverse(self):
        return (~self)*(1/self.squaredNorm())
    
    def __truediv__(self,other):
        if type(other) == type(self):
            return self*(other.inverse())
        else:
            return self*(1/other)
    
    def __rtruediv__(self,other):
        return other*(self.inverse())
    
    def __complex__(self):
        return complex(self.time,abs(self.space))
    
    def unit(self):
        return self/abs(self)
    
    def direction(self):
        return self.space.unit()
    
    def __pow__(self,k):
        return self.direction()(complex(self)**k)
    
    def __getitem__(self,other):
        """Commutator"""
        return Quaternion(0,2*(self.space @ other.space)) # equivalent to self*other - other*self
    
    def __call__(self,other,power=1):
        return (self*other)/self if power == 1 else (self**power)(other)
    
    def __mod__(self,other):
        """SO(3,1)+"""
        return lambda q: (self*q-other)/(other*q+self)
    
    

def Quat(t=0,x=0,y=0,z=0):
    return Quaternion(t,Vec3(x,y,z))

I = Quat(x=1)
J = Quat(y=1)
K = Quat(z=1)

U = Quat(1)

def rotate(direction,radAngle=pi):
    return direction(1j**(radAngle/pi))

def poly(*coeffs):
    return lambda x: sum(a*x**n for (n,a) in enumerate(coeffs)) #sum(coeffs[n]*x**n for n in range(len(coeffs)))

def pore(*coeffs):
    return lambda x: reduce(lambda ret, a: ret * x + a, coeffs)

def cyclot(n):
    return pore(*((1,)*n))


def ENorm(x,y):
    return x*x - x*y + y*y

def secure(x,other, attr="__call__"):
    if hasattr(x,attr):
        return x
    else:
        return x+other

class Moebius():
    def __init__(self,a,b,c,d):
        self.a = a
        self.b = b
        self.c = c
        self.d = d
    
    def abcd(self):
        return (self.a,self.b,self.c,self.d)
    
    def inv_axxd(self):
        """Inverse transformation, provided a,d!=0"""
        (a,b,c,d) = self.abcd()
        B = -(a.inverse())*b*a
        C = -(d.inverse())*c*d
        return Moebius(d,B,C,a)@Moebius((a*d+b*C).inverse(),0,0,(c*B+d*a).inverse())
    
    def inv_0bc0(self):
        """Inverse transformation, provided a,d = 0, b,c != 0"""
        (a,b,c,d) = self.abcd()
        return Moebius(0,c.inverse(),b.inverse(),0)
    
    def inv_abc0(self):
        """Inverse, provided a!=0 (and thus b,c!=0) and d = 0"""
        (a,b,c,d) = self.abcd()
        B = -(a.inverse())*b*a
        C = b.inverse()
        return Moebius(d,B,C,a)@Moebius((b*C).inverse(),0,0,(c*B).inverse())
    
    def inv_0bcd(self):
        """Inverse, provided d!=0 (and thus b,c !=0) and a = 0"""
        (a,b,c,d) = self.abcd()
        B = c.inverse()
        C = -(d.inverse())*c*d
        return Moebius(d,B,C,a)@Moebius((b*C).inverse(),0,0,(c*B).inverse())
    
    def inv_ax0d(self):
        """Inverse, provided c = 0"""
        (a,b,c,d) = self.abcd()
        return Moebius(a.inverse(),-(a.inverse())*b*(d.inverse()),0,d.inverse())
    
    def inv_xxcx(self):
        """Inverse, provided c != 0 (and b-ac^(-1)d != 0)"""
        (a,b,c,d) = self.abcd()
        A = b-a*(c.inverse())*d
        return Moebius(c.inverse(),0,0,1)@Moebius(1,-d,0,1)@Moebius(0,1,1,0)@Moebius(1,-(A.inverse())*a*(c.inverse()),0,1)@Moebius(A.inverse(),0,0,1)
    
    def inverse(self):
        """Inverse, general case, as long as invertible"""
        if abs(self.c) == 0:
            return self.inv_ax0d()
        else:
            return self.inv_xxcx()
    
    def oldinvert(self):
        """Test: inverting the actual Transformation (Not the Matrix); oops merde pour c=0"""
        (a,b,c,d) = self.abcd()
        alph = (~b)*c*(~c) - (~d)*c*(~a)
        coAlph = b*(~b)*c*(~c) + a*(~a)*d*(~d) - b*(~d)*c*(~a) - a*(~c)*d*(~b)
        return Moebius(~c,0,0,c*(~c))@Moebius(1,-d,0,1)@Moebius(0,1,1,0)@Moebius(coAlph*c*(~c),-alph*a*(~c),0,coAlph*c*(~c))@Moebius(alph,0,0,coAlph)
    
    def MatrixInverse(self):
        """Standard formula for 2x2 Quaternionic Matrix Inversion"""
        (a,b,c,d) = self.abcd()
        A=(a-b*(d.inverse())*c).inverse()
        B=(c-d*(b.inverse())*a).inverse()
        C=(b-a*(c.inverse())*d).inverse()
        D=(d-c*(a.inverse())*b).inverse()
        return Moebius(A,B,C,D)
    
    def __invert__(self):
        """Inverse Transformation (Not Matrix)"""
        (a,b,c,d) = self.abcd()
        O = 0*(a+b+c+d)
        (a,b,c,d) =  (O+x for x in (a,b,c,d)) #(secure(x,O) for x in (a,b,c,d)) #To avoid bugs for ill-defined `~`
        A = (~d)*d*(~a)-(~c)*d*(~b)
        B = (~b)*b*(~c)-(~a)*b*(~d)
        C = (~c)*c*(~b)-(~d)*c*(~a)
        D = (~a)*a*(~d)-(~b)*a*(~c)
        return Moebius(A,B,C,D)
        


    
    
    def __call__(self,q):
        (a,b,c,d) = self.abcd()
        return (a*q+b)/(c*q+d)
    
    def __matmul__(self,other):
        ((a,b,c,d),(A,B,C,D)) = (self.abcd(),other.abcd()) #(self.a,self.b,self.c,self.d,other.a,other.b,other.c,other.d)
        return Moebius(a*A+b*C, a*B+b*D, c*A+d*C, c*B+d*D)
    
    def __repr__(self):
        return "a: " + str(self.a) + ";\nb: " + str(self.b) + ";\nc: " + str(self.c) + ";\nd: " + str(self.d)+";"


def dualize(a,b):
    a.dual = b
    b.dual = a
