Ticket #14416: 14416_QQ_to_RDF.patch

File 14416_QQ_to_RDF.patch, 9.2 KB (added by jdemeyer, 8 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1365664937 -7200
    # Node ID 3e090bfa3d68818d468f890766152dfd6fab63c2
    # Parent  c8ae9724d3cfad1837a130d50ac5c83a811ec221
    Fix conversion QQ -> RDF
    
    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b  
    10191019            sage: B.rank()
    10201020            12
    10211021            sage: A = B.change_ring(RDF)
    1022             sage: A.condition() > 1.7 * 10^16
     1022            sage: A.condition() > 1.6e16 or A.condition()
    10231023            True
    10241024
    10251025            sage: sv = A.singular_values(eps=None)
  • sage/matrix/matrix_mod2e_dense.pyx

    diff --git a/sage/matrix/matrix_mod2e_dense.pyx b/sage/matrix/matrix_mod2e_dense.pyx
    a b  
    804804            sage: A = matrix(K,1000,1000)
    805805            sage: A.randomize(nonzero=False, density=0.1)
    806806            sage: float(A.density())
    807             0.0937679...
     807            0.093768
    808808           
    809809            sage: A.randomize(nonzero=False, density=0.05)
    810810            sage: float(A.density())
  • sage/plot/colors.py

    diff --git a/sage/plot/colors.py b/sage/plot/colors.py
    a b  
    13021302        sage: rainbow(7)
    13031303        ['#ff0000', '#ffda00', '#48ff00', '#00ff91', '#0091ff', '#4800ff', '#ff00da']
    13041304        sage: rainbow(7, 'rgbtuple')
    1305         [(1.0, 0.0, 0.0), (1.0, 0.8571428571428571, 0.0), (0.2857142857142858, 1.0, 0.0), (0.0, 1.0, 0.5714285714285712), (0.0, 0.5714285714285716, 1.0), (0.2857142857142847, 0.0, 1.0), (1.0, 0.0, 0.8571428571428577)]
     1305        [(1.0, 0.0, 0.0), (1.0, 0.8571428571428571, 0.0), (0.2857142857142858, 1.0, 0.0), (0.0, 1.0, 0.5714285714285712), (0.0, 0.5714285714285716, 1.0), (0.2857142857142856, 0.0, 1.0), (1.0, 0.0, 0.8571428571428577)]
    13061306
    13071307    AUTHORS:
    13081308
  • sage/rings/contfrac.py

    diff --git a/sage/rings/contfrac.py b/sage/rings/contfrac.py
    a b  
    633633            sage: a = CFF(-17/389); a
    634634            [-1, 1, 21, 1, 7, 2]
    635635            sage: float(a)
    636             -0.04370179948586118
     636            -0.043701799485861184
    637637        """
    638638        return float(self._rational_())
    639639
  • sage/rings/rational.pyx

    diff --git a/sage/rings/rational.pyx b/sage/rings/rational.pyx
    a b  
    19991999            -0.23529411764705882
    20002000            sage: float(-4/17)
    20012001            -0.23529411764705882
     2002            sage: float(1/3)
     2003            0.3333333333333333
     2004            sage: float(1/10)
     2005            0.1
     2006
     2007        TESTS:
     2008
     2009        Test that conversion agrees with `RR`::
     2010
     2011            sage: Q = [a/b for a in [-99..99] for b in [1..99]]
     2012            sage: all([RDF(q) == RR(q) for q  in Q])
     2013            True
     2014
     2015        At some point, the continued fraction and direct conversion
     2016        to ``RDF`` should agree::
     2017
     2018            sage: Q = continued_fraction(pi, bits=3000).convergents()[20:]
     2019            sage: RDFpi = RDF(pi)
     2020            sage: all([RDF(q) == RDFpi for q  in Q])
     2021            True
    20022022        """
    2003         return mpq_get_d(self.value)
     2023        return mpq_get_d_nearest(self.value)
    20042024
    20052025    def __hash__(self):
    20062026        """
     
    34293449        if neg: v = -v
    34303450        return v
    34313451
     3452
     3453# The except value is just some random double, it doesn't matter what it is.
     3454cdef double mpq_get_d_nearest(mpq_t x) except? -648555075988944.0:
     3455    """
     3456    Convert a ``mpq_t`` to a ``double``, with round-to-even. This
     3457    differs from ``mpq_get_d()`` which does round-to-zero.
     3458
     3459    TESTS::
     3460
     3461        sage: q= QQ(); float(q)
     3462        0.0
     3463        sage: q = 2^-10000; float(q)
     3464        0.0
     3465        sage: float(-q)
     3466        -0.0
     3467        sage: q = 2^10000/1; float(q)
     3468        inf
     3469        sage: float(-q)
     3470        -inf
     3471
     3472    ::
     3473
     3474        sage: q = 2^-1075; float(q)
     3475        0.0
     3476        sage: float(-q)
     3477        -0.0
     3478        sage: q = 2^52 / 2^1074; float(q)  # Smallest normal double
     3479        2.2250738585072014e-308
     3480        sage: float(-q)
     3481        -2.2250738585072014e-308
     3482        sage: q = (2^52 + 1/2) / 2^1074; float(q)
     3483        2.2250738585072014e-308
     3484        sage: float(-q)
     3485        -2.2250738585072014e-308
     3486        sage: q = (2^52 + 1) / 2^1074; float(q)  # Next normal double
     3487        2.225073858507202e-308
     3488        sage: float(-q)
     3489        -2.225073858507202e-308
     3490        sage: q = (2^52 - 1) / 2^1074; float(q)  # Largest denormal double
     3491        2.225073858507201e-308
     3492        sage: float(-q)
     3493        -2.225073858507201e-308
     3494        sage: q = 1 / 2^1074; float(q)  # Smallest denormal double
     3495        5e-324
     3496        sage: float(-q)
     3497        -5e-324
     3498        sage: q = (1/2) / 2^1074; float(q)
     3499        0.0
     3500        sage: float(-q)
     3501        -0.0
     3502        sage: q = (3/2) / 2^1074; float(q)
     3503        1e-323
     3504        sage: float(-q)
     3505        -1e-323
     3506        sage: q = (2/3) / 2^1074; float(q)
     3507        5e-324
     3508        sage: float(-q)
     3509        -5e-324
     3510        sage: q = (1/3) / 2^1074; float(q)
     3511        0.0
     3512        sage: float(-q)
     3513        -0.0
     3514        sage: q = (2^53 - 1) * 2^971/1; float(q)  # Largest double
     3515        1.7976931348623157e+308
     3516        sage: float(-q)
     3517        -1.7976931348623157e+308
     3518        sage: q = (2^53) * 2^971/1; float(q)
     3519        inf
     3520        sage: float(-q)
     3521        -inf
     3522        sage: q = (2^53 - 1/2) * 2^971/1; float(q)
     3523        inf
     3524        sage: float(-q)
     3525        -inf
     3526        sage: q = (2^53 - 2/3) * 2^971/1; float(q)
     3527        1.7976931348623157e+308
     3528        sage: float(-q)
     3529        -1.7976931348623157e+308
     3530
     3531    AUTHORS:
     3532
     3533    - Paul Zimmermann, Jeroen Demeyer (:trac:`14416`)
     3534    """
     3535    # a = mpq_numref(x)
     3536    # b = mpq_denref(x)
     3537    cdef int resultsign = mpz_sgn(mpq_numref(x))
     3538
     3539    if resultsign == 0:
     3540        return 0.0
     3541
     3542    cdef Py_ssize_t sa = mpz_sizeinbase(mpq_numref(x), 2)
     3543    cdef Py_ssize_t sb = mpz_sizeinbase(mpq_denref(x), 2)
     3544   
     3545    # Easy case: both numerator and denominator are exactly
     3546    # representable as doubles.
     3547    if sa <= 53 and sb <= 53:
     3548        return mpz_get_d(mpq_numref(x)) / mpz_get_d(mpq_denref(x))
     3549   
     3550    # General case
     3551   
     3552    # We should shift a right by this amount
     3553    cdef Py_ssize_t shift = sa - sb - 54
     3554
     3555    # At this point, we know that q0 = a/b / 2^shift satisfies
     3556    # 2^53 < q0 < 2^55.
     3557    # The end result d = q0 * 2^shift (rounded).
     3558
     3559    # Check for obvious overflow/underflow before shifting
     3560    if shift <= -1130:  # |d| < 2^-1075
     3561        if resultsign < 0:
     3562            return -0.0
     3563        else:
     3564            return 0.0
     3565    elif shift >= 971:  # |d| > 2^1024
     3566        if resultsign < 0:
     3567            return -1.0/0.0
     3568        else:
     3569            return 1.0/0.0
     3570
     3571    sig_on()
     3572
     3573    # Compute q = trunc(a / 2^shift) and let remainder_is_zero be True
     3574    # if and only if no truncation occured.
     3575    cdef mpz_t q, r
     3576    mpz_init(q)
     3577    mpz_init(r)
     3578    cdef bint remainder_is_zero
     3579    if shift >= 0:
     3580        mpz_tdiv_r_2exp(r, mpq_numref(x), shift)
     3581        remainder_is_zero = (mpz_cmp_ui(r, 0) == 0)
     3582        mpz_tdiv_q_2exp(q, mpq_numref(x), shift)
     3583    else:
     3584        mpz_mul_2exp(q, mpq_numref(x), -shift)
     3585        remainder_is_zero = True
     3586
     3587    # Now divide by b to get q = trunc(a/b / 2^shift).
     3588    # remainder_is_zero is True if and only if no truncation occured
     3589    # (in neither division).
     3590    mpz_tdiv_qr(q, r, q, mpq_denref(x))
     3591    if remainder_is_zero:
     3592        remainder_is_zero = (mpz_cmp_ui(r, 0) == 0)
     3593   
     3594    # The quotient q now has 54 or 55 bits, but we need exactly 54.
     3595    # If it's too large, shift it down.
     3596    cdef Py_ssize_t add_shift = mpz_sizeinbase(q, 2) - 54
     3597
     3598    if (shift + add_shift) < -1075:
     3599        # The result will be denormal, ensure the final shift is -1075
     3600        # to avoid a double rounding.
     3601        add_shift = -1075 - shift
     3602
     3603    # Add add_shift to shift and let q = trunc(a/b / 2^shift)
     3604    # for the new shift value.
     3605    if add_shift:
     3606        assert add_shift > 0
     3607        shift += add_shift
     3608        # We do an additional division of q by 2^add_shift.
     3609        if remainder_is_zero:
     3610            mpz_tdiv_r_2exp(r, q, add_shift)
     3611            remainder_is_zero = (mpz_cmp_ui(r, 0) == 0)
     3612        mpz_tdiv_q_2exp(q, q, add_shift)
     3613    mpz_clear(r)
     3614
     3615    if (mpz_tstbit(q, 0) == 0):
     3616        # Round towards zero
     3617        pass
     3618    else:
     3619        if not remainder_is_zero:
     3620            # Remainder is non-zero: round away from zero
     3621            if mpz_sgn(q) > 0:
     3622                mpz_add_ui(q, q, 1)
     3623            else:
     3624                mpz_sub_ui(q, q, 1)
     3625        else:
     3626            # Mid case: round to even
     3627            if mpz_tstbit(q, 1):
     3628                mpz_add_ui(q, q, 1)
     3629            else:
     3630                mpz_sub_ui(q, q, 1)
     3631    sig_off()
     3632
     3633    # Special case zero
     3634    if mpz_cmp_ui(q, 0) == 0:
     3635        mpz_clear(q)
     3636        if resultsign < 0:
     3637            return -0.0
     3638        else:
     3639            return 0.0
     3640   
     3641    # This conversion to double is *exact*.
     3642    # This is because q is even and satisfies abs(q) <= 2^54,
     3643    # with 2^53 <= abs(q) <= 2^54 in the normal case.
     3644    cdef double d = mpz_get_d(q)
     3645    mpz_clear(q)
     3646    return ldexp(d, shift)
     3647
     3648
    34323649def pyrex_rational_reconstruction(integer.Integer a, integer.Integer m):
    34333650    """
    34343651    Find the rational reconstruction of ``a mod m``, if it exists.