So I did some further investigation comparing the ULP error. With the formula that Wilco Dijkstra provided, there are cases where the substitution is super precise. With floats: with input : = 9.99999940395355224609375000000000000000000000000000e-01 sinh: before: = 2.89631005859375000000000000000000000000000000000000e+03 sinh: after : = 2.89630932617187500000000000000000000000000000000000e+03 sinh: mpfr : = 2.89630924626497842670468162463283783344599446025119e+03 ulp err befr: = 3 ulp err aftr: = 0
With doubles: with input : = 9.99999999999999888977697537484345957636833190917969e-01 sinh: before: = 6.71088640000000298023223876953125000000000000000000e+07 sinh: after : = 6.71088639999999925494194030761718750000000000000000e+07 sinh: mpfr : = 6.71088639999999944120645523071287770030292885894208e+07 ulp err befr: = 3 ulp err aftr: = 0 *However*, there are cases where some error shows up. The biggest ULP error that I could find was 2. With floats: with input : = 9.99968349933624267578125000000000000000000000000000e-01 sinh: before: = 1.25686134338378906250000000000000000000000000000000e+02 sinh: after : = 1.25686149597167968750000000000000000000000000000000e+02 sinh: mpfr : = 1.25686137592274042266452526368087062890399889097864e+02 ulp err befr: = 0 ulp err aftr: = 2 With doubles: with input : = 9.99999999999463651256803586875321343541145324707031e-01 sinh: before: = 9.65520209507428342476487159729003906250000000000000e+05 sinh: after : = 9.65520209507428109645843505859375000000000000000000e+05 sinh: mpfr : = 9.65520209507428288553227922831618987450806468855883e+05 ulp err befr: = 0 ulp err aftr: = 2 And with FMA we have the same results showed above. (super precise cases, and maximum ULP error equal 2). So maybe update the patch with the following rules? * If FMA is available, then compute 1 - x*x with it. * If FMA is not available, then do the dijkstra substitution when |x| > 0.5. The code I used for testing: https://pastebin.com/zxPeXmJB On Fri, Oct 19, 2018 at 11:32 AM Wilco Dijkstra <wilco.dijks...@arm.com> wrote: > > Hi, > > >> Maybe I am crazy, or the labels here are wrong, but that looks like the > >> error is three times as *big* after the patch. I.e. it worsened instead > >> of improving. > > This error is actually 1ULP, so just a rounding error. Can't expect any > better than that! > > > with input : = 9.99999880790710449218750000000000000000000000000000e-01 > > cosh: before: = 2.04800000000000000000000000000000000000000000000000e+03 > > cosh: after : = 2.04800024414062500000000000000000000000000000000000e+03 > > cosh: mpfr : = 2.04800006103515897848424084406334262726138617589463e+03 > > error before: = 6.10351589784842408440633426272613861758946325324235e-05 > > error after : = 1.83105466021515759155936657372738613824105367467577e-04 > > It may be less confusing to print relative error or ULP error... > > Wilco >
