diff --git a/test_primitive.py b/test_primitive.py
index 2c0376d..fe3c2ae 100644
--- a/test_primitive.py
+++ b/test_primitive.py
@@ -1,89 +1,55 @@
-from operator import itemgetter
+from sage.all import *
+import cProfile, pstats
+import itertools
 load('aplK.py')
 
-s3 = SimplicialAlgebraplN(3)
-s3.inject_variables()
-q = y1*y2*y3
-p = q.propagate([0,1,2,3],4)
-p = p.reduce()
-s4 = SimplicialAlgebraplN(4)
-s4.inject_variables()
-#p = -t3*y0*y1*y2 + t2*y0*y1*y3 - t1*y0*y2*y3 + t0*y1*y2*y3
-#p = p.reduce()
+def test():
+    gens = load("K3_aplKgens")
+    S = simplicial_complexes.K3Surface()
+    S.set_immutable()
+    #S = SimplicialComplex([[0,1,3,4,5],[0,2,6,7,8],[1,2,9,10,11,12],[11,13],[12,13]],is_mutable=False)
+    A = AlgebraAplK(S)
+    H4 = A.cohomology_group(4)
+    #gens_dim2 = gens[2]
+    """
+    dim = 2
+    g2_0 = gens[dim][(dim,0)]
+    g2_1 = gens[dim][(dim,1)]
+    g2_2 = gens[dim][(dim,2)]
+    """
+    g4_0 = gens[4][(4,0)]
+    products = {}
+    elements = {}
+    for i in range(len(gens[2].keys())):
+        for j in range(i,len(gens[2].keys())):
+            g1 = gens[2][(2,i)]
+            g2 = gens[2][(2,j)]
+            g = g1*g2
+            c = g.integral()
+            v = A.cochain_to_vector(c,4)
+            products[(i,j)] = H4(v)
+            elements[(i,j)] = g-H4(v)[0]*g4_0
 
-## def _get_primitive_basis(self):
-basis = []
-for monomial in p.dict().keys():
-    basis += p._get_primitive_basis_of_a_monomial(monomial)
-basis = list(set(basis))
-terms = []
-apl = p.parent()
-for b in basis:
-    dim = len(b)//2 - 1
-    terms.append(p._get_term(1, b, apl.gens()))
-terms.sort(reverse=False)
-"""
-with open("output/primitive_basis_without_kernel:terms.txt", 'w') as fp:
-    print("original polynomial: " + str(p),file=fp)
-    print(terms,file=fp)
-"""           
-#print("primitive basis without kernel terms")
-#print(terms)
+    save(elements,"K3_gens_prods_deg4")
+    primitives = {}
 
-kernel_terms = []
-kernel_terms += p._add_kernel_terms_to_primitive_basis(basis)
-kernel_terms = list(set(kernel_terms))
-terms = []
-apl = p.parent()
-for b in kernel_terms:
-    dim = len(b)//2 - 1
-    terms.append(p._get_term(1, b, apl.gens()))
-terms.sort(reverse=False)
-
-#print("kernel terms")
-#print(terms)
-
-def get_kernel_terms(p, kernel_basis, parent, degree):
-
-    kernel_terms = []
-    for k in kernel_basis:
-        dim = len(k)//2 - 1
-        term = p._get_term(1, k, parent.gens())
-        e = k[5:]
-        gens_idx = []
-        for idx,ex in enumerate(e):
-            if ex:
-                gens_idx.append(idx)
-        list_ring_gens = [itemgetter(*gens_idx)(parent.gens())]*degree
-        sets = list_ring_gens+[[term]]
-        cp = cartesian_product(sets)
-        kernel_terms += [prod(e) for e in cp.list()]
-    return list(set(kernel_terms))
+    for key in elements.keys():
+        print("Solving primitive of product " + str(key), flush=True)
+        g = elements[key]
+        f = g.solve_primitive()
+        primitives[key] = f
+        save(primitives,"K3_primitives_prods_deg4")
         
-        
-"""
-basis = []
+    for key in elements.keys():
+        print("Checking differential of element " + str(key),flush=True)
+        g = elements[key]
+        f = primitives[key]
+        print("f.differential() == g", flush=True)
+        print(f.differential() == g, flush=True)
 
-for degree in range(5):
-    basis = get_kernel_terms(p, kernel_terms, s4, degree)    
-    with open("output/kernel_terms_degree_"+str(degree)+".txt", 'w') as fp:
-        for b in basis:
-            print(b,file=fp)
-#implicit fp.close() here as "with" scope ends.
-"""
-  
-import itertools
-def findsubsets(S,m):
-    return set(itertools.combinations(S, m))
+    print("Check finished", flush=True)
     
-apl = p.parent()
-N = len(apl.gens())/2
-ygens = apl.gens()[N+1:] #N+1 because we remove y0 generator
-degree = 3
-listProds = [prod(s) for s in findsubsets(ygens,degree)]
-V = VectorSpace(GF(degree),N-1)
-listCoefs = [tuple([0]+list(v)+N*[0]) for v in V]
-listPolinomials = [p._get_term(1, c, parent.gens()) for c in listCoefs]
-cp = cartesian_product([listProds,listPolinomials])
-elements = [prod(c) for c in cp]
-
+#############################################
+cProfile.runctx("test()", globals(), locals(), "Profile.prof")
+s = pstats.Stats("Profile.prof")
+s.strip_dirs().sort_stats("cumtime").print_stats()