On 8/23/19 5:19 AM, Elen Kalda wrote: > Hi all, > > Libgcc currently implements the Smith's algorithm > (https://dl.acm.org/citation.cfm?id=368661) for complex division. It is > designed to prevent overflows and underflows and it does that well for around > 97.7% cases. However, there has been some interest in improving that result, > as > an example Baudin improves the robustness to ~99.5% by introducing an extra > underflow check and also gains in performance by replacing two divisions by > two > multiplications and one division (https://arxiv.org/abs/1210.4539). > > This patch proposes the implementation of Baudin algorithm to Libgcc (in case > the increase in rounding error is acceptable). > > In general, when it comes to deciding what is the best way of doing the > complex divisions, there are three things to consider: > - speed > - accuracy > - robustness (how often overflow/underflow happens, how many division results > are way off from the exact value) > > Improving the algorithms by looking at the failures of specific cases is not > simple. There are 4 real numbers to be manipulated to get the final result and > some minimum number of operations that needs to be done. The quite minimal > Smiths algorithm currently does 9 operations, of which 6 are multiplications > or divisions which are susceptible to overflow and underflow. It is not > feasible to check for both, underflow and overflow, in all 6 operations > without significantly impacting the performance. Most of the complex division > algorithms have been created to solve some special difficult cases, however, > because of the abundance of failure opportunities, algorithm that works well > for some type of divisions is likely to fail for other types of divisions. > > The simplest way to dodge overflow and underflow without impacting performance > much is to vary the order of operations. When naively expanding the complex > division x/y = (...)/(real(y)*real(y) + imag(y)*imag(y)), the squaring in the > denominator is the greatest source of overflow and underflow. Smiths > improvement was to introduce scaling factor r = real(y)/imag(y), such that the > denominator becomes real(y) + imag(y)*r. Baudin's improvement was to check > explicitly for underflow in r and change the order of operations if necessary. > > One way of comparing the algorithms is to generate lots of random complex > numbers and see how the different algorithms perform. That approach is closer > to the "real world" situation where the complex numbers are often random, not > powers of two (which oddly is the assumption many authors make, including > Baudin) and have an exponents in a less "dangerous" ranges. Smiths algorithm > was compared to two versions of Baudin. The difference is that in one version > (b2div), real and imaginary parts of the result are both divided through by > the same denominator d, in the other one (b1div) the real and imaginary > parts are multiplied through by the reciprocal of that d, effectively > replacing two divisions with one division and two multiplications. Since > division is significantly slower on AArch64, that swap introduces gains in > performance (though with the cost of rounding error, unfortunately). > > To compare the performance, I used a test that generates 1800 random complex > double pairs (which fit nicely into the 64 kb cache) and divide each pair 200 > 000 times, effectively doing 360 000 000 operations. > > | CPU time > ------------------ > smiths | 7 296 924 > b1div | 6 558 148 > b2div | 8 658 079 > > On AArch64, b1div is 10.13% faster than Smiths, b2div is 18.65% slower than > Smiths. > > For the accuracy, 1 000 000 pairs of complex doubles were divided and compared > to the exact results (assuming that the division of the same numbers as long > doubles gives the exact value). > > | >2ulp | >1ulp | >0.5ulp | <0.5ulp > --------------------------------------------- > smiths | 22 937 | 23 944 | 124 228 | 828 891 > b1div | 4 450 | 59 716 | 350 792 | 585 042 > b2div | 4 036 | 24 488 | 127 144 | 844 332 > > Same data, but showing the percentage of divisions falling into each > category: > > | >2ulp | >1ulp | >0.5ulp | <0.5ulp > --------------------------------------------- > smiths | 2.29% | 2.39% | 12.42% | 82.89% > b1div | 0.45% | 5.97% | 35.08% | 58.50% > b2div | 0.40% | 2.45% | 12.71% | 84.43% > > The rightmost column (<0.5ulp) shows the number of calculations for which the > ulp error of both, the real and imaginary parts, were inside 0.5 ulp, meaning > they were correctly rounded. >0.5ulp gives the number of calculations for > which the largest ulp error of the real and imaginary parts were between 0.5 > and 1.0 ulp, meaning that it was rounded to one of the nearest doubles, but > not THE nearest double. Anything less accurate is not great... > > FMA doesn't create any problems on AArch64. Bootstrapped and tested on > aarch64-none-linux-gnu -- no problems with bootstrap, but some regressions in > gcc/testsuite/gcc.dg/torture/builtin-math-7.c. The regressions are a result of > the loss in accuracy - for a division (3. + 4.i) / 5i = 0.8 - 0.6i, Smiths > gives: > 0.800000000000000044408920985006 - 0.599999999999999977795539507497 i > ulp error in real part: 0.4 > ulp error in imaginary part: 0.2 > b1div gives: > 0.800000000000000044408920985006 - 0.600000000000000088817841970013 i > ulp error in real part: 0.4 > ulp error in imaginary part: 0.8 > That means the imaginary part of the Baudin result is rounded to one of the > two > nearest floats but not the one which is the nearest. Similar thing happens to > another division in that test. > > > Some questions: > - Is the 10.13% improvement in performance with b1div worth the loss in > accuracy? > - In the case of b1div, most of the results that ceased to be correctly > rounded as a result of replacing the two divisions with multiplications, ended > up in being in 1.0ulp. Maybe that is acceptable? > - The improved algorithm reduces the number of bad fails from 2.3% to 0.5%. Is > that a significant/useful improvement? > > > Best wishes, > Elen > > libgcc/ChangeLog: > > 2019-07-31 Elen Kalda <elen.ka...@arm.com> > > * libgcc2.c CONCAT3(__div,MODE,3): Implement the Baudin's algorithm I'd really like to hear from Joseph on this. He's got *far* more experience in evaluating stuff in this space than I do and is almost certainly going to be able to give a more informed opinion than myself.
I believe he's digging out from some time off, so it might be a little while before he can respond. jeff >
