cvs commit: src/lib/msun/src e_rem_pio2.c e_rem_pio2f.c

[Bruce Evans]

bde         2008-02-28 16:22:36 UTC

  FreeBSD src repository

  Modified files:
    lib/msun/src         e_rem_pio2.c e_rem_pio2f.c 
  Log:
  Fix and improve some magic numbers for the "medium size" case.
  
  e_rem_pio2.c:
  This case goes up to about 2**20pi/2, but the comment about it said that
  it goes up to about 2**19pi/2.
  
  It went too far above 2**pi/2, giving a multiplier fn with 21 significant
  bits in some cases.  This would be harmful except for a numerical
  accident.  It happens that the terms of the approximation to pi/2,
  when rounded to 33 bits so that multiplications by 20-bit fn's are
  exact, happen to be rounded to 32 bits so multiplications by 21-bit
  fn's are exact too, so the bug only complicates the error analysis (we
  might lose a bit of accuracy but have bits to spare).
  
  e_rem_pio2f.c:
  The bogus comment in e_rem_pio2.c was copied and the code was changed
  to be bug-for-bug compatible with it, except the limit was made 90
  ulps smaller than necessary.  The approximation to pi/2 was not
  modified except for discarding some of it.
  
  The same rough error analysis that justifies the limit of 2**20pi/2
  for double precision only justifies a limit of 2**18pi/2 for float
  precision.  We depended on exhaustive testing to check the magic numbers
  for float precision.  More exaustive testing shows that we can go up
  to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision,
  with a the maximum error for cosf() and sinf() unchanged at 0.5009
  ulps despite the maximum error in rem_pio2f being ~0.25 ulps.  Implement
  this.
  
  Revision  Changes    Path
  1.18      +1 -1      src/lib/msun/src/e_rem_pio2.c
  1.29      +3 -3      src/lib/msun/src/e_rem_pio2f.c