# HG changeset patch
# User Jeroen Demeyer <jdemeyer@cage.ugent.be>
# Date 1340100027 7200
# Node ID 0a193a60c44d25f1c8e52a09f49752b3f9af5d1f
# Parent 840bb929a987235cfb6ccfd90b290a76f18bfefd
Fix ulp() in RealDoubleField
diff git a/sage/rings/real_double.pyx b/sage/rings/real_double.pyx
a

b


599  599  
600  600  EXAMPLES:: 
601  601  
602   sage: a = RDF(1) 
 602  sage: a = RDF(pi) 
 603  sage: a.ulp() 
 604  4.4408920985e16 
 605  sage: b = a + a.ulp() 
 606  
 607  Adding or subtracting an ulp always gives a different number:: 
 608  
 609  sage: a + a.ulp() == a 
 610  False 
603  611  sage: a  a.ulp() == a 
604  612  False 
605   sage: a  a.ulp()/2 == a 
 613  sage: b + b.ulp() == b 
 614  False 
 615  sage: b  b.ulp() == b 
 616  False 
 617  
 618  Adding or subtracting something less than half an ulp never 
 619  gives the same number (unless the number is exactly a power of 
 620  2 and subtracting an ulp decreases the ulp):: 
 621  
 622  sage: a  a.ulp()/3 == a 
606  623  True 
 624  sage: a + a.ulp()/3 == a 
 625  True 
 626  sage: b  b.ulp()/3 == b 
 627  True 
 628  sage: b + b.ulp()/3 == b 
 629  True 
 630  sage: c = RDF(1) 
 631  sage: c  c.ulp()/3 == c 
 632  False 
 633  sage: c.ulp() 
 634  2.22044604925e16 
 635  sage: (c  c.ulp()).ulp() 
 636  1.11022302463e16 
607  637  
608   sage: a = RDF.pi() 
609   sage: b = a + a.ulp() 
610   sage: (a+b)/2 in [a,b] 
611   True 
 638  The ulp is always positive:: 
 639  
 640  sage: RDF(1).ulp() 
 641  2.22044604925e16 
 642  
 643  The ulp of zero is the smallest positive number in RDF:: 
 644  
 645  sage: RDF(0).ulp() 
 646  4.94065645841e324 
 647  sage: RDF(0).ulp()/2 
 648  0.0 
612  649  
 650  Some special values:: 
 651  
613  652  sage: a = RDF(1)/RDF(0); a 
614  653  +infinity 
615  654  sage: a.ulp() 
… 
… 

620  659  sage: a.ulp() is a 
621  660  True 
622  661  """ 
623   cdef int e, v = gsl_isinf(self._value) 
 662  # First, check special values 
 663  if self._value == 0: 
 664  return RealDoubleElement(ldexp(1.0, 1074)) 
624  665  if gsl_isnan(self._value): 
625  666  return self 
626   elif self._value == 0: 
627   return RealDoubleElement(ldexp(1.0, e1082)) 
628   elif v == 1: 
629   return self 
630   elif v == 1: 
631   return self 
632   else: 
633   frexp(self._value, &e) 
634   return RealDoubleElement(ldexp(1.0, e54)) 
 667  if gsl_isinf(self._value): 
 668  return self.abs() 
 669  
 670  # Normal case 
 671  cdef int e 
 672  frexp(self._value, &e) 
 673  return RealDoubleElement(ldexp(1.0, e53)) 
635  674  
636  675  def real(self): 
637  676  """ 
… 
… 

1731  1770  Make sure that we can take the log of small numbers accurately 
1732  1771  and the fix doesn't break preexisting values (:trac:`12557`):: 
1733  1772  
 1773  sage: R = RealField(128) 
1734  1774  sage: def check_error(x): 
1735  1775  ... x = RDF(x) 
1736  1776  ... log_RDF = x.log() 
1737   ... log_exact = log(x.n(1000)) 
1738   ... return abs(log_RDFlog_exact) <= log_RDF.ulp() 
 1777  ... log_RR = R(x).log() 
 1778  ... diff = R(log_RDF)  log_RR 
 1779  ... if abs(diff) <= log_RDF.ulp(): 
 1780  ... return True 
 1781  ... print "logarithm check failed for %s (diff = %s ulp)"% \ 
 1782  ... (x, diff/log_RDF.ulp()) 
 1783  ... return False 
1739  1784  sage: all( check_error(2^x) for x in range(100,100) ) 
1740  1785  True 
1741  1786  sage: all( check_error(x) for x in sxrange(0.01, 2.00, 0.01) ) 