Ticket #10596: trac-10596.patch

File trac-10596.patch, 42.7 KB (added by spancratz, 10 years ago)
  • sage/rings/arith.py

    diff -r b6a68db0aac8 sage/rings/arith.py
    a b  
    22712271        sage: odd_part(factorial(31))
    22722272        122529844256906551386796875
    22732273    """
    2274     n = ZZ(n)
    2275     return n._shift_helper(n.valuation(2), -1)
     2274    if not isinstance(n, integer.Integer):
     2275        n = ZZ(n)
     2276    return n.odd_part()
    22762277
    22772278def prime_to_m_part(n,m):
    22782279    """
  • sage/rings/integer.pyx

    diff -r b6a68db0aac8 sage/rings/integer.pyx
    a b  
    9393    sage: 'sage'*Integer(3)
    9494    'sagesagesage'
    9595
    96 COERCIONS: Returns version of this integer in the multi-precision
    97 floating real field R.
    98 
    99 ::
     96COERCIONS:
     97
     98Returns version of this integer in the multi-precision floating
     99real field R::
    100100
    101101    sage: n = 9390823
    102102    sage: RR = RealField(200)
    103103    sage: RR(n)
    104104    9.3908230000000000000000000000000000000000000000000000000000e6
    105105
    106 
    107106"""
    108107#*****************************************************************************
    109108#       Copyright (C) 2004,2006 William Stein <wstein@gmail.com>
     
    287286    cdef Integer temp
    288287    cdef int v_int
    289288    if power_index < 5:
    290         # It turns out that simple repeated division is very fast for relatively few digits.
    291         # I don't think this is a real algorithmic statement, it's an annoyance introduced by memory allocation.
    292         # I think that manual memory management with mpn_* would make the divide & conquer approach even faster, but the code would be much more complicated.
     289        # It turns out that simple repeated division is very fast for
     290        # relatively few digits.  I don't think this is a real algorithmic
     291        # statement, it's an annoyance introduced by memory allocation.
     292        # I think that manual memory management with mpn_* would make the
     293        # divide & conquer approach even faster, but the code would be much
     294        # more complicated.
    293295        _digits_naive(v,l,offset,power_list[0],digits)
    294296    else:
    295297        mpz_init(mpz_quot)
     
    434436        The class ``Integer`` is implemented in Cython, as a
    435437        wrapper of the GMP ``mpz_t`` integer type.
    436438   
    437    
    438439    EXAMPLES::
    439440   
    440441        sage: Integer(010)
     
    566567            sage: ZZ(MyFloat(5))
    567568            5
    568569        """
    569 
    570570        # TODO: All the code below should somehow be in an external
    571571        # cdef'd function.  Then e.g., if a matrix or vector or
    572572        # polynomial is getting filled by mpz_t's, it can use the
     
    760760        if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):       
    761761            return (<Integer>x)._xor(y)
    762762        return bin_op(x, y, operator.xor)
    763        
    764        
     763
    765764    def __richcmp__(left, right, int op):
    766765        """
    767766        cmp for integers
     
    814813        set_mpz(z,self.value)
    815814        return z
    816815
    817 
    818816    def list(self):
    819817        """
    820818        Return a list with this integer in it, to be compatible with the
     
    828826        """
    829827        return [ self ]
    830828
    831 
    832829    def  __dealloc__(self):
    833830        mpz_clear(self.value)
    834831
     
    915912        ::
    916913       
    917914            sage: big = 10^5000000
    918             sage: s = big.str()                 # long time (> 20 seconds)
    919             sage: len(s)                        # long time (depends on above defn of s)
     915            sage: s = big.str()       # long time (> 20 seconds)
     916            sage: len(s)              # long time (depends on above defn of s)
    920917            5000001
    921             sage: s[:10]                        # long time (depends on above defn of s)
     918            sage: s[:10]              # long time (depends on above defn of s)
    922919            '1000000000'
    923920        """
    924921        if base < 2 or base > 36:
     
    11261123       
    11271124        INPUT:
    11281125       
    1129        
    11301126        -  ``base`` - integer (default: 10)
    11311127       
    11321128        -  ``digits`` - optional indexable object as source for
     
    12731269                z = the_integer_ring._zero_element
    12741270            l = [z]*(s if s >= padto else padto)
    12751271            for i from 0<= i < s:
    1276                 if mpz_tstbit(self_abs.value,i):  # mpz_tstbit seems to return 0 for the high-order bit of negative numbers?!
     1272                # mpz_tstbit seems to return 0 for the high-order bit of
     1273                # negative numbers?!
     1274                if mpz_tstbit(self_abs.value,i):
    12771275                    l[i] = o
    12781276        else:
    12791277            s = mpz_sizeinbase(self.value, 2)
    12801278            if s > 256:
    12811279                _sig_on
    12821280
    1283             #   We use a divide and conquer approach (suggested
    1284             # by the prior author, malb?, of the digits method)
    1285             # here: for base b, compute b^2, b^4, b^8,
    1286             # ... (repeated squaring) until you get larger
    1287             # than your number; then compute (n // b^256,
    1288             # n % b ^ 256) (if b^512 > number) to split
    1289             # the number in half and recurse
    1290 
    1291             #   Pre-computing the exact number of digits up-front is actually faster (especially for large
    1292             # values of self) than trimming off trailing zeros after the fact.  It also seems that it would
     1281            # We use a divide and conquer approach (suggested by the prior
     1282            # author, malb?, of the digits method) here: for base b, compute
     1283            # b^2, b^4, b^8, ... (repeated squaring) until you get larger
     1284            # than your number; then compute (n // b^256, n % b^256)
     1285            # (if b^512 > number) to split the number in half and recurse
     1286
     1287            # Pre-computing the exact number of digits up-front is actually
     1288            # faster (especially for large values of self) than trimming off
     1289            # trailing zeros after the fact.  It also seems that it would
    12931290            # avoid duplicating the list in memory with a list-slice.
    12941291            z = the_integer_ring._zero_element if digits is None else digits[0]
    12951292            s = self_abs.exact_log(_base)
    12961293            l = [z]*(s+1 if s+1 >= padto else padto)
    12971294
    1298             # set up digits for optimal access once we get inside the worker functions
     1295            # set up digits for optimal access once we get inside the worker
     1296            # functions
    12991297            if not digits is None:
    13001298                # list objects have fastest access in the innermost loop
    13011299                if not PyList_CheckExact(digits):
    13021300                    digits = [digits[i] for i in range(_base)]
    13031301            elif mpz_cmp_ui(_base.value,s) < 0 and mpz_cmp_ui(_base.value,10000):
    1304                 # We can get a speed boost by pre-allocating digit values in big cases
    1305                 # We do this we have more digits than the base and the base is not too
    1306                 # extremely large (currently, "extremely" means larger than 10000 --
    1307                 # that's very arbitrary.)
     1302                # We can get a speed boost by pre-allocating digit values in
     1303                # big cases.
     1304                # We do this we have more digits than the base and the base
     1305                # is not too extremely large (currently, "extremely" means
     1306                # larger than 10000 -- that's very arbitrary.)
    13081307                if mpz_sgn(self.value) > 0:
    13091308                    digits = [Integer(i) for i in range(_base)]
    13101309                else:
     
    13261325                for power_index from 1 <= power_index < i:
    13271326                    power_list[power_index] = power_list[power_index-1]**2
    13281327
    1329                 #   Note that it may appear that the recursive calls to _digit_internal would be assigning
    1330                 # list elements i in l for anywhere from 0<=i<(1<<power_index).  However, this is not
    1331                 # the case due to the optimization of skipping assigns assigning zero.
     1328                # Note that it may appear that the recursive calls to
     1329                # _digit_internal would be assigning list elements i in l for
     1330                # anywhere from 0<=i<(1<<power_index).  However, this is not
     1331                # the case due to the optimization of skipping assigns
     1332                # assigning zero.
    13321333                _digits_internal(self.value,l,0,i-1,power_list,digits)
    13331334
    13341335            if s > 256:
     
    13441345       
    13451346        INPUT:
    13461347       
    1347        
    13481348        -  ``base`` - integer (default: 10)
    13491349       
    1350        
    13511350        EXAMPLES::
    13521351       
    13531352            sage: n = 52
     
    13591358            sage: n = 15
    13601359            sage: n.ndigits(2)
    13611360            4
    1362             sage: n=1000**1000000+1
     1361            sage: n = 1000**1000000+1
    13631362            sage: n.ndigits()
    13641363            3000001
    1365             sage: n=1000**1000000-1
     1364            sage: n = 1000**1000000-1
    13661365            sage: n.ndigits()
    13671366            3000000
    1368             sage: n=10**10000000-10**9999990
     1367            sage: n = 10**10000000-10**9999990
    13691368            sage: n.ndigits()
    13701369            10000000
    13711370        """
    13721371        cdef Integer temp
    1373         if mpz_cmp_si(self.value,0) == 0:
     1372        cdef unsigned long b
     1373
     1374        if mpz_sgn(self.value) == 0:
    13741375            temp = PY_NEW(Integer)
    1375             mpz_set_ui(temp.value,0)
     1376            mpz_set_ui(temp.value, 1)
    13761377            return temp
    1377         elif mpz_cmp_si(self.value,0) > 0:
    1378             temp=self.exact_log(base)
    1379             mpz_add_ui(temp.value,temp.value,1)
     1378
     1379        if mpz_sgn(self.value) > 0:
     1380            temp = self.exact_log(base)
     1381            mpz_add_ui(temp.value, temp.value, 1)
    13801382            return temp
    13811383        else:
    13821384            return self.abs().exact_log(base) + 1
     
    15841586            sage: Integer(32) / Integer(32)
    15851587            1
    15861588        """
    1587         # this is vastly faster than doing it here, since here
     1589        # This is vastly faster than doing it here, since here
    15881590        # we can't cimport rationals.
    15891591        return the_integer_ring._div(self, right)
    15901592
     
    16541656        else:
    16551657            return bin_op(x, y, operator.floordiv)
    16561658
    1657 
    16581659    def __pow__(self, n, modulus):
    16591660        r"""
    16601661        Computes `\text{self}^n`
     
    17411742            from sage.rings.finite_rings.integer_mod import Mod
    17421743            return Mod(self, modulus) ** n
    17431744       
    1744        
    17451745        if not PY_TYPE_CHECK(self, Integer):
    17461746            if isinstance(self, str):
    17471747                return self * n
     
    17801780        else:
    17811781            return x
    17821782
    1783 
    17841783    def nth_root(self, int n, bint truncate_mode=0):
    17851784        r"""
    17861785        Returns the (possibly truncated) n'th root of self.
    17871786       
    17881787        INPUT:
    17891788       
    1790        
    17911789        -  ``n`` - integer >= 1 (must fit in C int type).
    17921790       
    17931791        -  ``truncate_mode`` - boolean, whether to allow truncation if
    17941792           self is not an n'th power.
    17951793       
    1796        
    1797         OUTPUT: If truncate_mode is 0 (default), then returns the
    1798         exact n'th root if self is an n'th power, or raises a
    1799         ValueError if it is not.
     1794        OUTPUT:
     1795
     1796        If truncate_mode is 0 (default), then returns the exact n'th root
     1797        if self is an n'th power, or raises a ValueError if it is not.
    18001798       
    18011799        If truncate_mode is 1, then if either n is odd or self is
    18021800        positive, returns a pair (root, exact_flag) where root is the
     
    18991897
    19001898    cpdef size_t _exact_log_log2_iter(self,Integer m):
    19011899        """
    1902         This is only for internal use only.  You should expect it to crash and burn for
    1903         negative or other malformed input.  In particular, if the base `2\leq m<4` the log2
    1904         approximation of m is 1 and certain input causes endless loops.  Along these lines,
    1905         it is clear that this function is most useful for m with a relatively large number
     1900        This is only for internal use only.  You should expect it to crash
     1901        and burn for negative or other malformed input.  In particular, if
     1902        the base `2 \leq m < 4` the log2 approximation of m is 1 and certain
     1903        input causes endless loops.  Along these lines, it is clear that
     1904        this function is most useful for m with a relatively large number
    19061905        of bits.
    19071906
    1908         For ``small`` values (which I'll leave quite ambiguous), this function is a fast
    1909         path for exact log computations.  Any integer division with such input tends to
    1910         dominate the runtime.  Thus we avoid division entirely in this function.
     1907        For ``small`` values (which I'll leave quite ambiguous), this function
     1908        is a fast path for exact log computations.  Any integer division with
     1909        such input tends to dominate the runtime.  Thus we avoid division
     1910        entirely in this function.
    19111911
    19121912        AUTHOR::
    19131913
     
    19291929        cdef mpz_t accum
    19301930        cdef mpz_t temp_exp
    19311931
    1932         if mpz_cmp_si(m.value,4)<0:
     1932        if mpz_cmp_si(m.value,4) < 0:
    19331933            raise ValueError, "This is undefined or possibly non-convergent with this algorithm."
    19341934
    19351935        n_log2=mpz_sizeinbase(self.value,2)-1
     
    19431943            mpz_set_ui(accum,1)
    19441944            l = 0
    19451945            while l_min != l_max:
    1946                 #print "self=...",m,l_min,l_max
     1946                # print "self=...",m,l_min,l_max
    19471947                if l_min + 1 == l_max:
    19481948                    mpz_pow_ui(temp_exp,m.value,l_min+1-l)
    1949                     mpz_mul(accum,accum,temp_exp) # this might over-shoot and make accum > self, but we'll know that it's only over by a factor of m^1.
     1949                    # This might over-shoot and make accum > self, but
     1950                    # we'll know that it's only over by a factor of m^1.
     1951                    mpz_mul(accum,accum,temp_exp)
    19501952                    if mpz_cmp(self.value,accum) >= 0:
    19511953                        l_min += 1
    19521954                    break
     
    19541956                mpz_mul(accum,accum,temp_exp)
    19551957                l = l_min
    19561958
    1957                 #Let x=n_log2-(mpz_sizeinbase(accum,2)-1) and y=m_log2.  Now, with x>0 and y>0, we have the
    1958                 #following observation.  If floor((x-1)/(y+1))=0, then x-1<y+1 which implies that x/y<1+2/y. 
    1959                 #So long as y>=2, this means that floor(x/y)<=1.  This shows that this iteration is forced to
    1960                 #converge for input m>=4.  If m=3, we can find input so that floor((x-1)/(y+1))=0 and floor(x/y)=2
    1961                 #which results in non-convergence.
    1962 
    1963                 # We need the additional '-1' in the l_min computation because mpz_sizeinbase(accum,2)-1 is smaller than the true log_2(accum)
     1959                # Let x=n_log2-(mpz_sizeinbase(accum,2)-1) and y=m_log2. 
     1960                # Now, with x>0 and y>0, we have the following observation. 
     1961                # If floor((x-1)/(y+1))=0, then x-1<y+1 which implies that
     1962                # x/y<1+2/y. 
     1963                # So long as y>=2, this means that floor(x/y)<=1.  This shows
     1964                # that this iteration is forced to converge for input m >= 4. 
     1965                # If m=3, we can find input so that floor((x-1)/(y+1))=0 and
     1966                # floor(x/y)=2 which results in non-convergence.
     1967
     1968                # We need the additional '-1' in the l_min computation
     1969                # because mpz_sizeinbase(accum,2)-1 is smaller than the
     1970                # true log_2(accum)
    19641971                l_min=l+(n_log2-(mpz_sizeinbase(accum,2)-1)-1)/(m_log2+1)
    19651972                l_max=l+(n_log2-(mpz_sizeinbase(accum,2)-1))/m_log2
    19661973            mpz_clear(temp_exp)
     
    19701977
    19711978    cpdef size_t _exact_log_mpfi_log(self,m):
    19721979        """
    1973         This is only for internal use only.  You should expect it to crash and burn for
    1974         negative or other malformed input.
    1975 
    1976         I avoid using this function until the input is large.  The overhead associated
    1977         with computing the floating point log entirely dominates the runtime for small values.
    1978         Note that this is most definitely not an artifact of format conversion.  Tricks
    1979         with log2 approximations and using exact integer arithmetic are much better for
    1980         small input.
     1980        This is only for internal use only.  You should expect it to crash
     1981        and burn for negative or other malformed input.
     1982
     1983        I avoid using this function until the input is large.  The overhead
     1984        associated with computing the floating point log entirely dominates
     1985        the runtime for small values.  Note that this is most definitely not
     1986        an artifact of format conversion.  Tricks with log2 approximations
     1987        and using exact integer arithmetic are much better for small input.
    19811988
    19821989        AUTHOR::
    19831990
     
    20112018            # I'm not sure what to do now
    20122019            upper = 0
    20132020        lower = rif_log.lower().floor()
    2014         # since the log function is monotonic increasing, lower and upper bracket our desired answer
     2021        # since the log function is monotonic increasing, lower
     2022        # and upper bracket our desired answer
    20152023
    20162024        # if upper - lower == 1: "we are done"
    20172025        if upper - lower == 2:
    2018             # You could test it by checking rif_m**(lower+1), but I think that's a waste of time since it won't be conclusive
    2019             # We must test with exact integer arithmetic which takes all the bits of self into account.
     2026            # You could test it by checking rif_m**(lower+1), but I think
     2027            # that's a waste of time since it won't be conclusive.
     2028            # We must test with exact integer arithmetic which takes all
     2029            # the bits of self into account.
    20202030            if self >= m**(lower+1):
    20212031                return lower + 1
    20222032            else:
     
    20602070
    20612071    def exact_log(self, m):
    20622072        r"""
    2063         Returns the largest integer `k` such that
    2064         `m^k \leq \text{self}`, i.e., the floor of
    2065         `\log_m(\text{self})`.
    2066        
     2073        Returns the largest integer `k` such that `m^k \leq \text{self}`,
     2074        i.e., the floor of `\log_m(\text{self})`.
     2075
    20672076        This is guaranteed to return the correct answer even when the usual
    20682077        log function doesn't have sufficient precision.
    20692078       
    20702079        INPUT:
    20712080       
    2072        
    20732081        -  ``m`` - integer >= 2
    20742082       
    2075        
    20762083        AUTHORS:
    20772084
    20782085        - David Harvey (2006-09-15)
     
    21402147        cdef size_t n_log2
    21412148        cdef size_t m_log2
    21422149        cdef size_t guess # this will contain the final answer
    2143         cdef bint guess_filled = 0  # this variable is only used in one branch below.
     2150        cdef bint guess_filled = 0  # this variable is only used in one branch below
    21442151        cdef mpz_t z
    21452152        if PY_TYPE_CHECK(m,Integer):
    21462153            _m=<Integer>m
     
    21552162        n_log2=mpz_sizeinbase(self.value,2)-1
    21562163        m_log2=mpz_sizeinbase(_m.value,2)-1
    21572164        if mpz_divisible_2exp_p(_m.value,m_log2):
    2158             # Here, m is a power of 2 and the correct answer is found by a log 2 approximation.
    2159             guess = n_log2/m_log2 # truncating division
     2165            # Here, m is a power of 2 and the correct answer is found
     2166            # by a log 2 approximation.
     2167            guess = n_log2/m_log2  # truncating division
    21602168        elif n_log2/(m_log2+1) == n_log2/m_log2:
    2161             # In this case, we have an upper bound and lower bound which give the same answer, thus, the correct answer.
     2169            # In this case, we have an upper bound and lower bound which
     2170            # give the same answer, thus, the correct answer.
    21622171            guess = n_log2/m_log2
    2163         elif m_log2 < 8: # i.e. m<256
     2172        elif m_log2 < 8:  # i.e. m<256
    21642173            # if the base m is at most 256, we can use mpz_sizeinbase
    21652174            # to get the following guess which is either the exact
    21662175            # log, or 1+ the exact log
     
    21832192                    guess =  guess - 1
    21842193            if not guess_filled:
    21852194                # At this point, either
    2186                 #  1)  self is close enough to a perfect power of m that we need an exact comparison, or
    2187                 #  2)  the numbers are small enough that converting to the interval field is more work than the exact comparison.
     2195                #  1)  self is close enough to a perfect power of m that we
     2196                #      need an exact comparison, or
     2197                #  2)  the numbers are small enough that converting to the
     2198                #      interval field is more work than the exact comparison.
    21882199                compare = _m**guess
    21892200                if self < compare:
    21902201                    guess = guess - 1
    21912202        elif n_log2 < 5000:
    2192             # for input with small exact log, it's very fast to work in exact integer arithmetic starting from log2 approximations
     2203            # for input with small exact log, it's very fast to work in exact
     2204            # integer arithmetic starting from log2 approximations
    21932205            guess = self._exact_log_log2_iter(_m)
    21942206        else:
    2195             # finally, we are out of easy cases
    2196             # this subroutine uses interval arithmetic to guess and check the exact log.
     2207            # finally, we are out of easy cases this subroutine uses interval
     2208            # arithmetic to guess and check the exact log.
    21972209            guess = self._exact_log_mpfi_log(_m)
    21982210
    2199         result = <Integer>PY_NEW(Integer)
     2211        result = PY_NEW(Integer)
    22002212        mpz_set_ui(result.value,guess)
    22012213        return result
    22022214
     
    22122224       
    22132225        INPUT:
    22142226       
    2215        
    22162227        -  ``m`` - default: natural log base e
    22172228       
    22182229        -  ``prec`` - integer (default: None): if None, returns
    22192230           symbolic, else to given bits of precision as in RealField
    22202231       
    2221        
    22222232        EXAMPLES::
    22232233       
    22242234            sage: Integer(124).log(5)
     
    23372347       
    23382348        INPUT:
    23392349       
    2340        
    23412350        -  ``m`` - Integer
    23422351       
    2343        
    23442352        OUTPUT: Integer
    23452353       
    23462354        EXAMPLES::
     
    24252433            raise ValueError, "n must be nonzero"
    24262434        f = self.factor()
    24272435
    2428         # All of the declarations below are for optimizing the word-sized case.
    2429         # Operations are performed in c as far as possible without overflow before
    2430         # moving to python objects.
     2436        # All of the declarations below are for optimizing the word-sized
     2437        # case.  Operations are performed in c as far as possible without
     2438        # overflow before moving to python objects.
    24312439        cdef long long p_c, pn_c, apn_c
    24322440        cdef long all_len, sorted_len, prev_len
    24332441        cdef long long* ptr
     
    24682476                sage_free(ptr)
    24692477                fits_c = False
    24702478
    2471             # The two cases below are essentially the same algorithm, one operating
    2472             # on Integers in Python lists, the other on long longs.
     2479            # The two cases below are essentially the same algorithm, one
     2480            # operating on Integers in Python lists, the other on long longs.
    24732481            if fits_c:
    24742482               
    24752483                pn_c = p_c = p
     
    25672575            sage: abs(z) == abs(1)
    25682576            True
    25692577        """
    2570         cdef Integer x
    2571         x = PY_NEW(Integer)
     2578        cdef Integer x = PY_NEW(Integer)
    25722579        mpz_abs(x.value, self.value)
    25732580        return x
    25742581
     
    26102617            sage: a % b
    26112618            59
    26122619         """
    2613         cdef Integer z = <Integer>PY_NEW(Integer)
     2620        cdef Integer z = PY_NEW(Integer)
    26142621        cdef long yy, res
    26152622
    26162623        # first case: Integer % Integer
     
    26502657            except ValueError:
    26512658                return bin_op(x, y, operator.mod)
    26522659
    2653 
    26542660    def quo_rem(Integer self, other):
    26552661        """
    26562662        Returns the quotient and the remainder of self divided by other.
     
    26592665       
    26602666        INPUT:
    26612667       
    2662        
    26632668        -  ``other`` - the integer the divisor
    26642669       
    2665        
    26662670        OUTPUT:
    26672671       
    26682672       
     
    26702674       
    26712675        -  ``r`` - the remainder of self/other
    26722676       
    2673        
    26742677        EXAMPLES::
    26752678       
    26762679            sage: z = Integer(231)
     
    33813384            return sage.rings.infinity.infinity
    33823385        if mpz_cmp_ui(p.value, 2) < 0:
    33833386            raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
    3384         cdef Integer v
    3385         v = PY_NEW(Integer)
     3387
     3388        cdef Integer v = PY_NEW(Integer)
    33863389        cdef mpz_t u
    33873390        mpz_init(u)
    33883391        _sig_on
     
    34233426       
    34243427        INPUT:
    34253428       
    3426        
    34273429        -  ``p`` - an integer at least 2.
    34283430       
    3429        
    34303431        EXAMPLE::
    34313432       
    34323433            sage: n = 60
     
    34483449        """
    34493450        return self._valuation(Integer(p))
    34503451
     3452    # Alias for valuation
     3453    ord = valuation
     3454
    34513455    def val_unit(self, p):
    34523456        r"""
    34533457        Returns a pair: the p-adic valuation of self, and the p-adic unit
     
    34553459       
    34563460        INPUT:
    34573461       
    3458        
    34593462        -  ``p`` - an integer at least 2.
    34603463       
    3461        
    34623464        OUTPUT:
    34633465       
    3464        
    3465         -  ``v_p(self)`` - the p-adic valuation of
    3466            ``self``
    3467        
    3468         -  ``u_p(self)`` - ``self`` /
    3469            `p^{v_p(\mathrm{self})}`
    3470        
     3466        -  ``v_p(self)`` - the p-adic valuation of ``self``
     3467       
     3468        -  ``u_p(self)`` - ``self`` / `p^{v_p(\mathrm{self})}`
    34713469       
    34723470        EXAMPLE::
    34733471       
     
    34853483        """
    34863484        return self._val_unit(Integer(p))
    34873485
    3488     def ord(self, p=None):
    3489         """
    3490         Synonym for valuation
    3491        
     3486    def odd_part(self):
     3487        r"""
     3488        The odd part of the integer `n`. This is `n / 2^v`,
     3489        where `v = \mathrm{valuation}(n,2)`.
     3490
     3491        IMPLEMENTATION:
     3492
     3493        Currently returns 0 when self is 0.  This behaviour is fairly arbitrary,
     3494        and in Sage 4.6 this special case was not handled at all, eventually
     3495        propagating a TypeError.  The caller should not rely on the behaviour
     3496        in case self is 0.
     3497
    34923498        EXAMPLES::
    34933499       
    3494             sage: n=12
    3495             sage: n.ord(3)
     3500            sage: odd_part(5)
     3501            5
     3502            sage: odd_part(4)
    34963503            1
    3497         """
    3498         return self.valuation(p)
     3504            sage: odd_part(factorial(31))
     3505            122529844256906551386796875
     3506        """
     3507        cdef Integer odd
     3508        cdef unsigned long bits
     3509
     3510        if mpz_cmpabs_ui(self.value, 1) <= 0:
     3511            return self
     3512
     3513        odd  = PY_NEW(Integer)
     3514        bits = mpz_scan1(self.value, 0)
     3515        mpz_tdiv_q_2exp(odd.value, self.value, bits)
     3516        return odd
    34993517
    35003518    cdef Integer _divide_knowing_divisible_by(Integer self, Integer right):
    35013519        r"""
     
    35633581            sage: n._lcm(150)
    35643582            300
    35653583        """
    3566         cdef Integer z
    3567         z = PY_NEW(Integer)
     3584        cdef Integer z = PY_NEW(Integer)
    35683585        _sig_on
    35693586        mpz_lcm(z.value, self.value, n.value)
    35703587        _sig_off
     
    35893606        """
    35903607        Return the numerator of this integer.
    35913608       
    3592         EXAMPLE::
     3609        EXAMPLES::
    35933610       
    35943611            sage: x = 5
    35953612            sage: x.numerator()
     
    36883705                sign = ONE
    36893706            return sign / Integer(-k-n).multifactorial(k)
    36903707               
    3691         # compute the actual product, optimizing the number of large multiplications
     3708        # compute the actual product, optimizing the number of large
     3709        # multiplications
    36923710        cdef int i,j
    36933711       
    36943712        # we need (at most) log_2(#factors) concurrent sub-products
     
    37273745        free(sub_prods)
    37283746
    37293747        return z
    3730 
    37313748       
    37323749    def gamma(self):
    37333750        r"""
    3734         The gamma function on integers is the factorial function (shifted
    3735         by one) on positive integers, and `\pm \infty` on
    3736         non-positive integers.
     3751        The gamma function on integers is the factorial function (shifted by
     3752        one) on positive integers, and `\pm \infty` on non-positive integers.
    37373753       
    37383754        EXAMPLES::
    37393755       
     
    38023818       
    38033819    def is_one(self):
    38043820        r"""
    3805         Returns ``True`` if the integer is `1`,
    3806         otherwise ``False``.
     3821        Returns ``True`` if the integer is `1`, otherwise ``False``.
    38073822       
    38083823        EXAMPLES::
    38093824       
     
    38163831
    38173832    def __nonzero__(self):
    38183833        r"""
    3819         Returns ``True`` if the integer is not `0`,
    3820         otherwise ``False``.
     3834        Returns ``True`` if the integer is not `0`, otherwise ``False``.
    38213835       
    38223836        EXAMPLES::
    38233837       
     
    38423856
    38433857    def is_unit(self):
    38443858        r"""
    3845         Returns ``true`` if this integer is a unit, i.e., 1 or
    3846         `-1`.
     3859        Returns ``true`` if this integer is a unit, i.e., 1 or `-1`.
    38473860       
    38483861        EXAMPLES::
    38493862       
     
    38553868            1 True
    38563869            2 False
    38573870        """
    3858         return mpz_cmp_si(self.value, -1) == 0 or mpz_cmp_si(self.value, 1) == 0
     3871        return mpz_cmpabs_ui(self.value, 1) == 0
    38593872
    38603873    def is_square(self):
    38613874        r"""
    3862         Returns ``True`` if self is a perfect square
     3875        Returns ``True`` if self is a perfect square.
    38633876       
    38643877        EXAMPLES::
    38653878       
     
    38723885
    38733886    def is_power(self):
    38743887        r"""
    3875         Returns ``True`` if self is a perfect power, ie if
    3876         there exist integers a and b, `b > 1` with
    3877         `self = a^b`.
     3888        Returns ``True`` if self is a perfect power, i.e. if there exist
     3889        integers `a` and `b`, `b > 1` with `self = a^b`.
    38783890       
    38793891        EXAMPLES::
    38803892       
     
    38873899
    38883900    cdef bint _is_power_of(Integer self, Integer n):
    38893901        r"""
    3890         Returns a non-zero int if there is an integer b with
     3902        Returns a non-zero int if there is an integer b with 
    38913903        `\mathtt{self} = n^b`.
    38923904       
    3893         For more documentation see ``is_power_of``
     3905        For more documentation see ``is_power_of``.
    38943906       
    38953907        AUTHORS:
    38963908
     
    41344146            return True
    41354147        if not self.is_perfect_power():
    41364148            return False
    4137         return len(self.factor()) == 1  # this is potentially very slow, but at least it is right!!
     4149
     4150        # This is potentially very slow, but at least it is right!!
     4151        return len(self.factor()) == 1
     4152
    41384153        # PARI has a major bug -- see trac #4777. It gives wrong answer on
    41394154        # input of "150607571^14".
    41404155        #k, g = self._pari_().ispower()
     
    41984213            sage: z.is_irreducible()
    41994214            True
    42004215        """
    4201         n = self if self>=0 else -self
     4216        n = self if self >= 0 else -self
    42024217        return bool(n._pari_().isprime())
    42034218
    42044219    def is_pseudoprime(self):
     
    42344249            sage: (-4).is_perfect_power()
    42354250            False
    42364251
    4237         This is a test to make sure we workaround a bug in GMP. (See
    4238         trac \#4612.)
     4252        This is a test to make sure we work around a bug in GMP,
     4253        see trac \#4612.
    42394254
    42404255        ::
    42414256
     
    42444259        """
    42454260        cdef mpz_t tmp
    42464261        cdef int res
    4247         if mpz_sgn(self.value) == -1:
     4262        if mpz_sgn(self.value) < 0:
    42484263            if mpz_cmp_si(self.value, -1) == 0:
    42494264                return True
    42504265            mpz_init(tmp)
    4251             mpz_mul_si(tmp, self.value, -1)
     4266            mpz_neg(tmp, self.value)
    42524267            while mpz_perfect_square_p(tmp):
    42534268                mpz_sqrt(tmp, tmp)
    42544269            res = mpz_perfect_power_p(tmp)
     
    42984313       
    42994314    def kronecker(self, b):
    43004315        r"""
    4301         Calculate the Kronecker symbol
    4302         `\left(\frac{self}{b}\right)` with the Kronecker extension
    4303         `(self/2)=(2/self)` when self odd, or `(self/2)=0`
    4304         when `self` even.
     4316        Calculate the Kronecker symbol `\left(\frac{self}{b}\right)`
     4317        with the Kronecker extension `(self/2)=(2/self)` when `self` is odd,
     4318        or `(self/2)=0` when `self` is even.
    43054319       
    43064320        EXAMPLES::
    43074321       
     
    43534367       
    43544368    def squarefree_part(self, long bound=-1):
    43554369        r"""
    4356         Return the square free part of `x` (=self), i.e., the
    4357         unique integer `z` that `x = z y^2`, with
    4358         `y^2` a perfect square and `z` square-free.
    4359        
    4360         Use ``self.radical()`` for the product of the primes
    4361         that divide self.
     4370        Return the square free part of `x` (=self), i.e., the unique integer
     4371        `z` that `x = z y^2`, with `y^2` a perfect square and `z` square-free.
     4372       
     4373        Use ``self.radical()`` for the product of the primes that divide self.
    43624374       
    43634375        If self is 0, just returns 0.
    43644376       
     
    43974409        cdef Integer z
    43984410        cdef long even_part, p, p2
    43994411        cdef char switch_p
    4400         if mpz_cmp_ui(self.value, 0) == 0:
     4412        if mpz_sgn(self.value) == 0:
    44014413            return self
    44024414        if 0 <= bound < 2:
    44034415            return self
    44044416        elif 2 <= bound <= 10000:
    4405             z = <Integer>PY_NEW(Integer)
     4417            z = PY_NEW(Integer)
    44064418            even_part = mpz_scan1(self.value, 0)
    44074419            mpz_fdiv_q_2exp(z.value, self.value, even_part ^ (even_part&1))
    44084420            _sig_on
     
    44584470       
    44594471        INPUT:
    44604472       
    4461        
    44624473        -  ``proof`` - bool or None (default: None, see
    44634474           proof.arithmetic or sage.structure.proof) Note that the global Sage
    44644475           default is proof=True
    44654476       
    4466        
    44674477        EXAMPLES::
    44684478       
    44694479            sage: Integer(100).next_prime()
     
    44964506            n += 2
    44974507        return integer_ring.ZZ(n)
    44984508   
    4499 
    45004509    def additive_order(self):
    45014510        """
    45024511        Return the additive order of self.
     
    46254634        """
    46264635        Return string that evaluates in Magma to this element.
    46274636
    4628         For small integers we just use base 10.  For large integers we use base 16,
    4629         but use Magma's StringToInteger command, which (for no good reason)
    4630         is much faster than 0x[string literal].  We only use base 16 for integers
    4631         with at least 10000 binary digits, since e.g., for a large list of
    4632         small integers the overhead of calling StringToInteger can be a killer.
     4637        For small integers we just use base 10.  For large integers we use
     4638        base 16, but use Magma's StringToInteger command, which (for no
     4639        good reason) is much faster than 0x[string literal].  We only use
     4640        base 16 for integers with at least 10000 binary digits, since e.g.,
     4641        for a large list of small integers the overhead of calling
     4642        StringToInteger can be a killer.
    46334643
    46344644        EXAMPLES:
    46354645            sage: (117)._magma_init_(magma)           # optional - magma
     
    47394749        """
    47404750        if mpz_sgn(self.value) < 0:
    47414751            raise ValueError, "square root of negative integer not defined."
    4742         cdef Integer x
    4743         x = PY_NEW(Integer)
     4752
     4753        cdef Integer x = PY_NEW(Integer)
    47444754
    47454755        _sig_on
    47464756        mpz_sqrt(x.value, self.value)
     
    47764786       
    47774787        INPUT:
    47784788       
    4779        
    47804789        -  ``prec`` - integer (default: None): if None, returns
    47814790           an exact square root; otherwise returns a numerical square root if
    47824791           necessary, to the given bits of precision.
     
    47894798        -  ``all`` - bool (default: False); if True, return all
    47904799           square roots of self, instead of just one.
    47914800       
    4792        
    47934801        EXAMPLES::
    47944802       
    47954803            sage: Integer(144).sqrt()
     
    48644872       
    48654873        INPUT:
    48664874       
    4867        
    48684875        -  ``self`` - integer
    48694876       
    48704877        -  ``n`` - integer
     
    48724879        -  ``minimal`` - boolean (default false), whether to
    48734880           compute minimal cofactors (see below)
    48744881       
    4875        
    4876         OUTPUT: A triple (g, s, t), where `g` is the non-negative
    4877         gcd of self and `n`, and `s` and `t` are
    4878         cofactors satisfying the Bezout identity
     4882        OUTPUT:
     4883
     4884        A triple (g, s, t), where `g` is the non-negative gcd of self
     4885        and `n`, and `s` and `t` are cofactors satisfying the Bezout identity
    48794886       
    48804887        .. math::
    48814888       
    48824889             g = s \cdot \mbox{\rm self} + t \cdot n.
    48834890       
    4884        
    4885        
    48864891        .. note::
    48874892
    4888            If minimal is False, then there is no guarantee that the
    4889            returned cofactors will be minimal in any sense; the only
    4890            guarantee is that the Bezout identity will be satisfied
    4891            (see examples below).
    4892        
    4893            If minimal is True, the cofactors will satisfy the
    4894            following conditions. If either self or `n` are zero,
    4895            the trivial solution is returned. If both self and
    4896            `n` are nonzero, the function returns the unique
    4897            solution such that `0 \leq s < |n|/g` (which then
    4898            must also satisfy
    4899            `0 \leq |t| \leq |\mbox{\rm self}|/g`).
     4893           If minimal is False, then there is no guarantee that the returned
     4894           cofactors will be minimal in any sense; the only guarantee is that
     4895           the Bezout identity will be satisfied (see examples below).
     4896           
     4897           If minimal is True, the cofactors will satisfy the following
     4898           conditions. If either self or `n` are zero, the trivial solution
     4899           is returned. If both self and `n` are nonzero, the function returns
     4900           the unique solution such that `0 \leq s < |n|/g` (which then must
     4901           also satisfy `0 \leq |t| \leq |\mbox{\rm self}|/g`).
    49004902       
    49014903        EXAMPLES::
    49024904       
     
    50705072        n *= sign
    50715073
    50725074        # Now finally call into MPIR to do the shifting.
    5073         cdef Integer z = <Integer>PY_NEW(Integer)
     5075        cdef Integer z = PY_NEW(Integer)
    50745076        if n < 0:
    50755077            mpz_fdiv_q_2exp(z.value, self.value, -n)
    50765078        else:
     
    51415143        return (<Integer>x)._shift_helper(y, -1)
    51425144
    51435145    cdef _and(Integer self, Integer other):
    5144         cdef Integer x
    5145         x = PY_NEW(Integer)
     5146        cdef Integer x = PY_NEW(Integer)
    51465147        mpz_and(x.value, self.value, other.value)
    51475148        return x
    51485149
     
    51625163            return (<Integer>x)._and(y)
    51635164        return bin_op(x, y, operator.and_)
    51645165
    5165 
    51665166    cdef _or(Integer self, Integer other):
    5167         cdef Integer x
    5168         x = PY_NEW(Integer)
     5167        cdef Integer x = PY_NEW(Integer)
    51695168        mpz_ior(x.value, self.value, other.value)
    51705169        return x
    51715170
     
    51835182            return (<Integer>x)._or(y)
    51845183        return bin_op(x, y, operator.or_)
    51855184
    5186 
    51875185    def __invert__(self):
    51885186        """
    51895187        Return the multiplicative inverse of self, as a rational number.
     
    51965194            sage: n.__invert__()
    51975195            1/10
    51985196        """
    5199         return one/self    # todo: optimize
    5200        
     5197        return one / self
     5198
    52015199    def inverse_of_unit(self):
    52025200        """
    52035201        Return inverse of self if self is a unit in the integers, i.e.,
     
    52185216            ...
    52195217            ZeroDivisionError: Inverse does not exist.
    52205218        """
    5221         if mpz_cmp_si(self.value, 1) == 0 or mpz_cmp_si(self.value, -1) == 0:
     5219        if mpz_cmpabs_ui(self.value, 1) == 0:
    52225220            return self
    52235221        else:
    52245222            raise ZeroDivisionError, "Inverse does not exist."
    52255223
    52265224    def inverse_mod(self, n):
    52275225        """
    5228         Returns the inverse of self modulo `n`, if this inverse
    5229         exists. Otherwise, raises a ``ZeroDivisionError``
    5230         exception.
     5226        Returns the inverse of self modulo `n`, if this inverse exists.
     5227        Otherwise, raises a ``ZeroDivisionError`` exception.
    52315228       
    52325229        INPUT:
    52335230       
    5234        
    52355231        -  ``self`` - Integer
    52365232       
    52375233        -  ``n`` - Integer
    52385234       
    5239        
    52405235        OUTPUT:
    52415236       
    5242        
    52435237        -  ``x`` - Integer such that x\*self = 1 (mod m), or
    52445238           raises ZeroDivisionError.
    52455239
     
    52985292            left, right = canonical_coercion(self, n)
    52995293            return left.gcd(right)
    53005294        cdef Integer m = as_Integer(n)
    5301         cdef Integer g = <Integer>PY_NEW(Integer)
     5295        cdef Integer g = PY_NEW(Integer)
    53025296        _sig_on
    53035297        mpz_gcd(g.value, self.value, m.value)
    53045298        _sig_off
     
    53065300
    53075301    def crt(self, y, m, n):
    53085302        """
    5309         Return the unique integer between `0` and `mn` that
    5310         is congruent to the integer modulo `m` and to `y`
    5311         modulo `n`. We assume that `m` and `n` are
    5312         coprime.
     5303        Return the unique integer between `0` and `mn` that is congruent to
     5304        the integer modulo `m` and to `y` modulo `n`. We assume that `m` and
     5305        `n` are coprime.
    53135306       
    53145307        EXAMPLES::
    53155308       
     
    53285321        if not g.is_one():
    53295322            raise ArithmeticError, "CRT requires that gcd of moduli is 1."
    53305323        # Now s*m + t*n = 1, so the answer is x + (y-x)*s*m, where x=self.
    5331         return (self + (_y-self)*s*_m) % (_m*_n)
     5324        return (self + (_y - self) * s * _m) % (_m * _n)
    53325325
    53335326    def test_bit(self, long index):
    53345327        r"""
    53355328        Return the bit at ``index``.
    53365329       
     5330        If the index is negative, returns 0.
     5331       
     5332        Although internally a sign-magnitude representation is used
     5333        for integers, this method pretends to use a two's complement
     5334        representation.  This is illustrated with a negative integer
     5335        below.
     5336       
    53375337        EXAMPLES::
    53385338       
    53395339            sage: w = 6
     
    53435343            1
    53445344            sage: w.test_bit(-1)
    53455345            0
     5346            sage: x = -20
     5347            sage: x.str(2)
     5348            '-10100'
     5349            sage: x.test_bit(4)
     5350            0
     5351            sage: x.test_bit(5)
     5352            1
     5353            sage: x.test_bit(6)
     5354            1
    53465355        """
    53475356        if index < 0:
    53485357            return 0
     
    53955404   
    53965405    INPUT:
    53975406   
    5398    
    53995407    -  ``v`` - list or tuple
    54005408   
    5401    
    54025409    OUTPUT: integer
    54035410   
    54045411    EXAMPLES::
     
    54405447       
    54415448    return z
    54425449
    5443 
    5444 
    54455450def GCD_list(v):
    54465451    r"""
    54475452    Return the GCD of a list v of integers. Elements of v are converted
     
    54515456   
    54525457    INPUT:
    54535458   
    5454    
    54555459    -  ``v`` - list or tuple
    54565460   
    5457    
    54585461    OUTPUT: integer
    54595462   
    54605463    EXAMPLES::
     
    55185521        ...
    55195522        TypeError: expected string or Unicode object, sage.rings.integer.Integer found
    55205523    """
    5521     cdef Integer r
    5522     r = PY_NEW(Integer)
     5524    cdef Integer r = PY_NEW(Integer)
    55235525    r._reduce_set(s)
    55245526    return r
    55255527   
     
    55635565    def _repr_type(self):
    55645566        return "Native"
    55655567
    5566 
    55675568############### INTEGER CREATION CODE #####################
    55685569
    55695570# We need a couple of internal GMP datatypes.
     
    58285829    integer_pool_size = 0
    58295830    integer_pool_count = 0
    58305831    sage_free(integer_pool)
     5832
  • sage/rings/integer_ring.pyx

    diff -r b6a68db0aac8 sage/rings/integer_ring.pyx
    a b  
    280280        raise TypeError, 'len() of unsized object'
    281281
    282282    def _div(self, integer.Integer left, integer.Integer right):
    283         cdef rational.Rational x
    284         x = PY_NEW(rational.Rational)
    285         if mpz_cmp_si(right.value, 0) == 0:
     283        cdef rational.Rational x = PY_NEW(rational.Rational)
     284        if mpz_sgn(right.value) == 0:
    286285            raise ZeroDivisionError, 'Rational division by zero'
    287         mpq_set_num(x.value, left.value)
    288         mpq_set_den(x.value, right.value)
     286        mpz_set(mpq_numref(x.value), left.value)
     287        mpz_set(mpq_denref(x.value), right.value)
    289288        mpq_canonicalize(x.value)
    290289        return x
    291290