template <typename T> T q(T x) {
      return fabs(x - sqrt(x*x + 1.)) - 1./(x + sqrt(x*x + 1.));
    }

The author doesn't seem to have a lot of numeric experience. He probably doesn't know that, even when T is float, the compiler uses the double version of sqrt.

Don't believe me? Read the assembler:

    0000000000000c3c <float q<float>(float)>:
     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 <std::piecewise_construct+0x8>
     c68:   00
     c69:   f2 0f 58 c1             addsd  xmm0,xmm1
     c6d:   e8 5e fc ff ff          call   8d0 <sqrt@plt>
     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 <std::piecewise_construct+0x10>
     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 <std::piecewise_construct+0x8>
     caf:   00
     cb0:   f2 0f 58 c1             addsd  xmm0,xmm1
     cb4:   e8 17 fc ff ff          call   8d0 <sqrt@plt>
     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 <std::piecewise_construct+0x8>
     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

And all you need to know is that the call instruction is calling sqrt@plt instead of sqrtf@plt.

And if you don't believe the compiler either, I'll quote this sentence from the C++11 standard, section 2.14.4 [lex.fcon]:

> 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.

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.