1
$\begingroup$

I have been using various Krylov subspace methods for a long time, especially conjugate gradient. I'm wondering if there is some way to modify the algorithms to postpone / avoid the divisions that occur so that one could build the whole algorithm using only multiplications and additions. Any reference or description of how to do it would be welcome.

$\endgroup$
  • 1
    $\begingroup$ What is the motivation behind the question? On some system platforms, floating-point division is basically implemented via floating-point multiplies and adds (often, using fused multiply-add operations, or FMA), and one could make that replacement on other platforms with a native division capability (possibly at some loss of performance). $\endgroup$ – njuffa Oct 11 '17 at 4:28
  • $\begingroup$ @njuffa : Yes the idea would be to try and reduce the latency and pipeline clogging in modern CPUs when doing a hardware DIV (which as far as I know) is using the euclidean algorithm. What you describe sounds very interesting. Feel free to refer me to any place where describes these division algorithms. $\endgroup$ – mathreadler Oct 11 '17 at 6:34
  • 1
    $\begingroup$ See Wikipedia as a starter. I worked on the AMD Athlon processor, whose FPU uses Goldschmidt's algorithm for division. I later worked on NVIDIA's parallel programming technology CUDA, which implements division in software using an iteration for the reciprocal with cubic convergence (Halley iteration). Earlier work in this direction was done for Intel's Itanium processor line (see publications by Peter Markstein). It's fairly easy to implement this almost correctly rounded, challenging to do with correct IEEE-754 rounding $\endgroup$ – njuffa Oct 11 '17 at 7:38
  • $\begingroup$ That's funny, I was just gonna ask a question about something very similar to Goldschmidt, it would be very fast now that we have the FFS instructions since a few years back. So it is actually used? I was under the impression that euclidean algorithm was still being used mostly. $\endgroup$ – mathreadler Oct 11 '17 at 7:47
  • $\begingroup$ When you talk about Euclidean algorithm and FFS it seems to suggest you are asking about integer division, whereas I originally assumed you are interested in floating-point division (context of conjugate gradient methods). Could you clarify the question please? I was just in the process of pulling together a few reference for floating-point division based on multiply/add. $\endgroup$ – njuffa Oct 11 '17 at 7:52
1
$\begingroup$

As noted by Wikipedia the primary algorithms for fast floating-point division algorithms based on existing floating-point multiply and add capabilities are Newton-Raphson iteration and the related Goldschmidt algorithm. The latter usually provides more parallelism, but is not self-correcting like NR-iterations , making error analysis more difficult. It is relatively easy to build implementations that deliver approximate results, but designing implementations that provide proper rounding as defined by the floating-point standard IEEE-754 is a task for experts.

Various hardware platforms have no native floating-point division capability and use sequences of multiplies and adds, or – more commonly these days – sequences of fused multiply-add operations (FMA). The AMD Athlon processor used Goldschmidt's algorithm, the Intel Itanium processor used NR iterations, and NVIDIA GPUs using CUDA use Halley iteration on the reciprocal (related to NR-iterations but with cubic convergence). This latter work is not documented in a paper, but one can easily examine the relevant code through disassembly of a CUDA binary containing floating-point division with cuobjdump --dump-sass.

Stuart F. Oberman, "Floating Point Division and Square Root Algorithms and Implementation in the AMD-K7TM Microprocessor". In Proceedings 14th IEEE Symposium on Computer Arithmetic, IEEE 1999, pp. 106-115 (online)

Peter Markstein, "Software division and square root using Goldschmidt’s algorithms." Proceedings of the 6th Conference on Real Numbers and Computers (RNC’6), Vol. 123, Nov. 2004, pp. 146-157 (online)

Peter Markstein, IA-64 and Elementary Functions. Speed and Precision, Prentice Hall 2000 (see chapter 8)

Marius Cornea, John Harrison, and Ping Tak Peter Tang, Scientific Computing on Itanium®-based Systems. Intel Press 2002 (see chapter 8)

Marius Cornea, "Software implementations of division and square root operations for Intel® Itanium® processors." In Proceedings of the 2004 workshop on Computer architecture education, ACM 2004. Article No. 17 (online)


Integer division can be implemented in much the same way, see this question on Stackoverflow for a worked example of 64-bit integer division via Halley iteration.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.