Ticket #10771: trac10771-lcm_gcd_fraction_fields.patch

File trac10771-lcm_gcd_fraction_fields.patch, 19.7 KB (added by SimonKing, 9 years ago)

Implement gcd and lcm for general fields, with the special case of the fraction field of a unique factorization domain

  • sage/categories/fields.py 

    # HG changeset patch
# User Simon King <simon.king@uni-jena.de>
# Date 1297947839 -3600
# Node ID 6289126bdd0bdd30ed9b7467572af6bf5e80030b
# Parent  fb00ec75853019eb9799fd863b193fe82ee97c74
#10771: gcd/lcm for fraction fields of UFDs and for general fields.
In the case of fraction fields, it restricts to gcd/lcm in the base ring.
For general fields, properly restrict to ZZ, but otherwise return 0/1.

diff --git a/sage/categories/fields.py b/sage/categories/fields.py
    a b  
    137137            return True
    138138
    139139    class ElementMethods:
    140         pass
     140        # Fields are unique factorization domains, so, there is gcd and lcm
     141        # Of course, in general gcd and lcm in a field are not very interesting.
     142        # However, they should be implemented!
     143        def gcd(self,other):
     144            """
     145            Greatest common divisor.
     146
     147            NOTE:
     148
     149            Since we are in a field and the greatest common divisor is
     150            only determined up to a unit, it is correct to either return
     151            zero or one. Note that fraction fields of unique factorization
     152            domains provide a more sophisticated gcd.
     153
     154            EXAMPLES::
     155
     156                sage: GF(5)(1).gcd(GF(5)(1))
     157                1
     158                sage: GF(5)(1).gcd(GF(5)(0))
     159                1
     160                sage: GF(5)(0).gcd(GF(5)(0))
     161                0
     162
     163            For fields of characteristic zero (i.e., containing the
     164            integers as a sub-ring), evaluation in the integer ring is
     165            attempted. This is for backwards compatibility::
     166
     167                sage: gcd(6.0,8); gcd(6.0,8).parent()
     168                2
     169                Integer Ring
     170
     171            If this fails, we resort to the default we see above::
     172
     173                sage: gcd(6.0*CC.0,8*CC.0); gcd(6.0*CC.0,8*CC.0).parent()
     174                1.00000000000000
     175                Complex Field with 53 bits of precision
     176
     177            AUTHOR:
     178
     179            - Simon King (2011-02): Trac ticket #10771
     180
     181            """
     182            P = self.parent()
     183            try:
     184                other = P(other)
     185            except (TypeError, ValueError):
     186                raise ArithmeticError, "The second argument can not be interpreted in the parent of the first argument. Can't compute the gcd"
     187            from sage.rings.integer_ring import ZZ
     188            if ZZ.is_subring(P):
     189                try:
     190                    return ZZ(self).gcd(ZZ(other))
     191                except TypeError:
     192                    pass
     193            # there is no custom gcd, so, we resort to something that always exists
     194            # (that's new behaviour)
     195            if self==0 and other==0:
     196                return P.zero()
     197            return P.one()           
     198
     199        def lcm(self,other):
     200            """
     201            Least common multiple.
     202
     203            NOTE:
     204
     205            Since we are in a field and the least common multiple is
     206            only determined up to a unit, it is correct to either return
     207            zero or one. Note that fraction fields of unique factorization
     208            domains provide a more sophisticated lcm.
     209
     210            EXAMPLES::
     211
     212                sage: GF(2)(1).lcm(GF(2)(0))
     213                0
     214                sage: GF(2)(1).lcm(GF(2)(1))
     215                1
     216
     217            If the field contains the integer ring, it is first
     218            attempted to compute the gcd there::
     219
     220                sage: lcm(15.0,12.0); lcm(15.0,12.0).parent()
     221                60
     222                Integer Ring
     223
     224            If this fails, we resort to the default we see above::
     225
     226                sage: lcm(6.0*CC.0,8*CC.0); lcm(6.0*CC.0,8*CC.0).parent()
     227                1.00000000000000
     228                Complex Field with 53 bits of precision
     229                sage: lcm(15.2,12.0)
     230                1.00000000000000
     231
     232            AUTHOR:
     233
     234            - Simon King (2011-02): Trac ticket #10771
     235
     236            """
     237            P = self.parent()
     238            try:
     239                other = P(other)
     240            except (TypeError, ValueError):
     241                raise ArithmeticError, "The second argument can not be interpreted in the parent of the first argument. Can't compute the lcm"
     242            from sage.rings.integer_ring import ZZ
     243            if ZZ.is_subring(P):
     244                try:
     245                    return ZZ(self).lcm(ZZ(other))
     246                except TypeError:
     247                    pass
     248            # there is no custom lcm, so, we resort to something that always exists
     249            if self==0 or other==0:
     250                return P.zero()
     251            return P.one()

  • sage/categories/quotient_fields.py 

    diff --git a/sage/categories/quotient_fields.py b/sage/categories/quotient_fields.py
    a b  
    5252        def denominator(self):
    5353            pass
    5454
     55        def gcd(self,other):
     56            """
     57            Greatest common divisor
     58
     59            NOTE:
     60
     61            In a field, the greatest common divisor is not very
     62            informative, as it is only determined up to a unit. But in
     63            the fraction field of an integral domain that provides
     64            both gcd and lcm, it is possible to be a bit more specific
     65            and define the gcd uniquely up to a unit of the base ring
     66            (rather than in the fraction field).
     67
     68            AUTHOR:
     69
     70            - Simon King (2011-02): See trac ticket #10771
     71
     72            EXAMPLES::
     73
     74                sage: R.<x>=QQ[]
     75                sage: p = (1+x)^3*(1+2*x^2)/(1-x^5)
     76                sage: q = (1+x)^2*(1+3*x^2)/(1-x^4)
     77                sage: factor(p)
     78                (-2) * (x - 1)^-1 * (x + 1)^3 * (x^2 + 1/2) * (x^4 + x^3 + x^2 + x + 1)^-1
     79                sage: factor(q)
     80                (-3) * (x - 1)^-1 * (x + 1) * (x^2 + 1)^-1 * (x^2 + 1/3)
     81                sage: gcd(p,q)
     82                (x + 1)/(x^7 + x^5 - x^2 - 1)
     83                sage: factor(gcd(p,q))
     84                (x - 1)^-1 * (x + 1) * (x^2 + 1)^-1 * (x^4 + x^3 + x^2 + x + 1)^-1
     85                sage: factor(gcd(p,1+x))
     86                (x - 1)^-1 * (x + 1) * (x^4 + x^3 + x^2 + x + 1)^-1
     87                sage: factor(gcd(1+x,q))
     88                (x - 1)^-1 * (x + 1) * (x^2 + 1)^-1
     89
     90            TESTS:
     91
     92            The following tests that the fraction field returns a correct gcd
     93            even if the base ring does not provide lcm and gcd::
     94
     95                sage: R = ZZ.extension(x^2+5,names='q'); R
     96                Order in Number Field in q with defining polynomial x^2 + 5
     97                sage: R.1
     98                q
     99                sage: gcd(R.1,R.1)
     100                Traceback (most recent call last):
     101                ...
     102                TypeError: unable to find gcd of q and q
     103                sage: (R.1/1).parent()
     104                Number Field in q with defining polynomial x^2 + 5
     105                sage: gcd(R.1/1,R.1)
     106                1
     107                sage: gcd(R.1/1,0)
     108                1
     109                sage: gcd(R.zero(),0)
     110                0           
     111
     112            """
     113            try:
     114                other = self.parent()(other)
     115            except (TypeError, ValueError):
     116                raise ArithmeticError, "The second argument can not be interpreted in the parent of the first argument. Can't compute the gcd"
     117            try:
     118                selfN = self.numerator()
     119                selfD = self.denominator()
     120                selfGCD = selfN.gcd(selfD)
     121                otherN = other.numerator()
     122                otherD = other.denominator()
     123                otherGCD = otherN.gcd(otherD)
     124                selfN = selfN // selfGCD
     125                selfD = selfD // selfGCD
     126                otherN = otherN // otherGCD
     127                otherD = otherD // otherGCD
     128                return selfN.gcd(otherN)/selfD.lcm(otherD)
     129            except (AttributeError, NotImplementedError, TypeError, ValueError):
     130                if self==0 and other==0:
     131                    return self.parent().zero()
     132                return self.parent().one()
     133
     134        def lcm(self,other):
     135            """
     136            Least common multiple
     137
     138            NOTE:
     139
     140            In a field, the least common multiple is not very
     141            informative, as it is only determined up to a unit. But in
     142            the fraction field of an integral domain that provides
     143            both gcd and lcm, it is reasonable to be a bit more
     144            specific and to define the least common multiple so that
     145            it restricts to the usual least common multiple in the
     146            base ring and is unique up to a unit of the base ring
     147            (rather than up to a unit of the fraction field).
     148
     149            AUTHOR:
     150
     151            - Simon King (2011-02): See trac ticket #10771
     152
     153            EXAMPLES::
     154
     155                sage: R.<x>=QQ[]
     156                sage: p = (1+x)^3*(1+2*x^2)/(1-x^5)
     157                sage: q = (1+x)^2*(1+3*x^2)/(1-x^4)
     158                sage: factor(p)
     159                (-2) * (x - 1)^-1 * (x + 1)^3 * (x^2 + 1/2) * (x^4 + x^3 + x^2 + x + 1)^-1
     160                sage: factor(q)
     161                (-3) * (x - 1)^-1 * (x + 1) * (x^2 + 1)^-1 * (x^2 + 1/3)
     162                sage: factor(lcm(p,q))
     163                (x - 1)^-1 * (x + 1)^3 * (x^2 + 1/3) * (x^2 + 1/2)
     164                sage: factor(lcm(p,1+x))
     165                (x + 1)^3 * (x^2 + 1/2)
     166                sage: factor(lcm(1+x,q))
     167                (x + 1) * (x^2 + 1/3)
     168
     169            TESTS:
     170
     171            The following tests that the fraction field returns a correct lcm
     172            even if the base ring does not provide lcm and gcd::
     173
     174                sage: R = ZZ.extension(x^2+5,names='q'); R
     175                Order in Number Field in q with defining polynomial x^2 + 5
     176                sage: R.1
     177                q
     178                sage: lcm(R.1,R.1)
     179                Traceback (most recent call last):
     180                ...
     181                TypeError: unable to find lcm of q and q
     182                sage: (R.1/1).parent()
     183                Number Field in q with defining polynomial x^2 + 5
     184                sage: lcm(R.1/1,R.1)
     185                1
     186                sage: lcm(R.1/1,0)
     187                0
     188                sage: lcm(R.zero(),0)
     189                0           
     190
     191            """
     192            try:
     193                other = self.parent()(other)
     194            except (TypeError, ValueError):
     195                raise ArithmeticError, "The second argument can not be interpreted in the parent of the first argument. Can't compute the lcm"
     196            try:
     197                selfN = self.numerator()
     198                selfD = self.denominator()
     199                selfGCD = selfN.gcd(selfD)
     200                otherN = other.numerator()
     201                otherD = other.denominator()
     202                otherGCD = otherN.gcd(otherD)
     203                selfN = selfN // selfGCD
     204                selfD = selfD // selfGCD
     205                otherN = otherN // otherGCD
     206                otherD = otherD // otherGCD
     207                return selfN.lcm(otherN)/selfD.gcd(otherD)
     208            except (AttributeError, NotImplementedError, TypeError, ValueError):
     209                if self==0 or other==0:
     210                    return self.parent().zero()
     211                return self.parent().one()
     212
    55213        def factor(self, *args, **kwds):
    56214            """
    57215            Return the factorization of ``self`` over the base ring.
    … …  
    214372            Returns the derivative of this rational function with respect to the
    215373            variable ``var``.
    216374           
    217             Over an ring with a working GCD implementation, the derivative of a
     375            Over an ring with a working gcd implementation, the derivative of a
    218376            fraction `f/g`, supposed to be given in lowest terms, is computed as
    219377            `(f'(g/d) - f(g'/d))/(g(g'/d))`, where `d` is a greatest common
    220378            divisor of `f` and `g`.

  • sage/rings/arith.py 

    diff --git a/sage/rings/arith.py b/sage/rings/arith.py
    a b  
    14121412   
    14131413    Additional keyword arguments are passed to the respectively called
    14141414    methods.
     1415
     1416    OUTPUT:
     1417
     1418    The given elements are first coerced into a common parent. Then,
     1419    their greatest common divisor *in that common parent* is returned.
    14151420   
    14161421    EXAMPLES::
    14171422   
    … …  
    14191424        1
    14201425        sage: GCD(97*10^15, 19^20*97^2)
    14211426        97
    1422         sage: GCD(2/3, 4/3)
    1423         1
     1427        sage: GCD(2/3, 4/5)
     1428        2/15
    14241429        sage: GCD([2,4,6,8])
    14251430        2
    14261431        sage: GCD(srange(0,10000,10))  # fast  !!
    … …  
    14511456        0
    14521457        sage: type(gcd([]))
    14531458        <type 'sage.rings.integer.Integer'>
     1459
     1460    TESTS:
     1461
     1462    The following shows that indeed coercion takes place before computing
     1463    the gcd. This behaviour was introduced in trac ticket #10771::
     1464
     1465        sage: R.<x>=QQ[]
     1466        sage: S.<x>=ZZ[]
     1467        sage: p = S.random_element()
     1468        sage: q = R.random_element()
     1469        sage: parent(gcd(1/p,q))
     1470        Fraction Field of Univariate Polynomial Ring in x over Rational Field
     1471        sage: parent(gcd([1/p,q]))
     1472        Fraction Field of Univariate Polynomial Ring in x over Rational Field
     1473
     1474
    14541475    """
    14551476    # Most common use case first:
    14561477    if b is not None:
    1457         if hasattr(a, "gcd"):
    1458             return a.gcd(b, **kwargs)
    1459         else:
     1478        from sage.structure.element import get_coercion_model
     1479        cm = get_coercion_model()
     1480        a,b = cm.canonical_coercion(a,b)
     1481        try:
     1482            GCD = a.gcd
     1483        except AttributeError:
    14601484            try:
    14611485                return ZZ(a).gcd(ZZ(b))
    14621486            except TypeError:
    14631487                raise TypeError, "unable to find gcd of %s and %s"%(a,b)
    1464    
     1488        return GCD(b, **kwargs)
     1489
    14651490    from sage.structure.sequence import Sequence
    14661491    seq = Sequence(a)
    14671492    U = seq.universe()
    … …  
    15281553   
    15291554    -  ``a`` - a list or tuple of elements of a ring with
    15301555       lcm
    1531    
     1556
     1557    OUTPUT:
     1558
     1559    First, the given elements are coerced into a common parent. Then,
     1560    their least common multiple *in that parent* is returned.
    15321561   
    15331562    EXAMPLES::
    15341563   
    … …  
    15451574        sage: v = LCM(range(1,10000))   # *very* fast!
    15461575        sage: len(str(v))
    15471576        4349
     1577
     1578    TESTS:
     1579
     1580    The following tests against a bug that was fixed in trac
     1581    ticket #10771::
     1582
     1583        sage: lcm(4/1,2)
     1584        4
     1585
     1586    The following shows that indeed coercion takes place before
     1587    computing the least common multiple::
     1588
     1589        sage: R.<x>=QQ[]
     1590        sage: S.<x>=ZZ[]
     1591        sage: p = S.random_element()
     1592        sage: q = R.random_element()
     1593        sage: parent(lcm([1/p,q]))
     1594        Fraction Field of Univariate Polynomial Ring in x over Rational Field
     1595
    15481596    """
    15491597    # Most common use case first:
    15501598    if b is not None:
    1551         if hasattr(a, "lcm"):
    1552             return a.lcm(b)
    1553         else:
     1599        from sage.structure.element import get_coercion_model
     1600        cm = get_coercion_model()
     1601        a,b = cm.canonical_coercion(a,b)
     1602        try:
     1603            LCM = a.lcm
     1604        except AttributeError:
    15541605            try:
    15551606                return ZZ(a).lcm(ZZ(b))
    15561607            except TypeError:
    1557                 raise TypeError, "unable to find gcd of %s and %s"%(a,b)
    1558    
     1608                raise TypeError, "unable to find lcm of %s and %s"%(a,b)
     1609        return LCM(b)
     1610
    15591611    from sage.structure.sequence import Sequence
    15601612    seq = Sequence(a)
    15611613    U = seq.universe()

  • sage/rings/polynomial/multi_polynomial.pyx 

    diff --git a/sage/rings/polynomial/multi_polynomial.pyx b/sage/rings/polynomial/multi_polynomial.pyx
    a b  
    918918            sage: f.content().parent()
    919919            Integer Ring
    920920
    921         TESTS::
     921        TESTS:
     922
     923        Since trac ticket #10771, the gcd in QQ restricts to the
     924        gcd in ZZ.
    922925
    923926            sage: R.<x,y> = QQ[]
    924927            sage: f = 4*x+6*y
    925             sage: f.content()
    926             1
     928            sage: f.content(); f.content().parent()
     929            2
     930            Rational Field
    927931
    928932        """
    929933        from sage.rings.arith import gcd
    930934        from sage.rings.all import ZZ
    931 
    932         return gcd(self.coefficients(),integer=self.parent() is ZZ)
     935        return gcd(self.coefficients())
    933936
    934937    def is_generator(self):
    935938        r"""

  • sage/rings/rational.pyx 

    diff --git a/sage/rings/rational.pyx b/sage/rings/rational.pyx
    a b  
    896896        from sage.rings.arith import gcd, lcm
    897897        return gcd(nums) / lcm(denoms)
    898898       
    899     def gcd(self, other, **kwds):
    900         """
    901         Return a gcd of the rational numbers self and other.
    902 
    903         If self = other = 0, this is by convention 0.  In all other
    904         cases it can (mathematically) be any nonzero rational number,
    905         but for simplicity we choose to always return 1.
    906        
    907         EXAMPLES::
    908 
    909             sage: gcd(1/3, 2/1)
    910             1
    911             sage: gcd(1/1, 0/1)
    912             1
    913             sage: gcd(0/1, 0/1)
    914             0
    915         """
    916         if self == 0 and other == 0:
    917             return Rational(0)
    918         else:
    919             return Rational(1)
     899#    def gcd_rational(self, other, **kwds):
     900#        """
     901#        Return a gcd of the rational numbers self and other.
     902#
     903#        If self = other = 0, this is by convention 0.  In all other
     904#        cases it can (mathematically) be any nonzero rational number,
     905#        but for simplicity we choose to always return 1.
     906#       
     907#        EXAMPLES::
     908#
     909#            sage: (1/3).gcd_rational(2/1)
     910#            1
     911#            sage: (1/1).gcd_rational(0/1)
     912#            1
     913#            sage: (0/1).gcd_rational(0/1)
     914#            0
     915#        """
     916#        if self == 0 and other == 0:
     917#            return Rational(0)
     918#        else:
     919#            return Rational(1)
    920920
    921921    def valuation(self, p):
    922922        r"""

  • sage/rings/ring.pyx 

    diff --git a/sage/rings/ring.pyx b/sage/rings/ring.pyx
    a b  
    16121612            2^3 * 11
    16131613
    16141614        In a field, any nonzero element is a GCD of any nonempty set
    1615         of elements.  For concreteness, Sage returns 1 in these cases::
     1615        of nonzero elements. In previous versions, Sage used to return
     1616        1 in the case of the rational field. However, since trac
     1617        ticket #10771, the rational field is considered as the
     1618        *fraction field* of the integer ring. For the fraction field
     1619        of an integral domain that provides both GCD and LCM, it is
     1620        possible to pick a GCD that is compatible with the GCD of the
     1621        base ring::
    16161622
    16171623            sage: QQ.gcd(ZZ(42), ZZ(48)); type(QQ.gcd(ZZ(42), ZZ(48)))
    1618             1
     1624            6
    16191625            <type 'sage.rings.rational.Rational'>
    16201626            sage: QQ.gcd(1/2, 1/3)
    1621             1
     1627            1/6
    16221628
    16231629        Polynomial rings over fields are GCD domains as well. Here is a simple
    16241630        example over the ring of polynomials over the rationals as well as

  • sage/stats/basic_stats.py 

    diff --git a/sage/stats/basic_stats.py b/sage/stats/basic_stats.py
    a b  
    165165        sage: std([])
    166166        NaN
    167167        sage: std([I, sqrt(2), 3/5])
    168         sqrt(1/450*(5*sqrt(2) + 5*I - 6)^2 + 1/450*(5*sqrt(2) - 10*I + 3)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2)
     168        sqrt(1/450*(-5*sqrt(2) + 10*I - 3)^2 + 1/450*(-5*sqrt(2) - 5*I + 6)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2)
    169169        sage: std([RIF(1.0103, 1.0103), RIF(2)])
    170170        0.6998235813403261?
    171171        sage: import numpy
    … …  
    230230        sage: variance([])                     
    231231        NaN
    232232        sage: variance([I, sqrt(2), 3/5])       
    233         1/450*(5*sqrt(2) + 5*I - 6)^2 + 1/450*(5*sqrt(2) - 10*I + 3)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2
     233        1/450*(-5*sqrt(2) + 10*I - 3)^2 + 1/450*(-5*sqrt(2) - 5*I + 6)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2
    234234        sage: variance([RIF(1.0103, 1.0103), RIF(2)])
    235235        0.4897530450000000?
    236236        sage: import numpy