Square root metafunction?

Question

Is it possible to compute the square root of an integer with a metafunction with the following signature :

template<unsigned int N> inline double sqrt();

(or maybe using the constexpr keyword, I don't know what is the best). With that, sqrt<2>() would be replaced by 1.414... at compile-time.

What would be the best implementation for a such function ?

Sudarium

I googled "template metaprogramming sqrt" and found this informit.com/articles/article.aspx?p=30667&seqNum=3 – Storer

I've already seen that but it is only for the integral part of sqrt of an integer. I would like to have the floating-point result at compile-time. – Sudarium

Since sqrt is a standard function I'd use sqrtt instead. – Tetragonal

Can you do any meaningful manipulation on floating point numbers at compile time? – Ravo

my metaprogramming library includes a square root function implemented with the newton root finding method. (Check the README examples). – Chopping

Answer 1

This may not be what you are looking for, but I wanted to make sure you realized that typically with optimization the compiler will calculate the result at compile time anyway. For example, if you have this code:

void g()
{
  f(sqrt(42));
}

With g++ 4.6.3 with optimization -O2, the resulting assembly code is:

   9 0000 83EC1C                subl    $28, %esp
  11 0003 DD050000              fldl    .LC0
  12 0009 DD1C24                fstpl   (%esp)
  13 000c E8FCFFFF              call    _Z1fd
  14 0011 83C41C                addl    $28, %esp
  16 0014 C3                    ret
  73                    .LC0:
  74 0000 6412264A              .long   1244009060
  75 0004 47EC1940              .long   1075440711

The sqrt function is never actually called, and the value is just stored as part of the program.

Therefore to create a function that technically meets your requirements, you simply would need:

template<unsigned int N> inline double meta_sqrt() { return sqrt(N); }

Uzia

It may be true that most compilers calculate this at compile time, but that is not a guarantee. The sqrt function does not return a constexpr value; the output to this function cannot be used as an input to other metafunctions. This is not truly a metafunction and does not answer the question. – Fascine

Answer 2

Eigen contains this meta_sqrt which uses binary search:

template<int Y,
         int InfX = 0,
         int SupX = ((Y==1) ? 1 : Y/2),
         bool Done = ((SupX-InfX)<=1 ? true : ((SupX*SupX <= Y) && ((SupX+1)*(SupX+1) > Y))) >
                                // use ?: instead of || just to shut up a stupid gcc 4.3 warning
class meta_sqrt
{
    enum {
  MidX = (InfX+SupX)/2,
  TakeInf = MidX*MidX > Y ? 1 : 0,
  NewInf = int(TakeInf) ? InfX : int(MidX),
  NewSup = int(TakeInf) ? int(MidX) : SupX
};
  public:
    enum { ret = meta_sqrt<Y,NewInf,NewSup>::ret };
};

template<int Y, int InfX, int SupX>
class meta_sqrt<Y, InfX, SupX, true>
{
    public:  enum { ret = (SupX*SupX <= Y) ? SupX : InfX };
};

Rhynd

Answer 3

The issue I see is that metaP effectively abuses enums into variables. The problem is that enums are treated internally as integers which precludes trying to get a floating point value out of them. However, you may be able to create your own floating point format that creates two results, an integer portion and an exponent. You'll still have to process this into a float, as Out = Sqrt<42>::mantissa * pow(10,Sqrt<42>::exponent);. Actually determining the values is left as an exercise for the reader, but you will probably have to scale the input upwards (by an even power of 10), computing the root, and storing the -power/2 you used earlier.

To compute sqrt<42>, you would first set the exponent enum to a suitable power, such as '-4' (the lower, the more decimals, but watch for overflow). You then multiply the input by '10^(-2*exponent)'. In this case you get 42*10^8 = 4200000000. You then take the root of this value getting '64807' as the final value. At runtime, you calculate the "val * 10 ^ exponent" = "64807 * 10 ^ -4" = 64807 * 0.0001 = 6.4807m and store it to a float.

The extra conversion work kinda defeats the purpose, but you can reduce it a bit by storing the exponent as 10^k (ie 10^4) then doing out=sqrt<x>::mantissa/sqrt<x>::exponent.

edit I just noticed that with the mantissa/exponent method, the choice of exponent is arbitrary as long as it's bigger than the integer portion of the final root. It can even be a constant, which eases the design of your meta functions. In the case of 42 for example you could select the 'exponent' to always be 6000. You then multiply the input by 6000^2, take the integer root of the product, then at runtime, divide the result by 6000 to get the root. Instead of treating the output as a*10^b, it uses the relation sqr(x*b^2)=sqr(x)*b. The math checks out:

• 42*6000*6000=1512000000
• sqr(1512000000)=38884
• 38884/6000=6.4806 (squared is 41.999)

Impletion

You know the computer doesn't use powers of 10 for storing floating-point values, right? Also, you've got the concepts of exponent and mantissa swapped. – Pacify

Fixed :) I know about the power of ten thing, but since I'm working with decimals anyway, it doesn't hurt. Also added that the choice of exponent is arbitrary. – Impletion

But if you use powers of two, you can use ldexp to combine the exponent and mantissa without requiring any multiplication or division. – Pacify