diff --git a/aplK.py b/aplK.py
new file mode 100644
index 0000000..be84089
--- /dev/null
+++ b/aplK.py
@@ -0,0 +1,383 @@
+reset()
+from sage.structure.parent import Parent
+from sage.structure.element import Element
+from sage.structure.unique_representation import UniqueRepresentation
+from sage.categories.rings import Rings
+from sage.structure.richcmp import rich_to_bool, op_EQ, op_NE, op_LT, op_LE, op_GT, op_GE
+
+#####################################################################
+class Apln(UniqueRepresentation):
+    
+  def __init__(self,n):
+      self.n = n
+      ts = ['t'+str(t) for t in range(n+1)]
+      ys = ['y'+str(t) for t in range(n+1)]
+      degrees = tuple((n+1)*[(0,2)])+tuple((n+1)*[(1,2)])
+      self.gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
+      t_relation = -1
+      y_relation = 0
+      #t_relation = sum(self.gca.gens()[:n])-1
+      #y_relation = sum(self.gca.gens()[n:])
+      for i in range(n+1):
+          t_relation += self.gca.gens()[i]
+          y_relation += self.gca.gens()[n+1+i]
+      self.gca_with_relations = self.gca.quotient(self.gca.ideal(t_relation,y_relation),check=False)
+      dfs = {self.gca.gens_dict()[ts[i]]:self.gca.gens_dict()[ys[i]] for i in range(n+1)}
+      self.cdga = self.gca.cdg_algebra(dfs)
+  
+  def _repr_(self):
+    return "Algebra Apln of dimension %s" % self.n
+    
+  def inject_variables(self):
+      self.cdga.inject_variables()
+      
+  def generators(self):
+      return self.cdga.gens()
+      
+  def reduce_with_relations(self,pol):
+      return self.gca_with_relations(pol)
+      
+  def homogeneousPart(self,pol,degree):
+      parts = [pol.homogeneous_parts()[k] for k in pol.homogeneous_parts() if list(k)[0]==degree]
+      hompart = 0
+      for p in parts:
+          hompart += p
+      return hompart
+
+  def face(self,i):
+      if(self.n==0):
+          return 0
+      domain = self.cdga
+      codomain = self.cdga
+      tImages = []
+      yImages = []
+      n = len(domain.gens())//2
+      for k in range(n):
+          if(k<i):
+              tImages.append(codomain.gens()[k])
+              yImages.append(codomain.gens()[n+k])
+          elif(k==i):
+              tImages.append(0)
+              yImages.append(0)
+          elif(k>i):
+              tImages.append(codomain.gens()[k-1])
+              yImages.append(codomain.gens()[n+k-1])
+      images = tImages+yImages
+      return domain.hom(codomain=codomain,im_gens=images)
+      
+  def _reduce_polynomial(self,p,n):
+    ## Input:
+    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
+    #  - n: dimension of the algebra (apl)_n where p is defined
+    ## Output:
+    #  - The polynomial after applying the corresponding relations to reduce the variables t_n, y_n
+    #  - For practical reasons, we change the nomenclature in this way: t_i ==> xs_i, y_i ==> zs_i
+    
+    xs = ['x'+str(x) for x in range(n+1)]
+    zs = ['z'+str(z) for z in range(n+1)]
+    R = PolynomialRing(QQ,names=xs+zs)
+    r_gens = R.gens()
+    gca_gens = self.gca.gens()
+    pol = R(0)
+    for monomial_exponents in list(p.dict().keys()):
+      monomial = R(1)
+      coef = p.dict()[monomial_exponents]
+      for i in range(n+1):
+        monomial *= r_gens[i]**monomial_exponents[i]
+        monomial *= r_gens[n+1+i]**monomial_exponents[self.n+1+i]
+      pol += coef*monomial
+    xn_relation = 1-sum(r_gens[:n])
+    zn_relation = -sum(r_gens[n+1:-1])
+    pol = pol.subs({r_gens[n]:xn_relation,r_gens[2*n+1]:zn_relation})
+    return pol
+    
+  def _convert_integrable_polynomial(self,p,n):
+    ## Input:
+    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
+    #  - n: dimension of the algebra (apl)_n where p is defined
+    ## Output:
+    #  - A polynomial in the ring Q<xs_{0},..,xs_{n-1}>.
+    #    (After reducing the variables t_{n},y_{n} by the corresponding relations, 
+    #     the algorithm removes the variables y_{0},..,y_{n-1}). 
+
+    xs = ['x'+str(x) for x in range(n)]
+    zs = ['z'+str(z) for z in range(n)]
+    R = PolynomialRing(QQ,names=xs+zs)
+    pol = R(self._reduce_polynomial(p,n))
+    r_gens = R.gens()
+    relations = {}
+    for i in range(n):
+      relations[r_gens[n+i]]=1
+    return pol.subs(relations)
+        
+  def integral(self,p,n):
+    ## Input:
+    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
+    #  - n: dimension of the algebra (apl)_n where p is defined
+    ## Output:
+    #  - the integral of the polynomial after reducing the variables t_n, y_n by the corresponding relations 
+    #
+    xs = ['x'+str(x) for x in range(n)]
+    R = PolynomialRing(QQ,names=xs)
+    r_gens = R.gens()
+    pol = R(self._convert_integrable_polynomial(p,n))
+    result = pol
+    for i in range(n):
+      result = result.integral(r_gens[n-1-i])
+      result = result.subs({r_gens[n-1-i]:1-sum(r_gens[:n-1-i])})
+    return R(((-1)**(n*(n-1)/2)))*result
+            
+#####################################################################
+class AplK(UniqueRepresentation, Parent):
+    
+  def __init__(self,K, category=None):
+    self.K = K
+    self.n = K.dimension()
+    Parent.__init__(self,category=category )
+   
+  def _repr_(self):
+    return "Algebra AplK over %s" % self.K
+    
+  def _element_constructor_(self,values):
+    return self.element_class(self,values)
+    
+  def dimension(self):
+    return self.K.dimension()
+    
+  def simplicialComplex(self):
+    return self.K
+    
+  def _n_cells_sorted(self,n):
+    return self.K._n_cells_sorted(n)
+    
+#####################################################################
+class AplKElement(Element):
+
+  def __init__(self, parent,values):
+    self.elements = values
+    self.apl = Apln(parent.dimension())
+    Element.__init__(self,parent)
+    self.check()
+    
+  def _repr_(self):
+    return str(self.elements)
+    
+  def dimension(self):
+    return self.parent().dimension()
+    
+  def nullElement(self):
+    C = self.__class__
+    values = {}
+    for k in self.elements.keys():
+      values[k]=0
+    return C(self.parent(),values)    
+    
+  def _add_(self,other):
+    C = self.__class__
+    if(self.elements.keys()!=other.elements.keys()):
+      raise ValueError("error: keys do not match")
+    newValues = {}
+    for k in self.elements.keys():
+      newValues[k] = self.elements[k]+other.elements[k]
+    return C(self.parent(),newValues)
+      
+  def _sub_(self,other):
+    C = self.__class__
+    if(self.elements.keys()!=other.elements.keys()):
+      raise ValueError("error: keys do not match")
+    newValues = {}
+    for k in self.elements.keys():
+      newValues[k] = self.elements[k]-other.elements[k]
+    return C(self.parent(),newValues)
+      
+  def _mul_(self,other):
+    C = self.__class__
+    if(self.elements.keys()!=other.elements.keys()):
+      raise ValueError("error: keys do not match")
+    newValues = {}
+    for k in self.elements.keys():
+      newValues[k] = self.elements[k]*other.elements[k]
+    return C(self.parent(),newValues)
+    
+  def _richcmp_(self, other, op):
+    if(op==op_EQ):
+      if(self.elements.keys()!=other.elements.keys()):
+        return False
+      for k in self.elements.keys():
+        if(self.elements[k]!=other.elements[k]):
+          return False
+      return True
+    elif(op==op_NE):
+      return not self._richcmp_(other,op_EQ)
+    else:
+      raise NotImplementedError("operation not implemented for this element class")
+      
+  def differential(self):
+    C = self.__class__
+    newValues = {}
+    for k in self.elements.keys():
+      newValues[k] = self.apl.cdga(self.elements[k]).differential()
+    return C(self.parent(),newValues)
+      
+  def check(self):
+    S = self.parent().simplicialComplex()
+    dim = S.dimension()
+    for n in range(1,dim+1):
+      print("Checking dimension " + str(n))
+      for simplex in S._n_cells_sorted(n):
+        for f in range(n+1):
+          sFace = tuple(S.face(simplex,f))
+          aFace = self.apl.face(f)(self.elements[tuple(simplex)])
+          if(self.elements[sFace] != aFace):
+            raise ValueError("the assignation is not valid")
+            
+  def homogeneous_component(self,degree):
+    C = self.__class__
+    component = {}
+    for k in self.elements.keys():
+      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
+    return C(self.parent(),component)
+    
+  def homogeneous_components(self):
+    components = {}
+    for n in range(self.dimension()+1):
+      components[n] = self.homogeneous_component(n)
+    return components
+               
+  def is_homogeneous(self):
+    isHomogeneous = True
+    for k in self.elements.keys():
+      isHomogeneous = isHomogeneous and self.apl.cdga(self.elements[k]).is_homogeneous()
+    return isHomogeneous
+  
+  def is_in_kernel(self):
+    diff = self.differential()
+    if(diff == self.nullElement()):
+      return True
+    else:
+      return False
+      
+  def _integrable_component(self):
+    C = self.__class__
+    component = {}
+    for k in self.elements.keys():
+      degree = len(k)-1
+      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
+      #component[k] = self.homogeneous_component(dim).elements[k]
+    return component
+    
+  def integral(self):
+    component = self._integrable_component()
+    cochain = {}
+    #integration of each polynomial
+    for k in component.keys():
+      print("Integrating simplex " + str(k))
+      n = len(k)
+      pol = self.apl.cdga(component[k])
+      cochain[k] = self.apl.integral(pol,n)
+    return cochain
+
+  
+      
+#####################################################################
+#S = SimplicialComplex([[0,1,2]],is_mutable=False)  
+S = SimplicialComplex([[0,1],[0,2],[1,2]],is_mutable=False)                                                                                 
+apl = Apln(S.dimension())
+apl.inject_variables()
+
+params1 = {}
+params1[(0,)] = t0
+params1[(1,)] = t0
+params1[(2,)] = 3*t0
+params1[(0,1)] = t0+t1
+params1[(0,2)] = t0+3*t1
+params1[(1,2)] = t0+3*t1
+
+params2 = {}
+params2[(0,)] = 3*t0
+params2[(1,)] = 4*t0
+params2[(2,)] = 5*t0
+params2[(0,1)] = 3*t0+4*t1
+params2[(0,2)] = 3*t0+5*t1
+params2[(1,2)] = 4*t0+5*t1
+
+params3 = {}
+for d in range(S.dimension()+1):
+  for c in S._n_cells_sorted(d):
+    params3[tuple(c)]=randint(1,100)*t0
+    
+params4 = {}
+params4[(0,)] = 3*t0+y0
+params4[(1,)] = 4*t0
+params4[(2,)] = 5*t0
+params4[(0,1)] = 3*t0+4*t1+y0
+params4[(0,2)] = 3*t0+5*t1+y0
+params4[(1,2)] = 4*t0+5*t1
+
+A = AplK(S)
+e1 = AplKElement(A,params1)
+e2 = AplKElement(A,params2)
+#e3 = AplKElement(A,params3)
+e4 = AplKElement(A,params4)
+
+
+class morphism_from_QQ_to_AplK(Morphism):
+
+  def __init__(self, codomain):
+    self.simplicialComplex= codomain.K
+    Morphism.__init__(self, QQ, codomain)
+
+  def _call_(self, x):
+    S = self.simplicialComplex
+    params1 = {}
+    for d in range(S.dimension()+1):
+      for c in S._n_cells_sorted(d):
+        params1[tuple(c)]=x
+    A = AplK(S)
+    return AplKElement(A,params1)
+    
+f = morphism_from_QQ_to_AplK(A)
+A.register_coercion(f)
+
+#p = t0*t1**2*y0+t0*y1
+
+
+S = SimplicialComplex([[0,1,2,3,4,5,6]],is_mutable=False)                                                                                 
+apl = Apln(S.dimension())
+apl.inject_variables()
+p = y0*y1*y2*y3*y4*y5+y0*y1*y2*y3*y6+t4*t5
+
+initial_assignation = {}
+initial_assignation[(0,1,2,3,4,5,6)] = p
+
+def generateElement(S,initial_assignation):
+  apl = Apln(S.dimension())
+  aplk_structure = initial_assignation.copy()
+  
+  for n in reversed(range(1,S.dimension()+1)):
+    print("Generating dimension "+ str(n))
+    for simplex in S._n_cells_sorted(n):
+      p = aplk_structure[tuple(simplex)]
+      for f in range(n+1):
+        sFace = tuple(S.face(simplex,f))
+        aFace = apl.face(f)(p)
+        if(sFace not in aplk_structure.keys()):
+          aplk_structure[tuple(sFace)] = aFace
+        elif(aplk_structure[tuple(sFace)] != aFace):
+          raise ValueError("Error during the checking. The assignation is not valid")
+  A = AplK(S)
+  return AplKElement(A,aplk_structure)
+          
+          
+g = generateElement(S,initial_assignation)
+#g.check()
+#cochain = g.integral()
+  
+#n = 1
+#a = Apln(n)
+#a.inject_variables()
+#p = t0*y1
+##print(p)
+#print(a._reduce_polynomial(p,n))
+#print(a._convert_integrable_polynomial(p,n))
+#print(a.integral(p,n))