Ticket #11725: trac_11725-random-algebraic-numbers.patch

File trac_11725-random-algebraic-numbers.patch, 5.3 KB (added by rbeezer, 10 years ago)
  • sage/rings/qqbar.py

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1314143349 25200
    # Node ID e341c0e6d80313deff7930d714d411004ebb7d25
    # Parent  2a2abbcad325ccca9399981ceddf5897eb467e64
    11725: random elements of QQbar
    
    diff -r 2a2abbcad325 -r e341c0e6d803 sage/rings/qqbar.py
    a b  
    938938        """
    939939        return AlgebraicNumber(ANRoot(poly, interval, multiplicity))
    940940
     941    def random_element(self, poly_degree=2, *args, **kwds):
     942        r"""
     943        Returns a random algebraic number.
     944       
     945        INPUT:
     946       
     947        - ``poly_degree`` - default: 2 - degree of the random
     948          polynomial over the integers that the algebraic number
     949          is a root of.  This is not the degree of the minimal
     950          polynomial of the number.  Increase this parameter to
     951          achieve a greater diversity of algebraic numbers, at a
     952          cost of greater compuation time.  You can also vary the
     953          distribution of the coefficients but that will not vary
     954          the degree of the extension containing the element.
     955         
     956        - ``args``, ``kwds`` - arguments and keywords passed to
     957          the random number generator for elements of ``ZZ``, the
     958          integers.  See
     959          :meth:`~sage.rings.integer_ring.IntegerRing_class.random_element`
     960          for details, or see example below.
     961         
     962        OUTPUT:
     963       
     964        An element of ``QQbar``, the field of algebraic
     965        numbers (see :mod:`sage.rings.qqbar`).
     966       
     967        ALGORITHM:
     968       
     969        A polynomial with degree between 1 and ``poly_degree``,
     970        with random integer coefficents is created.  A root of this
     971        polynomial is chosen at random.  The default degree is
     972        2 and the integer coefficents come from a distribution
     973        heavily weighted towards $0, \pm 1, \pm 2$.
     974       
     975        EXAMPLES::
     976       
     977            sage: a = QQbar.random_element()
     978            sage: a                         # random
     979            0.2626138748742799? + 0.8769062830975992?*I
     980            sage: a in QQbar
     981            True
     982       
     983            sage: b = QQbar.random_element(poly_degree=20)
     984            sage: b                         # random
     985            -0.8642649077479498? - 0.5995098147478391?*I
     986            sage: b in QQbar
     987            True
     988       
     989        Parameters for the distribution of the integer coefficients
     990        of the polynomials can be passed on to the random element method
     991        for integers.  For example if we do not include zero as a possible
     992        coefficient, there will never be a zero constant term, and
     993        thus never a zero root. ::
     994       
     995            sage: z = [QQbar.random_element(x=1, y=10) for _ in range(20)]
     996            sage: QQbar(0) in z
     997            False
     998           
     999         If you just want real algebraic numbers you can filter them out. 
     1000         Using an odd degree for the polynomials will insure some degree of
     1001         success.  ::
     1002         
     1003            sage: r = []
     1004            sage: while len(r) < 3:
     1005            ...     x = QQbar.random_element(poly_degree=3)
     1006            ...     if x in AA:
     1007            ...       r.append(x)
     1008            sage: (len(r) == 3) and all([z in AA for z in r])
     1009            True
     1010           
     1011        TESTS:
     1012       
     1013            sage: QQbar.random_element('junk')
     1014            Traceback (most recent call last):
     1015            ...
     1016            TypeError: polynomial degree must be an integer, not junk
     1017            sage: QQbar.random_element(poly_degree=0)
     1018            Traceback (most recent call last):
     1019            ...
     1020            ValueError: polynomial degree must be greater than zero, not 0
     1021           
     1022        Random vectors already have a 'degree' keyword, so
     1023        we cannot use that for the polynomial's degree.  ::
     1024       
     1025            sage: v = random_vector(QQbar, degree=2, poly_degree=3)
     1026            sage: v                                 # random
     1027            (0.4694381338921299?, -0.500000000000000? + 0.866025403784439?*I)
     1028        """
     1029        import sage.rings.all
     1030        import sage.misc.prandom
     1031        try:
     1032            poly_degree = sage.rings.all.ZZ(poly_degree)
     1033        except TypeError:
     1034            msg = "polynomial degree must be an integer, not {0}"
     1035            raise TypeError(msg.format(poly_degree))
     1036        if poly_degree < 1:
     1037            msg = "polynomial degree must be greater than zero, not {0}"
     1038            raise ValueError(msg.format(poly_degree))
     1039        R = PolynomialRing(sage.rings.all.ZZ, 'x')
     1040        p = R.random_element(degree=poly_degree, *args, **kwds)
     1041        # degree zero polynomials have no roots
     1042        # totally zero poly has degree -1
     1043        # add a random leading term
     1044        if p.degree() < 1:   
     1045            g = R.gen(0)
     1046            m = sage.misc.prandom.randint(1, poly_degree)
     1047            p = p + g**m
     1048        roots = p.roots(ring=QQbar, multiplicities=False)
     1049       
     1050        if len(roots) == 0:
     1051            print "FUBAR", p
     1052        # pick one at random
     1053        # could we instead just compute one root "randomly"?
     1054        m = sage.misc.prandom.randint(0, len(roots)-1)
     1055        return roots[m]
     1056
    9411057def is_AlgebraicField(F):
    9421058    return isinstance(F, AlgebraicField)
    9431059