Changeset 8201:4e9b1bf4cf99

Timestamp:

11/28/07 08:24:54 (5 years ago)

Author:

dfdeshom@…

Branch:

default

Message:

Better random mpoly generation as proposed by malb and Steffen
<sreidt@…>

File:

• sage/rings/polynomial/multi_polynomial_ring_generic.pyx (modified) (3 diffs)

sage/rings/polynomial/multi_polynomial_ring_generic.pyx

@@ -1 +1 @@
 include '../../ext/stdsage.pxi'
+
 
 from sage.structure.parent_gens cimport ParentWithGens
 import sage.misc.latex
-import multi_polynomial_ideal
+import multi_polynomial_ideal 
 from term_order import TermOrder
 from sage.rings.integer_ring import ZZ
@@ -9 +10 @@
 import multi_polynomial_element
 import polynomial_ring
+
 
 def is_MPolynomialRing(x):
@@ -423 +425 @@
         return unpickle_MPolynomialRing_generic_v1,(base_ring, n, names, order)
 
-
-    def random_element(self, degree=2, terms=5, *args, **kwds):
-        r"""
-        Return a random polynomial in this polynomial ring.
-        
+    def _precomp_counts(self,n, d):
+        """
+        Given a number of variable n and a degree d return a tuple (C,t)
+        such that C is a list of the cardinalities of the sets of
+        monomials up to degree d (including) in n variables and t is the
+        sum of these cardinalities.
+
         INPUT:
