cvs commit: src/lib/msun/src s_cbrtf.c

[Bruce Evans]

bde         2005-12-11 13:22:01 UTC

  FreeBSD src repository

  Modified files:
    lib/msun/src         s_cbrtf.c 
  Log:
  Fixed all 502518670 errors of more than 1 ulp for cbrtf() on amd64.
  The maximum error was 3.56 ulps.
  
  The bug was another translation error.  The double precision version
  has a comment saying "new cbrt to 23 bits, may be implemented in
  precision".  This means exactly what it says -- that the 23 bit second
  approximation for the double precision cbrt() may be implemented in
  single (i.e., float) precision.  It doesn't mean what the translation
  assumed -- that this approximation, when implemented in float precision,
  is good enough for the the final approximation in float precision.
  First, float precision needs a 24 bit approximation.  The "23 bit"
  approximation is actually good to 24 bits on float precision args, but
  only if it is evaluated in double precision.  Second, the algorithm
  requires a cleanup step to ensure its error bound.
  
  In float precision, any reasonable algorithm works for the cleanup
  step.  Use the same algorithm as for double precision, although this
  is much more than enough and is a significant pessimization, and don't
  optimize or simplify anything using double precision to implement the
  float case, so that the whole double precision algorithm can be verified
  in float precision.  A maximum error of 0.667 ulps is claimed for cbrt()
  and the max for cbrtf() using the same algorithm shouldn't be different,
  but the actual max for cbrtf() on amd64 is now 0.9834 ulps.  (On i386
  -O1 the max is 0.5006 (down from < 0.7) due to extra precision.)
  
  Revision  Changes    Path
  1.8       +13 -1     src/lib/msun/src/s_cbrtf.c