Ticket #14082: 14082_pari_complex_double.patch

File 14082_pari_complex_double.patch, 12.8 KB (added by jdemeyer, 7 years ago)
  • sage/rings/complex_double.pxd

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1360321201 -3600
    # Node ID a47af7899aa7979ae6534397bfdc8736e209a2cb
    # Parent  107c810587f2667c634e921deeba9e8376e83ffb
    Fix PARI calls for ComplexDoubleElement
    
    diff --git a/sage/rings/complex_double.pxd b/sage/rings/complex_double.pxd
    a b  
    1717        cdef gsl_complex _complex
    1818        cdef GEN _gen(self)
    1919        cdef ComplexDoubleElement _new_c(self, gsl_complex x)
    20         cdef _new_from_gen_c(self, GEN g, pari_sp sp)
     20        cdef _new_from_gen(self, sage.libs.pari.gen.gen g)
    2121
    2222cdef ComplexDoubleElement new_ComplexDoubleElement()
  • sage/rings/complex_double.pyx

    diff --git a/sage/rings/complex_double.pyx b/sage/rings/complex_double.pyx
    a b  
    4848- William Stein (2006-09): first version
    4949
    5050- Travis Scrimshaw (2012-10-18): Added doctests to get full coverage
     51
     52- Jeroen Demeyer (2013-02-27): fixed all PARI calls (:trac:`14082`)
    5153"""
    5254
    53 #################################################################################
     55#*****************************************************************************
    5456#       Copyright (C) 2006 William Stein <wstein@gmail.com>
     57#       Copyright (C) 2013 Jeroen Demeyer <jdemeyer@cage.ugent.be>
    5558#
    5659#  Distributed under the terms of the GNU General Public License (GPL)
    57 #
     60#  as published by the Free Software Foundation; either version 2 of
     61#  the License, or (at your option) any later version.
    5862#                  http://www.gnu.org/licenses/
    5963#*****************************************************************************
    6064
     
    6569include '../ext/interrupt.pxi'
    6670include '../ext/stdsage.pxi'
    6771
    68 cdef extern from "math.h":
    69     double modf (double value, double *integer_part)
    70 
    71 # The M_PI_4 constant is not available on cygwin in "math.h" (though
    72 # it is on most other platforms).
    73 cdef double M_PI_4 = 0.785398163397448309615660845819875721
    74 
    75 cdef extern from "complex.h":
     72cdef extern from "<complex.h>":
    7673    double complex csqrt(double complex)
    7774    double cabs(double complex)
    7875
     
    9895
    9996from real_double import RealDoubleElement, RDF
    10097
    101 # PREC is the precision (in decimal digits) that all PARI computations with doubles
    102 # are done with in this module.  A double is by definition 8 bytes or 64 bits.  Since
    103 #              log(2^64,10) = 19.265919722494793  < 28
    104 # 28 *digits* should be way more than enough for correctness to 64 bits.   We choose
    105 # 28 since 20 gets rounded up to 28 digits automatically by PARI anyways.   In theory
    106 # 19 digits might also be OK, but I've noticed that last digit or two of PARI output
    107 # can be questionable -- it varies from version to version, so...
    108 cdef int PREC
    109 PREC = 28
    11098
    11199from sage.structure.parent_gens import ParentWithGens
    112100from sage.categories.morphism cimport Morphism
     
    309297            4.5
    310298            sage: CDF(1+I) # indirect doctest
    311299            1.0 + 1.0*I
     300            sage: CDF(pari(1))
     301            1.0
     302            sage: CDF(pari("I"))
     303            1.0*I
     304            sage: CDF(pari("x^2 + x + 1").polroots()[0])
     305            -0.5 - 0.866025403784*I
    312306
    313307        A ``TypeError`` is raised if the coercion doesn't make sense::
    314308
     
    346340            sage: CDF((1,2)) # indirect doctest
    347341            1.0 + 2.0*I
    348342        """
    349         cdef pari_sp sp
     343        cdef sage.libs.pari.gen.gen g
    350344        if PY_TYPE_CHECK(x, ComplexDoubleElement):
    351345            return x
    352346        elif PY_TYPE_CHECK(x, tuple):
     
    358352        elif isinstance(x, complex_number.ComplexNumber):
    359353            return ComplexDoubleElement(x.real(), x.imag())
    360354        elif isinstance(x, sage.libs.pari.gen.gen):
    361             # It seems we should get a speed increase by
    362             # using _new_from_gen_c instead; I wasn't
    363             # able to get this to work.
    364             return ComplexDoubleElement(x.real(), x.imag())
     355            g = x
     356            if typ(g.g) == t_COMPLEX:
     357                return ComplexDoubleElement(gtodouble(gel(g.g, 1)), gtodouble(gel(g.g, 2)))
     358            else:
     359                return ComplexDoubleElement(gtodouble(g.g), 0.0)
    365360        elif isinstance(x, str):
    366361            t = cdf_parser.parse_expression(x)
    367362            if isinstance(t, float):
     
    690685        z._complex = x
    691686        return z
    692687
    693     cdef _new_from_gen_c(self, GEN g, pari_sp sp):
     688    cdef _new_from_gen(self, sage.libs.pari.gen.gen g):
    694689        """
    695690        C-level code for creating a :class:`ComplexDoubleElement` from a
    696691        PARI gen.
    697692
    698693        INPUT:
    699694
    700         -  ``g`` -- GEN
    701 
    702         -  ``sp`` -- stack pointer; if nonzero resets avma to sp.
     695        -  ``g`` -- A PARI ``gen``.
     696
     697        OUTPUT: A ``ComplexDoubleElement``.
    703698        """
    704699        cdef gsl_complex x
    705         sig_on()
    706         # Signal handling is important, since rtodbl can easily overflow
    707         x = gsl_complex_rect(  rtodbl(greal(g)), rtodbl(gimag(g))  )
    708         sig_off()
    709         z = self._new_c(x)
    710         if sp:
    711             avma = sp
    712         return z
     700        if typ(g.g) == t_COMPLEX:
     701            x = gsl_complex_rect(gtodouble(gel(g.g, 1)), gtodouble(gel(g.g, 2)))
     702        else:
     703            x = gsl_complex_rect(gtodouble(g.g), 0.0)
     704        return self._new_c(x)
    713705
    714706    def __hash__(self):
    715707        """
     
    10271019            sage: pari(CDF(1,2))
    10281020            1.00000000000000 + 2.00000000000000*I
    10291021        """
    1030         # The explicit declaration and assignment is necessary in
    1031         # order to get access to the internal C structure of gen.pari.
    1032         cdef pari_sp sp
    1033         global avma
    1034         sp = avma
    10351022        cdef sage.libs.pari.gen.PariInstance P
    10361023        P = sage.libs.pari.gen.pari
    1037         x = P.new_gen_noclear(self._gen())
    1038         avma = sp
    1039         return x
     1024        sig_on()
     1025        return P.new_gen(self._gen())
    10401026
    10411027    #######################################################################
    10421028    # Arithmetic
     
    19421928    #######################################################################
    19431929    def eta(self, int omit_frac=0):
    19441930        r"""
    1945         Return the value of the Dedekind `\eta` function on self,
    1946         intelligently computed using `\mathbb{SL}(2,\ZZ)`
    1947         transformations.
     1931        Return the value of the Dedekind `\eta` function on self.
    19481932
    19491933        INPUT:
    19501934
     
    19561940
    19571941        OUTPUT: a complex double number
    19581942
    1959         ALGORITHM: Uses the PARI C library, but with some modifications so
    1960         it always works instead of failing on easy cases involving large
    1961         input (like PARI does).
     1943        ALGORITHM: Uses the PARI C library.
    19621944
    19631945        The `\eta` function is
    19641946
     
    19791961
    19801962        :meth:`eta()` works even if the inputs are large::
    19811963
    1982             sage: CDF(0,10^15).eta()
     1964            sage: CDF(0, 10^15).eta()
    19831965            0.0
    1984             sage: CDF(10^15,0.1).eta()     # slightly random-ish arch dependent output
    1985             -0.121339721991 - 0.19619461894*I
     1966            sage: CDF(10^15, 0.1).eta()  # abs tol 1e-10
     1967            -0.115342592727 - 0.19977923088*I
    19861968
    19871969        We compute a few values of :meth:`eta()`, but with the fractional power
    19881970        of `e` omitted::
     
    20011983
    20021984        The optional argument allows us to omit the fractional part::
    20031985
    2004             sage: z = CDF(1,1)
    2005             sage: z.eta(omit_frac=True)
    2006             0.998129069926
     1986            sage: z.eta(omit_frac=True)  # abs tol 1e-12
     1987            0.998129069926 - 8.12769318782e-22*I
    20071988            sage: pi = CDF(pi)
    2008             sage: prod([1-exp(2*pi*i*n*z) for n in range(1,10)])      # slightly random-ish arch dependent output
    2009             0.998129069926 + 4.5908467128e-19*I
     1989            sage: prod([1-exp(2*pi*i*n*z) for n in range(1,10)])  # abs tol 1e-12
     1990            0.998129069926 + 4.59084695545e-19*I
    20101991
    20111992        We illustrate what happens when `z` is not in the upper half plane::
    20121993
     
    20232004        """
    20242005        cdef GEN a, b, c, y, t
    20252006
    2026         # Turn on Sage C interrupt handling.  There must
    2027         # be no Python code between here and sig_off().
    2028         #sig_on()  # we're not using this since it would dominate the runtime
    2029 
    20302007        if self._complex.dat[1] <= 0:
    20312008            raise ValueError, "value must be in the upper half plane"
    20322009
    2033         if self._complex.dat[1] > 100000:
     2010        if self._complex.dat[1] > 100000 and not omit_frac:
    20342011            # To the precision of doubles for such large imaginary
    2035             # part the answer is automatically 0.  If we don't do
    2036             # this PARI can easily die, which will cause this function
    2037             # to bomb unless we use sig_on() and sig_off().  But
    2038             # I don't want to use those, since they take more time
    2039             # than this entire function!  Moreover, I don't want Sage's
    2040             # eta to every bomb -- this function should work on all input; there's
    2041             # no excuse for having it fail on easy edge cases (like PARI does).
     2012            # part, the answer is automatically 0. If we don't do
     2013            # this, PARI can easily underflow.
    20422014            return ComplexDoubleElement(0,0)
    20432015
    2044         # save position on PARI stack
    2045         cdef pari_sp sp
    2046         global avma
    2047         sp = avma
    2048 
    2049         # If omit_frac is True, then eta(z) = eta(z+24*n) for any integer n,
    2050         # so we can replace the real part of self by its fractional part.  We
    2051         # do this for two reasons:
    2052         #   (1) Because when the real part is sufficiently large, PARI overflows
    2053         #       and crashes, and I don't want to use sig_on()/sig_off() to catch this.
    2054         #   (2) Because this is easy and it makes it so our eta always works,
    2055         #       unlike PARI's.
    2056         # convert our complex number to a PARI number
    2057         # If omit_frac is False, then eta(z) = eta(z+24*n) for any integer n,
    2058         # and similar remarks apply.
    2059         cdef double dummy
    2060         a = dbltor(modf(self._complex.dat[0], &dummy))    # see big comment above
    2061         b = dbltor(self._complex.dat[1])
    2062 
    2063         y = cgetg(3, t_COMPLEX)    # allocate space for a complex number
    2064         set_gel(y,1,a); set_gel(y,2,b)
    2065 
    2066         c = eta(y, PREC)
    2067 
    2068         # Convert the PARI complex number c back to a GSL complex number
    2069         cdef gsl_complex w
    2070         w = gsl_complex_rect(rtodbl(greal(c)), rtodbl(gimag(c)))
    2071 
    2072         # put the pari stack back; this frees all memory used
    2073         # by PARI above.
    2074         avma = sp
    2075 
    2076         if not omit_frac:
    2077             # multiply z by e^{\pi i z / 12} = exp(pi * i * z / 12), where z = self.
    2078             w = gsl_complex_mul(w, gsl_complex_exp(gsl_complex_mul(\
    2079                             gsl_complex_rect(0,M_PI_4/(<double>3)), self._complex)))
    2080 
    2081         return self._new_c(w)
     2016        cdef int flag = 0 if omit_frac else 1
     2017        return self._new_from_gen(self._pari_().eta(flag))
    20822018
    20832019    def agm(self, right, algorithm="optimal"):
    20842020        r"""
     
    21402076        cdef double complex a, b, a1, b1, r
    21412077        cdef double d, e, eps = 2.0**-51
    21422078
    2143         cdef pari_sp sp
    2144        
    2145         if algorithm=="pari":
    2146             sp = avma
    2147             return self._new_from_gen_c( agm(self._gen(), complex_gen(right), PREC), sp)
     2079        if algorithm == "pari":
     2080            return self._new_from_gen(self._pari_().agm(right))
    21482081
    21492082        if not isinstance(right, ComplexDoubleElement):
    21502083            right = CDF(right)
     
    21942127            sage: CDF(10000000,10000000).dilog()
    21952128            -134.411774491 + 38.793962999*I
    21962129        """
    2197         cdef pari_sp sp
    2198         sp = avma
    2199         return self._new_from_gen_c(  dilog(self._gen(), PREC),   sp)
     2130        return self._new_from_gen(self._pari_().dilog())
    22002131
    22012132    def gamma(self):
    22022133        r"""
     
    22242155                    return CC(self).gamma()
    22252156            except TypeError:
    22262157                pass
    2227         cdef pari_sp sp
    2228         sp = avma
    2229         return self._new_from_gen_c(  ggamma(self._gen(), PREC),   sp)
     2158        return self._new_from_gen(self._pari_().gamma())
    22302159
    22312160    def gamma_inc(self, t):
    22322161        r"""
     
    22412170            sage: CDF(2,0).gamma_inc(CDF(1,1))
    22422171            0.707092096346 - 0.42035364096*I
    22432172        """
    2244         cdef pari_sp sp
    2245         sp = avma
    2246         sig_on()    # because it is somewhat slow and/or unreliable in corner cases.
    2247         x = self._new_from_gen_c( incgam(self._gen(), complex_gen(t), PREC),   sp )
    2248         sig_off()
    2249         return x
     2173        return self._new_from_gen(self._pari_().incgam(t))
    22502174
    22512175    def zeta(self):
    22522176        """
     
    22652189        if self._complex.dat[0] == 1 and self._complex.dat[1] == 0:
    22662190            import infinity
    22672191            return infinity.unsigned_infinity
    2268         cdef pari_sp sp
    2269         sp = avma
    2270         return self._new_from_gen_c(  gzeta(self._gen(), PREC),   sp)
     2192        return self._new_from_gen(self._pari_().zeta())
    22712193
    22722194    def algdep(self, long n):
    22732195        """