-            degree -- maximum total degree of resulting polynomial
-            terms  -- maximum number of terms to generate
-
-        OUTPUT: a random polynomial of total degree \code{degree}
-                and with \code{term} terms in it.
+            n -- number of variables
+            d -- degree
+
+        EXAMPLE:
+        sage: C,t = _precomp_counts(10,2)
+        sage: C[0]
+        1
+        sage: C[1]
+        10
+        sage: C[2]
+        55
+        sage: t
+        66
+
+        TESTS:
+        sage: C,t = _precomp_counts(1,2)
+        sage: C[0]
+        1
+        sage: C[1]
+        1
+        sage: C[2]
+        1
+        sage: t
+        3
+        
+        """
+        from sage.rings.arith import binomial
+        C = [1]  #d = 0
+        for dbar in xrange(1, d+1):
+            C.append(binomial(n+dbar-1, dbar))
+        t = sum(C)
+        return C, t
+            
+    def _to_monomial(self,i, n, d):
+        """
+        Given an index i, a number of variables n and a degree d return
+        the i-th monomial of degree d in n variables.
+        
+        INPUT:
+            i -- index: 0 <= i < binom(n+d-1,n-1)
+            n -- number of variables
+            d -- degree
+
+        EXAMPLE:
+        sage: _to_monomial(0,10,2)
+        (0, 0, 0, 0, 0, 0, 0, 0, 0, 2)
+        sage: _to_monomial(8,10,2)
+        (0, 0, 0, 0, 0, 0, 1, 1, 0, 0)
+        sage: _to_monomial(54,10,2)
+        (2, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+
+        NOTE: We do not check if the provided index/rank is within the
+        allowed range. If it is not an infinite loop will occure.
+        """
+        from sage.combinat import choose_nk
+        comb = choose_nk.from_rank(i, n+d-1, n-1)
+        if comb == []:
+            return (d,)
+        monomial = [ comb[0] ]
+        res = []
+        for j in range(n-2):
+            res.append(comb[j+1]-comb[j]-1)
+        monomial += res
+        monomial.append( n+d-1-comb[-1]-1 )
+        return tuple(monomial)
+
+    def _random_monomial_upto_degree_class(self,n, degree, counts=None, total=None):
+        """
+        Choose a random exponent tuple for $n$ variables with a random
+        degree $d$, i.e. choose the degree uniformly at random first
+        before choosing a random monomial.
+        
+        INPUT:
+            n -- number of variables
+            degree -- degree of monomials
+            counts -- ignored
+            total -- ignored
+            """
+        # bug: doesn't handle n=1
+        from sage.rings.arith import binomial
+        #Select random degree
+        d = ZZ.random_element(0,degree+1)
+        total = binomial(n+d-1, d)
+        
+        #Select random monomial of degree d
+        random_index = ZZ.random_element(0, total-1)
+        #Generate the corresponding monomial
+        return self._to_monomial(random_index, n, d)
+
+    def _random_monomial_upto_degree_uniform(self,n, degree, counts=None, total=None):
+        """
+        Choose a random exponent tuple for $n$ variables with a random
+        degree up to $d$, i.e. choose a random monomial uniformly random
+        from all monomials up to degree $d$. This discriminates against
+        smaller degrees because there are more monomials of bigger
+        degrees.
+    
+        INPUT:
+            n -- number of variables
+            degree -- degree of monomials
+            counts -- ignored
+            total -- ignored
+            """
+        if counts is not None or total is not None:
+            counts, total = self._precomp_counts(n, degree)
+            
+        #Select a random one
+        random_index = ZZ.random_element(0, total-1)
+        #Figure out which degree it corresponds to
+        d = 0
+        while random_index >= counts[d]:
+            random_index -= counts[d]
+            d += 1
+        #Generate the corresponding monomial
+        return self._to_monomial(random_index, n, d)
+
+    def random_element(self, d, t, choose_degree=False,*args, **kwargs):
+        """
+        Return a random polynomial of at most degree $d$ and at most $t$
+        terms.
+
+        First monomials are chosen uniformly random from the set of all
+        possible monomials of degree up to $d$ (inclusive). This means
+        that it is more likely that a monomial of degree $d$ appears than
+        a monomial of degree $d-1$ because the former class is bigger.
+
+        Exactly $t$ \emph{distinct} monomials are chosen this way and each
+        one gets a random coefficient (possibly zero) from the base ring
+        assigned. 
+
+        The returned polynomial is the sum of this list of terms.
+
+        INPUT:
+            d -- maximal degree (likely to be reached)
+            t -- number of terms requested
+            choose_degree -- choose degrees of monomials randomly first rather
+                             than monomials uniformly random.
+            **kwargs -- passed to the random element generator of the base ring
 
         EXAMPLES:
-            sage: [QQ['x,y'].random_element() for _ in range(5)]
-            [-1/14*x*y + 1/2*x, x*y + x - y + 1, 3*x*y + x - 1/2, 1/3*x*y - 5*x + 1/2*y + 7/6, 2*x*y + 1/2*x + 1]
-            sage: R = MPolynomialRing(ZZ, 'x,y',2 );
-            sage: R.random_element(2)          # random
-            -1*x*y + x + 15*y - 2
-            sage: R.random_element(12)         # random
-            x^4*y^5 + x^3*y^5 + 6*x^2*y^2 - x^2
-            sage: R.random_element(12,3)       # random
-            -3*x^4*y^2 - x^5 - x^4*y
-            sage: R.random_element(3)          # random
-            2*y*z + 2*x + 2*y
-
-            sage: R.<x,y> = MPolynomialRing(RR)
-            sage: R.random_element(2)          # random
-            -0.645358174399450*x*y + 0.572655401740132*x + 0.197478565033010
-
-            sage: R.random_element(41)         # random
-            -4*x^6*y^4*z^4*a^6*b^3*c^6*d^5 + 1/2*x^4*y^3*z^5*a^4*c^5*d^6 - 5*x^3*z^3*a^6*b^4*c*d^5 + 10*x^2*y*z^5*a^4*b^2*c^3*d^4 - 5*x^3*y^5*z*b^2*c^5*d
-
-        AUTHOR:
-            -- didier deshommes
-        """
-        # General strategy:
-        # generate n-tuples of numbers with each element in the tuple
-        # not greater than  (degree/n) so that the degree 
-        # (ie, the sum of the elements in the tuple) does not exceed
-        # their total degree
-        
-        n = self.__ngens         # length of the n-tuple
-        max_deg = int(degree/n)  # max degree for each term
-        R = self.base_ring()
-        
-        # restrict exponents to positive integers only
-        exponents = [ tuple([ZZ.random_element(0,max_deg+1) for _ in range(n)])
-                       for _ in range(terms) ]
-        coeffs = []
-        for _ in range(terms):
-            c = R.random_element(*args,**kwds)
-            while not c:
-                c = R.random_element(*args,**kwds)
-            coeffs.append(c) # allow only nonzero coefficients
-
-        d = dict( zip(tuple(exponents), coeffs) ) 
-        return self(multi_polynomial_element.MPolynomial_polydict(self, PolyDict(d)))
-
+        sage: P.<x,y,z> = PolynomialRing(QQ)
+        sage: random_element(P, 2, 5)
+        1/2*y^2 + x*z - x + 2*y + 19/2*z
+
+        sage: random_element(P, 2, 5, choose_degree=True)
+        28*x*z + y*z - z^2 - 1/2*x - 44/39
+  
+        sage: random_element(P, 2, 15)
+        Traceback (most recent call last):
+        ...
+        Cannot compute polynomial with more terms than exist.
+
+        sage: random_element(P, 0, 1)
+        1
+
+        sage: random_element(P, 2, 0)
+        0
+
+        """
+        from sage.combinat.integer_vector import IntegerVectors
+        from sage.rings.arith import binomial
+        
+        k = self.base_ring()
+        n = self.ngens()
+
+        #total or d is 0. Just return
+        if 0 == total or 0 == d:
+            return tuple([0]*n)
+
+        counts, total = self._precomp_counts(n, d)
+        
+
+        if t < 0:
+            raise TypeError, "Cannot compute polynomial with a negative number of terms."
+
+        elif t < total/2:
+            # we choose random monomials if t < total/2 because then we
+            # expect the algorithm to be faster than generating all
+            # monomials and picking a random index from the list. if t ==
+            # total/2 we expect every second random monomial to be a
+            # double such that our runtime is doubled in the worst case.
+            M = set()
+            if not choose_degree:
+                while t:
+                    m = self._random_monomial_upto_degree_uniform(n, d, counts, total)
+                    if not m in M:
+                        M.add(m)
+                        t -= 1
+            else:
+                while t:
+                    m = self._random_monomial_upto_degree_class(n, d)
+                    if not m in M:
+                        M.add(m)
+                        t -= 1
+        elif t <= total:
+            # generate a list of all monomials and choose among them
+            if not choose_degree:
+                M = sum([list(IntegerVectors(_d,n)) for _d in xrange(d+1)],[])
+                for mi in xrange(total - t): # we throw away those we don't need
+                    M.pop( ZZ.random_element(0,len(M)-1) )
+                M = map(tuple, M)
+            else:
+                M = [list(IntegerVectors(_d,n)) for _d in xrange(d+1)]
+                Mbar = []
+                for mi in xrange(t): 
+                    d = ZZ.random_element(0,len(M)) #choose degree at random
+                    m = ZZ.random_element(0,len(M[d])) # choose monomial at random
+                    Mbar.append( M[d].pop(m) ) # remove and insert
+                    if len(M[d]) == 0:
+                        M.pop(d) # bookkeeping
+                M = map(tuple, Mbar)
+                
+        else:
+            raise TypeError, "Cannot compute polynomial with more terms than exist."
+
+        C = [k.random_element(*args,**kwargs) for _ in range(len(M))]
+    
+        return self(dict(zip(M,C)))
+        
     def new_ring(self, names=None, order=None):
         """