On Wed, Nov 20, 2013 at 09:13:14PM +0100, FX wrote: > > There is a missed optimization in > > > > + if (x < 12) > > + { > > + /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2). > > + This is not perfect, but much better than netlib. */ > > + return erfcq(x) * expq(x*x); > > + } > > > > If x is less than approximately -8192, then erfc(x)*exp(x*x) > > overflows. > > Hum, I get roughly -106 where you have -8192, so I?ll not > commit immediately with the value I have, and let you check it first. > > OK to commit? >
The optimization can come later. Your current version is fine, because it fixes the ULP issue. As to the optimization, for large negative x, one has efrc(x)*exp(x*x) --> 2 * exp(x*x) = 2**emax exp(x*x) = 2**(emax-1) x*x = (emax-1)*log(2) x = - sqrt((emax-1)*log(2)) ! Choosing minus branch x = -106.56... ! emax=16384 (could be off-by-1, did not check) So, yeah, you're correct. My suggestion was based on the not so careful mistake of replacing x*x by x+x and dropping log(2). That is, I had x+x = -emax --> x = - emax / 2. -- Steve