Floating point oddity
<p><pre><code>    template &lt;typename T&gt; T q(T x) {
      return fabs(x - sqrt(x*x + 1.)) - 1.&#x2F;(x + sqrt(x*x + 1.));
    }
</code></pre>
The author doesn&#x27;t seem to have a lot of numeric experience. He probably doesn&#x27;t know that, even when T is float, the compiler uses the double version of sqrt.<p>Don&#x27;t believe me? Read the assembler:<p><pre><code>    0000000000000c3c &lt;float q&lt;float&gt;(float)&gt;:
     c3c:   55                      push   rbp
     c3d:   48 89 e5                mov    rbp,rsp
     c40:   48 83 ec 20             sub    rsp,0x20
     c44:   f3 0f 11 45 fc          movss  DWORD PTR [rbp-0x4],xmm0
     c49:   f3 0f 5a 5d fc          cvtss2sd xmm3,DWORD PTR [rbp-0x4]
     c4e:   f2 0f 11 5d f0          movsd  QWORD PTR [rbp-0x10],xmm3
     c53:   f3 0f 10 45 fc          movss  xmm0,DWORD PTR [rbp-0x4]
     c58:   f3 0f 59 45 fc          mulss  xmm0,DWORD PTR [rbp-0x4]
     c5d:   f3 0f 5a c0             cvtss2sd xmm0,xmm0
     c61:   f2 0f 10 0d 4f 02 00    movsd  xmm1,QWORD PTR [rip+0x24f]        # eb8 &lt;std::piecewise_construct+0x8&gt;
     c68:   00
     c69:   f2 0f 58 c1             addsd  xmm0,xmm1
     c6d:   e8 5e fc ff ff          call   8d0 &lt;sqrt@plt&gt;
     c72:   f2 0f 10 5d f0          movsd  xmm3,QWORD PTR [rbp-0x10]
     c77:   f2 0f 5c d8             subsd  xmm3,xmm0
     c7b:   66 0f 28 c3             movapd xmm0,xmm3
     c7f:   f3 0f 7e 0d 39 02 00    movq   xmm1,QWORD PTR [rip+0x239]        # ec0 &lt;std::piecewise_construct+0x10&gt;
     c86:   00
     c87:   66 0f 54 c1             andpd  xmm0,xmm1
     c8b:   f2 0f 11 45 f0          movsd  QWORD PTR [rbp-0x10],xmm0
     c90:   f3 0f 5a 65 fc          cvtss2sd xmm4,DWORD PTR [rbp-0x4]
     c95:   f2 0f 11 65 e8          movsd  QWORD PTR [rbp-0x18],xmm4
     c9a:   f3 0f 10 45 fc          movss  xmm0,DWORD PTR [rbp-0x4]
     c9f:   f3 0f 59 45 fc          mulss  xmm0,DWORD PTR [rbp-0x4]
     ca4:   f3 0f 5a c0             cvtss2sd xmm0,xmm0
     ca8:   f2 0f 10 0d 08 02 00    movsd  xmm1,QWORD PTR [rip+0x208]        # eb8 &lt;std::piecewise_construct+0x8&gt;
     caf:   00
     cb0:   f2 0f 58 c1             addsd  xmm0,xmm1
     cb4:   e8 17 fc ff ff          call   8d0 &lt;sqrt@plt&gt;
     cb9:   f2 0f 58 45 e8          addsd  xmm0,QWORD PTR [rbp-0x18]
     cbe:   f2 0f 10 0d f2 01 00    movsd  xmm1,QWORD PTR [rip+0x1f2]        # eb8 &lt;std::piecewise_construct+0x8&gt;
     cc5:   00
     cc6:   f2 0f 5e c8             divsd  xmm1,xmm0
     cca:   66 0f 28 c1             movapd xmm0,xmm1
     cce:   f2 0f 10 55 f0          movsd  xmm2,QWORD PTR [rbp-0x10]
     cd3:   f2 0f 5c d0             subsd  xmm2,xmm0
     cd7:   66 0f 28 c2             movapd xmm0,xmm2
     cdb:   f2 0f 5a c0             cvtsd2ss xmm0,xmm0
     cdf:   c9                      leave
     ce0:   c3                      ret
</code></pre>
And all you need to know is that the call instruction is calling sqrt@plt instead of sqrtf@plt.<p>And if you don&#x27;t believe the compiler either, I&#x27;ll quote this sentence from the C++11 standard, section 2.14.4 [lex.fcon]:<p>&gt; The type of a floating literal is double unless explicitly specified by a suffix. The suffixes f and F specify float, the suffixes l and L specify long double.<p>So if you do `sqrt(x*x + 1.)` the right operand of plus is double, therefore the argument to sqrt is double, and you get the double version.