# HG changeset patch
# User D. S. McNeil <dsm054@gmail.com>
# Date 1337924458 14400
# Node ID 23d577ebaa57001b83288e202c4411b13c173dc2
# Parent 6285dcda57540e4d34f60379f9f87ccd92b3b4ed
Trac #12557: Fix RDF logarithms of small numbers
diff git a/sage/rings/real_double.pyx b/sage/rings/real_double.pyx
a

b


1681  1681  
1682  1682  
1683  1683  cdef _log_base(self, double log_of_base): 
1684   if self._value < 2: 
1685   if self._value == 0: 
1686   return RDF(1)/RDF(0) 
1687   if self._value < 0: 
1688   return RDF(0)/RDF(0) 
1689   sig_on() 
1690   a = self._new_c(gsl_sf_log_1plusx(self._value  1) / log_of_base) 
1691   sig_off() 
1692   return a 
 1684  if self._value == 0: 
 1685  return RDF(1)/RDF(0) 
 1686  elif self._value < 0: 
 1687  return RDF.NaN() 
1693  1688  sig_on() 
1694  1689  a = self._new_c(gsl_sf_log(self._value) / log_of_base) 
1695  1690  sig_off() 
… 
… 

1697  1692  
1698  1693  def log(self, base=None): 
1699  1694  """ 
 1695  Return the logarithm. 
 1696  
 1697  INPUT: 
 1698  
 1699   ``base``  integer or ``None`` (default). The base of the 
 1700  logarithm. If none is specified, the base is `e` (the socalled 
 1701  natural logarithm). 
 1702  
 1703  OUTPUT: 
 1704  
 1705  The logarithm of ``self``. If ``self`` is positive, a double 
 1706  floating point number. Infinity if ``self`` is zero. A 
 1707  imaginary complex floating point number if ``self`` is 
 1708  negative. 
 1709  
1700  1710  EXAMPLES:: 
1701  1711  
1702  1712  sage: RDF(2).log() 
… 
… 

1715  1725  3.14159265359*I 
1716  1726  sage: RDF(1).log(2) 
1717  1727  4.53236014183*I 
 1728  
 1729  TESTS: 
 1730  
 1731  Make sure that we can take the log of small numbers accurately 
 1732  and the fix doesn't break preexisting values (:trac:`12557`):: 
 1733  
 1734  sage: def check_error(x): 
 1735  ... x = RDF(x) 
 1736  ... log_RDF = x.log() 
 1737  ... log_exact = log(x.n(1000)) 
 1738  ... return abs(log_RDFlog_exact) <= log_RDF.ulp() 
 1739  sage: all( check_error(2^x) for x in range(100,100) ) 
 1740  True 
 1741  sage: all( check_error(x) for x in sxrange(0.01, 2.00, 0.01) ) 
 1742  True 
 1743  sage: all( check_error(x) for x in sxrange(0.99, 1.01, 0.001) ) 
 1744  True 
 1745  sage: RDF(1.000000001).log() 
 1746  1.00000008224e09 
 1747  sage: RDF(1e17).log() 
 1748  39.1439465809 
 1749  sage: RDF(1e50).log() 
 1750  115.12925465 
1718  1751  """ 
1719  1752  if self < 0: 
1720  1753  from sage.rings.complex_double import CDF 