Ticket #14416: 14416_QQ_to_RDF_v2.patch

File 14416_QQ_to_RDF_v2.patch, 12.4 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 c25873d7a2d5df935d780947f80391ad57198129
    # Parent  830fd81e30622731d69ee04882e3ca9c53339a5a
    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
    1025             sage: sv = A.singular_values(eps=None)
    1026             sage: [round(sv[i],15) for i in range(12)]
    1027             [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789,
    1028              0.000233089089022, 1.1163357483e-05, 4.08237611e-07,
    1029              1.1228611e-08, 2.25196e-10, 3.111e-12, 2.6e-14, 0.0]
    1030             sage: (10^-17 < sv[11]) and (sv[11] < 10^-15)
    1031             True
    1032 
    1033             sage: sv = A.singular_values(eps='auto')
    1034             sage: [round(sv[i],15) for i in range(12)]
    1035             [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789,
    1036              0.000233089089022, 1.1163357483e-05, 4.08237611e-07,
    1037              1.1228611e-08, 2.25196e-10, 3.111e-12, 2.6e-14, 0.0]
    1038             sage: sv[11] == 0.0
    1039             True
    1040 
    1041             sage: sv = A.singular_values(eps=1e-4)
    1042             sage: [round(sv[i],15) for i in range(12)]
    1043             [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789,
    1044              0.000233089089022, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    1045             sage: all([sv[i] == 0.0 for i in range(5, 12)])
    1046             True
     1025            sage: A.singular_values(eps=None)  # abs tol 1e-16
     1026            [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789, 0.000233089089022, 1.11633574832e-05, 4.08237611034e-07, 1.12286106622e-08, 2.25196451406e-10, 3.11134627616e-12, 2.65006299814e-14, 1.00050293715e-16]
     1027
     1028            sage: A.singular_values(eps='auto')  # abs tol 1e-16
     1029            [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789, 0.000233089089022, 1.11633574832e-05, 4.08237611034e-07, 1.12286106622e-08, 2.25196451406e-10, 3.11134627616e-12, 2.65006299814e-14, 0.0]
     1030
     1031            sage: A.singular_values(eps=1e-4)
     1032            [1.79537205956, 0.380275245955, 0.0447385487522, 0.00372231223789, 0.000233089089022, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    10471033
    10481034        With Sage's "verbose" facility, you can compactly see the cutoff
    10491035        at work.  In any application of this routine, or those that build upon
     
    19481934        over the same base ring as the coefficient matrix.  ::
    19491935
    19501936            sage: A=matrix(CDF, 5, [1/(i+j+1) for i in range(5) for j in range(5)])
    1951             sage: A.solve_right([1]*5)
    1952             (5.0..., -120.0..., 630.0..., -1120.0..., 630.0...)
     1937            sage: A.solve_right([1]*5)  # rel tol 1e-11
     1938            (5.0, -120.0, 630.0, -1120.0, 630.0)
    19531939
    19541940        TESTS:
    19551941
     
    20882074        over the same base ring as the coefficient matrix.  ::
    20892075
    20902076            sage: A=matrix(CDF, 5, [1/(i+j+1) for i in range(5) for j in range(5)])
    2091             sage: A.solve_left([1]*5)
    2092             (5.0..., -120.0..., 630.0..., -1120.0..., 630.0...)
     2077            sage: A.solve_left([1]*5)  # rel tol 1e-11
     2078            (5.0, -120.0, 630.0, -1120.0, 630.0)
    20932079
    20942080        TESTS:
    20952081
  • sage/matrix/matrix_mod2e_dense.pyx

    diff --git a/sage/matrix/matrix_mod2e_dense.pyx b/sage/matrix/matrix_mod2e_dense.pyx
    a b  
    801801           
    802802            sage: A.randomize(nonzero=False, density=0.05)
    803803            sage: float(A.density())
    804             0.1358539...
     804            0.135854
    805805        """
    806806        if self._ncols == 0 or self._nrows == 0:
    807807            return
  • 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  
    7373
    7474import sage.rings.real_mpfr
    7575import sage.rings.real_double
     76from libc.stdint cimport uint64_t
    7677
    7778cimport sage.rings.fast_arith
    7879import  sage.rings.fast_arith
     
    19992000            -0.23529411764705882
    20002001            sage: float(-4/17)
    20012002            -0.23529411764705882
     2003            sage: float(1/3)
     2004            0.3333333333333333
     2005            sage: float(1/10)
     2006            0.1
     2007
     2008        TESTS:
     2009
     2010        Test that conversion agrees with `RR`::
     2011
     2012            sage: Q = [a/b for a in [-99..99] for b in [1..99]]
     2013            sage: all([RDF(q) == RR(q) for q  in Q])
     2014            True
     2015
     2016        Test that the conversion has correct rounding on simple rationals::
     2017
     2018            sage: for p in [-100..100]:
     2019            ....:   for q in [1..100]:
     2020            ....:       r = RDF(p/q)
     2021            ....:       assert (RR(r).exact_rational() - p/q) <= r.ulp()/2
     2022
     2023        Test larger rationals::
     2024
     2025            sage: Q = continued_fraction(pi, bits=3000).convergents()
     2026            sage: all([RDF(q) == RR(q) for q  in Q])
     2027            True
     2028
     2029        At some point, the continued fraction and direct conversion
     2030        to ``RDF`` should agree::
     2031
     2032            sage: RDFpi = RDF(pi)
     2033            sage: all([RDF(q) == RDFpi for q  in Q[20:]])
     2034            True
    20022035        """
    2003         return mpq_get_d(self.value)
     2036        return mpq_get_d_nearest(self.value)
    20042037
    20052038    def __hash__(self):
    20062039        """
     
    34293462        if neg: v = -v
    34303463        return v
    34313464
     3465
     3466# The except value is just some random double, it doesn't matter what it is.
     3467cdef double mpq_get_d_nearest(mpq_t x) except? -648555075988944.5:
     3468    """
     3469    Convert a ``mpq_t`` to a ``double``, with round-to-nearest-even.
     3470    This differs from ``mpq_get_d()`` which does round-to-zero.
     3471
     3472    TESTS::
     3473
     3474        sage: q= QQ(); float(q)
     3475        0.0
     3476        sage: q = 2^-10000; float(q)
     3477        0.0
     3478        sage: float(-q)
     3479        -0.0
     3480        sage: q = 2^10000/1; float(q)
     3481        inf
     3482        sage: float(-q)
     3483        -inf
     3484
     3485    ::
     3486
     3487        sage: q = 2^-1075; float(q)
     3488        0.0
     3489        sage: float(-q)
     3490        -0.0
     3491        sage: q = 2^52 / 2^1074; float(q)  # Smallest normal double
     3492        2.2250738585072014e-308
     3493        sage: float(-q)
     3494        -2.2250738585072014e-308
     3495        sage: q = (2^52 + 1/2) / 2^1074; float(q)
     3496        2.2250738585072014e-308
     3497        sage: float(-q)
     3498        -2.2250738585072014e-308
     3499        sage: q = (2^52 + 1) / 2^1074; float(q)  # Next normal double
     3500        2.225073858507202e-308
     3501        sage: float(-q)
     3502        -2.225073858507202e-308
     3503        sage: q = (2^52 - 1) / 2^1074; float(q)  # Largest denormal double
     3504        2.225073858507201e-308
     3505        sage: float(-q)
     3506        -2.225073858507201e-308
     3507        sage: q = 1 / 2^1074; float(q)  # Smallest denormal double
     3508        5e-324
     3509        sage: float(-q)
     3510        -5e-324
     3511        sage: q = (1/2) / 2^1074; float(q)
     3512        0.0
     3513        sage: float(-q)
     3514        -0.0
     3515        sage: q = (3/2) / 2^1074; float(q)
     3516        1e-323
     3517        sage: float(-q)
     3518        -1e-323
     3519        sage: q = (2/3) / 2^1074; float(q)
     3520        5e-324
     3521        sage: float(-q)
     3522        -5e-324
     3523        sage: q = (1/3) / 2^1074; float(q)
     3524        0.0
     3525        sage: float(-q)
     3526        -0.0
     3527        sage: q = (2^53 - 1) * 2^971/1; float(q)  # Largest double
     3528        1.7976931348623157e+308
     3529        sage: float(-q)
     3530        -1.7976931348623157e+308
     3531        sage: q = (2^53) * 2^971/1; float(q)
     3532        inf
     3533        sage: float(-q)
     3534        -inf
     3535        sage: q = (2^53 - 1/2) * 2^971/1; float(q)
     3536        inf
     3537        sage: float(-q)
     3538        -inf
     3539        sage: q = (2^53 - 2/3) * 2^971/1; float(q)
     3540        1.7976931348623157e+308
     3541        sage: float(-q)
     3542        -1.7976931348623157e+308
     3543
     3544    AUTHORS:
     3545
     3546    - Paul Zimmermann, Jeroen Demeyer (:trac:`14416`)
     3547    """
     3548    cdef __mpz_struct* a = <__mpz_struct*>(mpq_numref(x))
     3549    cdef __mpz_struct* b = <__mpz_struct*>(mpq_denref(x))
     3550    cdef int resultsign = mpz_sgn(a)
     3551
     3552    if resultsign == 0:
     3553        return 0.0
     3554
     3555    cdef Py_ssize_t sa = mpz_sizeinbase(a, 2)
     3556    cdef Py_ssize_t sb = mpz_sizeinbase(b, 2)
     3557   
     3558    # Easy case: both numerator and denominator are exactly
     3559    # representable as doubles.
     3560    if sa <= 53 and sb <= 53:
     3561        return mpz_get_d(a) / mpz_get_d(b)
     3562   
     3563    # General case
     3564   
     3565    # We should shift a right by this amount
     3566    cdef Py_ssize_t shift = sa - sb - 54
     3567
     3568    # At this point, we know that q0 = a/b / 2^shift satisfies
     3569    # 2^53 < q0 < 2^55.
     3570    # The end result d = q0 * 2^shift (rounded).
     3571
     3572    # Check for obvious overflow/underflow before shifting
     3573    if shift <= -1130:  # |d| < 2^-1075
     3574        if resultsign < 0:
     3575            return -0.0
     3576        else:
     3577            return 0.0
     3578    elif shift >= 971:  # |d| > 2^1024
     3579        if resultsign < 0:
     3580            return -1.0/0.0
     3581        else:
     3582            return 1.0/0.0
     3583
     3584    sig_on()
     3585
     3586    # Compute q = trunc(a / 2^shift) and let remainder_is_zero be True
     3587    # if and only if no truncation occurred.
     3588    cdef mpz_t q, r
     3589    mpz_init(q)
     3590    mpz_init(r)
     3591    cdef bint remainder_is_zero
     3592    if shift > 0:
     3593        mpz_tdiv_r_2exp(r, a, shift)
     3594        remainder_is_zero = (mpz_cmp_ui(r, 0) == 0)
     3595        mpz_tdiv_q_2exp(q, a, shift)
     3596    else:
     3597        mpz_mul_2exp(q, a, -shift)
     3598        remainder_is_zero = True
     3599
     3600    # Now divide by b to get q = trunc(a/b / 2^shift).
     3601    # remainder_is_zero is True if and only if no truncation occurred
     3602    # (in neither division).
     3603    mpz_tdiv_qr(q, r, q, b)
     3604    if remainder_is_zero:
     3605        remainder_is_zero = (mpz_cmp_ui(r, 0) == 0)
     3606
     3607    # Convert q to a 64-bit integer.
     3608    cdef mp_limb_t* q_limbs = (<__mpz_struct*>q)._mp_d
     3609    cdef uint64_t q64
     3610    if sizeof(mp_limb_t) >= 8:
     3611        q64 = q_limbs[0]
     3612    else:
     3613        assert sizeof(mp_limb_t) == 4
     3614        q64 = q_limbs[1]
     3615        q64 = (q64 << 32) + q_limbs[0]
     3616
     3617    mpz_clear(q)
     3618    mpz_clear(r)
     3619    sig_off()
     3620   
     3621    # The quotient q64 has 54 or 55 bits, but we need exactly 54.
     3622    # Shift it down by 1 one if needed.
     3623    cdef Py_ssize_t add_shift
     3624    if q64 < (1ULL << 54):
     3625        add_shift = 0
     3626    else:
     3627        add_shift = 1
     3628
     3629    if (shift + add_shift) < -1075:
     3630        # The result will be denormal, ensure the final shift is -1075
     3631        # to avoid a double rounding.
     3632        add_shift = -1075 - shift
     3633
     3634    # Add add_shift to shift and let q = trunc(a/b / 2^shift)
     3635    # for the new shift value.
     3636    cdef uint64_t mask
     3637    if add_shift:
     3638        assert add_shift > 0
     3639        assert add_shift < 64
     3640        shift += add_shift
     3641        # We do an additional division of q by 2^add_shift.
     3642        if remainder_is_zero:
     3643            mask = ((1ULL << add_shift)-1)
     3644            remainder_is_zero = ((q64 & mask) == 0)
     3645        q64 = q64 >> add_shift
     3646
     3647    if ((q64 & 1) == 0):
     3648        # Round towards zero
     3649        pass
     3650    else:
     3651        if not remainder_is_zero:
     3652            # Remainder is non-zero: round away from zero
     3653            q64 += 1
     3654        else:
     3655            q64 += (q64 & 2) - 1
     3656
     3657    # The conversion of q64 to double is *exact*.
     3658    # This is because q64 is even and satisfies q64 <= 2^54,
     3659    # (with 2^53 <= q64 <= 2^54 unless in the denormal case).
     3660    cdef double d = <double>q64
     3661    if resultsign < 0:
     3662        d = -d
     3663    return ldexp(d, shift)
     3664
     3665
    34323666def pyrex_rational_reconstruction(integer.Integer a, integer.Integer m):
    34333667    """
    34343668    Find the rational reconstruction of ``a mod m``, if it exists.