Index: sagePolynomialOperations.sage
===================================================================
--- sagePolynomialOperations.sage	(revision 107)
+++ sagePolynomialOperations.sage	(revision 108)
@@ -497,6 +497,112 @@
     return((knownMonomials, polynomialsList))
 # End spo_polynomial_to_proto_matrix_2
 
+def spo_polynomial_to_polynomials_list_3(p, alpha, N, columnsWidth=0):
+    """
+    From p, alpha, N build a list of polynomials...
+    TODO: more in depth rationale...
+    
+    Our goal is to introduce each monomial with the smallest coefficient.
+     
+    
+        
+    Parameters
+    ----------
+    p: the (bivariate) polynomial;
+    pRing: the ring over which p is defined;
+    alpha:
+    N:
+    columsWidth: if == 0, no information is displayed, otherwise data is 
+                 printed in colums of columnsWitdth width.
+    """
+    pRing = p.parent()
+    knownMonomials = []
+    polynomialsList = []
+    pVariables = p.variables()
+    iVariable = pVariables[0]
+    tVariable = pVariables[1]
+    polynomialAtPower = pRing(1)
+    currentPolynomial = pRing(1)
+    pIdegree = p.degree(iVariable)
+    pTdegree = p.degree(tVariable)
+    currentIdegree = currentPolynomial.degree(iVariable)
+    nAtAlpha = N^alpha
+    nAtPower = nAtAlpha
+    polExpStr = ""
+    # We work from p^0 * N^alpha to p^alpha * N^0
+    for pPower in xrange(0, alpha + 1):
+        # pPower == 0 is a special case. We introduce all the monomials in i 
+        # up to i^pIdegree.
+        if pPower == 0:
+            # Notice who iPower runs up to i^pIdegree.
+            for iPower in xrange(0, pIdegree + 1): 
+                # No t power is taken into account as we limit our selves to
+                # degree 1 in t and make no "t-shifts".
+                if columnsWidth != 0:
+                    polExpStr = spo_expression_as_string(iPower, 
+                                                         0, 
+                                                         0, 
+                                                         alpha)
+                    print "->", polExpStr
+                currentExpression = iVariable^iPower * nAtAlpha
+                # polynomialAtPower == 1 here. Next line should be commented
+                # out but it does not work! Some conversion problem?
+                currentPolynomial = pRing(currentExpression)
+                polynomialsList.append(currentPolynomial) 
+            # End for iPower.
+        else: # pPower > 0: (p^1..p^alpha)
+            # This where we introduce the p^pPower * N^(alpha-pPower)
+            # polynomial. This is also where the t^pPower monomials shows up for
+            # the first time. It app
+            if columnsWidth != 0:
+                polExpStr = spo_expression_as_string(0, 0, pPower, alpha-pPower)
+                print "->", polExpStr
+            currentPolynomial = polynomialAtPower * nAtPower
+            polynomialsList.append(currentPolynomial) 
+            # Exit when pPower == alpha
+            if pPower == alpha:
+                return((knownMonomials, polynomialsList))
+            # This is where the "i-shifts" take place. Mixed terms, i^k * t^l
+            # (that were introduced at a previous
+            # stage or are introduced now) are only shifted to already existing 
+            # ones with the notable exception of i^iPower * t^pPower, which
+            # must be manually introduced.
+            # p^pPower is "shifted" to higher degrees in i as far as 
+            # possible, up to of the degree in i of p^(pPower+1).
+            # These "pure" i^k monomials can only show up with i multiplications.
+            for iPower in xrange(1, pIdegree + 1):
+                # The i^iPower * t^pPower monomial. Notice the alpha exponent
+                # for N.
+                internalIpower = iPower
+                for tPower in xrange(pPower,0,-1):
+                    if columnsWidth != 0:
+                        polExpStr = spo_expression_as_string(internalIpower, \
+                                                             tPower,  \
+                                                             0, \
+                                                             alpha)
+                        print "->", polExpStr
+                    currentExpression = i^internalIpower * t^tPower * nAtAlpha
+                    currentPolynomial = pRing(currentExpression)
+                    polynomialsList.append(currentPolynomial) 
+                    internalIpower += pIdegree
+                # End for tPower
+                # The i^iPower * p^pPower * N^(alpha-pPower) i-shift.
+                if columnsWidth != 0:
+                    polExpStr = spo_expression_as_string(iPower, \
+                                                         0,      \
+                                                         pPower, \
+                                                         alpha-pPower)
+                    print "->", polExpStr
+                currentExpression = i^iPower * nAtPower
+                currentPolynomial = pRing(currentExpression) * polynomialAtPower
+                polynomialsList.append(currentPolynomial) 
+            # End for iPower
+        polynomialAtPower *= p  
+        nAtPower /= N
+    # End for pPower loop
+    return((knownMonomials, polynomialsList))
+# End spo_polynomial_to_proto_matrix_3
+
 def spo_proto_to_column_matrix(protoMatrixColumns, \
                                knownMonomialsList, \
                                boundVar1, \
