Take the 2-minute tour ×
Programmers Stack Exchange is a question and answer site for professional programmers interested in conceptual questions about software development. It's 100% free, no registration required.

I have a cross platform math library that I am working on and I want to make sure that some common operations are performed in an optimized manner, so I wish to use some intrinsic functions wrapped in inline functions to compute things such as sqrt(float/double). I should mention that at this point I only plan on targeting a select few platforms which are mentioned in the source below.

I have started a single prototype of the a sqrt() function but before I proceed any further I want to know if I am doing anything blatantly wrong.

I am following the documentation at: http://msdn.microsoft.com/en-us/library/windows/desktop/ff471376%28v=vs.85%29.aspx
http://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html

Also, you will note I have a blank line under the pound define for macintosh. I know the newer versions are using clang/llvm. How do you get access to intrinsic functions under the newer versions of mac os and its associated compilers?

#ifndef INTRINSICS_H_INCLUDED
#define INTRINSICS_H_INCLUDED

#if defined(_WIN32)
#elif defined(__gnu_linux__) || defined(__linux__)
#elif defined(__ANDROID__)
#elif defined(__CYGWIN__)
#elif defined(macintosh)
#else
#include <cmath>
#endif

#if defined(_WIN32)
#pragma function( sqrt )
#endif
template<typename T>
inline float sqrt(T &val)
{
#if defined(_WIN32)
  return sqrt(val);
#elif defined(__gnu_linux__) || defined(__linux__)
  return __builtin_sqrt(val);
#elif defined(__ANDROID__)
  return __builtin_sqrt(val);
#elif defined(__CYGWIN__)
  return __builtin_sqrt(val);
#elif defined(macintosh)
  // How do I do this
#else
  return std::sqrt(val);
#endif
}
#if defined(_WIN32)
#pragma function( sqrt )
#endif
#endif
share|improve this question
    
I would expect that <cmath> would include proper definitions in the first place. First check whether id does by looking at the generated optimized assembly. As for using the __builtin_* function, you should simply use them #if defined(__GNUC__) as they are not platform-specific but compiler-specific. –  Jan Hudec May 30 '13 at 6:07
6  
Unless you have concrete proof (as in, actual measurements) that your wrapper function around the intrinsics performs better than just using std::sqrt directly, you are committing the sin of premature optimisation. Where intrinsics are available, I would expect std::sqrt to be just as thin a wrapper around them. –  Bart van Ingen Schenau May 30 '13 at 6:52
1  
I agree indeed. It is a sin of premature optimization. And it seems to be a not invented here syndrome too. That is, is it really necessary you do this yourself? There are many math libraries, well tested and carefully tuned. What are you doing so special that it does not exist yet? –  Jan Hudec May 30 '13 at 7:36
    
    
My apologies because I originally answered that the simple approach should be taken (just using std::sqrt), only to be astonished by the difficulty of the topic. The matter is further complicated on ARM processor, where support for sqrt was phased in (starting with pure-software implementation, then provided as a slow CPU instruction, then an alternate algorithm was implemented in NEON, and so on.) This makes the original premise of cross-platform manner hard to achieve. –  rwong May 31 '13 at 7:03
show 1 more comment

2 Answers

up vote 3 down vote accepted

Let me begin with a rhetoric.

If there is a superior implementation that doesn't have any negative trade-off, then why aren't the standard libraries already using it?


It is not a question of premature optimization, if the function is used in a mathematical or large-data processing library where the performance gains will be multiplied by millions or billions. A prime example is the OpenCV library, where performance gains - even if just a few percents here and there - translate into a reciprocal increase in the number of video frames processed per seconds. So let's take a look at the square-root implementation in OpenCV.

https://github.com/Itseez/opencv/blob/master/modules/core/src/mathfuncs.cpp

Summary:

  • If the caller is interested in the reciprocal of square root, go straight for that.
  • If the caller is interested in comparing the magnitudes of two values after taking square root on each one, try to eliminate that square root step (with some additional logic).
  • If the caller is interested in calculating the square roots (or their reciprocals) for an array of values, use the SIMD equivalents.
  • Otherwise, fall back to std::sqrt.

  • (Edited)

    • If the caller only needs the crude approximate of the square root, see the Wikipedia article on this topic. (This topic is so huge that it's not possible to summarize here.)

Edit:

I stand corrected (by dan04's answer) by some important omissions from my original answer.

  • Compiler needs to be told the target CPU's architecture.
  • Compiler needs to be told to favor speed over precision
  • Compiler needs to be given permission to omit input range checks
    • Example: sqrt being compiled into a library call
      • In this example, the library function validates that the input value is not negative. If it is, it calls an error handler, which then throws a C++ exception. These additional behaviors makes it impossible to omit the library call.
  • The CPU needs to be told to flush-to-zero, and to treat denormals-as-zeros
    • This is done at runtime.

(I'm giving examples in MSVC, but other compilers would also need to be explicitly told to generate fast code.)

Since these are compiler options, they need to be specified in the makefile, not in the source code.


Suggestion:

It may be necessary to use a makefile configuration utility such as CMake to verify the compiler's support and optimality of the various intrinsics (and for using the fastest version available on that platform).


Mathematical details:

There's a long article on Wikipedia on methods of computing square roots and its rough approximations, including the "Carmack trick" invented by Walsh and Tarolli.

Some of these techniques are suitable for VLSI implementation. If so, it is typical that they are already employed in modern general purpose CPUs. Some of these accelerated functions require the caller to make one or several Newton's method or Goldschmidt iteration in order to converge from an initial estimate.

share|improve this answer
add comment

Instead of cluttering up your source code with all these macros, you could use optimization flags, such as /Oi in MSVC++, to have the compilers do this optimization for you.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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