How is Math.Pow() implemented in .NET Framework?
Asked Answered
E

5

447

I was looking for an efficient approach for calculating ab (say a = 2 and b = 50). To start things up, I decided to take a look at the implementation of Math.Pow() function. But in .NET Reflector, all I found was this:

[MethodImpl(MethodImplOptions.InternalCall), SecuritySafeCritical]
public static extern double Pow(double x, double y);

What are some of the resources wherein I can see as what's going on inside when I call Math.Pow() function?

Exarch answered 15/1, 2012 at 14:31 Comment(4)
Just as an FYI, if you're confused about the whole InternalCall with an extern modifier (as they seem to be conflicting), please see the question (and the resulting answers) that I posted about this very same thing.Playgoer
For a 2^x operation if x is integer the result is a shift operation. So maybe you could construct the result using a mantissa of 2 and an exponent of x.Aachen
@SurajJain your comment is actually a question you need to post separately.Aachen
@SurajJain I agree with you. I am not a moderator so I can't do much here. Maybe the downvote question can be asked at meta.stackoverflow.comAachen
M
888

MethodImplOptions.InternalCall

That means that the method is actually implemented in the CLR, written in C++. The just-in-time compiler consults a table with internally implemented methods and compiles the call to the C++ function directly.

Having a look at the code requires the source code for the CLR. You can get that from the SSCLI20 distribution. It was written around the .NET 2.0 time frame, I've found the low-level implementations, like Math.Pow() to be still largely accurate for later versions of the CLR.

The lookup table is located in clr/src/vm/ecall.cpp. The section that's relevant to Math.Pow() looks like this:

FCFuncStart(gMathFuncs)
    FCIntrinsic("Sin", COMDouble::Sin, CORINFO_INTRINSIC_Sin)
    FCIntrinsic("Cos", COMDouble::Cos, CORINFO_INTRINSIC_Cos)
    FCIntrinsic("Sqrt", COMDouble::Sqrt, CORINFO_INTRINSIC_Sqrt)
    FCIntrinsic("Round", COMDouble::Round, CORINFO_INTRINSIC_Round)
    FCIntrinsicSig("Abs", &gsig_SM_Flt_RetFlt, COMDouble::AbsFlt, CORINFO_INTRINSIC_Abs)
    FCIntrinsicSig("Abs", &gsig_SM_Dbl_RetDbl, COMDouble::AbsDbl, CORINFO_INTRINSIC_Abs)
    FCFuncElement("Exp", COMDouble::Exp)
    FCFuncElement("Pow", COMDouble::Pow)
    // etc..
FCFuncEnd()

Searching for "COMDouble" takes you to clr/src/classlibnative/float/comfloat.cpp. I'll spare you the code, just have a look for yourself. It basically checks for corner cases, then calls the CRT's version of pow().

The only other implementation detail that's interesting is the FCIntrinsic macro in the table. That's a hint that the jitter may implement the function as an intrinsic. In other words, substitute the function call with a floating point machine code instruction. Which is not the case for Pow(), there is no FPU instruction for it. But certainly for the other simple operations. Notable is that this can make floating point math in C# substantially faster than the same code in C++, check this answer for the reason why.

By the way, the source code for the CRT is also available if you have the full version of Visual Studio vc/crt/src directory. You'll hit the wall on pow() though, Microsoft purchased that code from Intel. Doing a better job than the Intel engineers is unlikely. Although my high-school book's identity was twice as fast when I tried it:

public static double FasterPow(double x, double y) {
    return Math.Exp(y * Math.Log(x));
}

But not a true substitute because it accumulates error from 3 floating point operations and doesn't deal with the weirdo domain problems that Pow() has. Like 0^0 and -Infinity raised to any power.

Milomilon answered 15/1, 2012 at 14:52 Comment(16)
Great answer, StackOverflow needs more of this sort of thing, instead of 'Why would you want to know that?' that happens all too often.Salo
@Hans I'm not sure I can follow the last paragraph: If I look into crt/src/math.h I do actually find _Pow_int (~line 493 on my vs2010 install) which at least it seems to me is used for all the pow calls. And that seems like the obvious shift and multiply implementation. Am I missing something here? (I assume an optimized intel implementation would use some sort of SSE sorcery? It seems like we could parallelize it a bit, but not sure if that would be an improvement)Rarefy
@Rarefy - that's a template specialization for pow() calls that pass an integer power. Like pow(x,2) = x * x. Yes, the Intel version uses quite convoluted SSE2 machine code.Milomilon
@Hans Ups I completely missed the part where the power was a double. Mea culpa.Rarefy
@Blue - I don't know, short from making fun of Intel engineers. My high school book does have a problem raising something to the power of a negative integral. Pow(x, -2) is perfectly computable, Pow(x, -2.1) is undefined. Domain problems are a bitch to deal with.Milomilon
@BlueRaja-DannyPflughoeft: A lot of effort is spent trying to ensure that floating-point operations are as close as possible to the correctly-rounded value. pow is notoriously hard to implement accurately, being a transcendental function (see Table-Maker's Dilemma). It's a lot easier with an integral power.Clypeus
@Hans Passant: Why would Pow(x,-2.1) be undefined? Mathematically pow is defined everywhere for all x and y. You do tend to get complex numbers for negative x and non integer y.Viens
@Viens pow(0, 0) is not defined.Degraw
@emboss: right, that's the only case. Though some mathematicians define it as 1 or 0, pow does have an essential discontinuity at (0,0) so there is no unique way to patch it (so the best thing is probably indeed not to patch it at all).Viens
@Viens Complex numbers just got added to the Framework for 4.0, msdn.microsoft.com/en-us/library/…. It now defines Pow(Complex, Double) and Pow(Complex, Complex).Thorma
@Thorma Did they? That would have been so handy five years ago when I was writing geospatial functions for DIY tiling in a desktop app. Thankfully Google obsoleted it before I went mad. Madder.Spoken
@HansPassant - I think you meant Pow(x, -2.1) is undefined (or complex) for negative values of x, right? The sign on the exponent doesn't make it any more or less computable unless x is 0.Without
Is this still true as of roslyn? I see coreclr on Github still shows this as an internal call: github.com/dotnet/coreclr/blob/…Mcnelly
Look at the question, note that it is identical.Milomilon
I think this is the actual C++ CLR implementation (at least for coreclr) github.com/dotnet/coreclr/blob/…Mcnelly
github.com/dotnet/coreclr/blob/…Milomilon
C
115

Hans Passant's answer is great, but I would like to add that if b is an integer, then a^b can be computed very efficiently with binary decomposition. Here's a modified version from Henry Warren's Hacker's Delight:

public static int iexp(int a, uint b) {
    int y = 1;

    while(true) {
        if ((b & 1) != 0) y = a*y;
        b = b >> 1;
        if (b == 0) return y;
        a *= a;
    }    
}

He notes that this operation is optimal (does the minimum number of arithmetic or logical operations) for all b < 15. Also there is no known solution to the general problem of finding an optimal sequence of factors to compute a^b for any b other than an extensive search. It's an NP-Hard problem. So basically that means that the binary decomposition is as good as it gets.

Cimah answered 10/7, 2012 at 18:15 Comment(3)
This algorithm(square and multiply) also applies if a is a floating point number.Mcmillin
In practice it is possible to do quite a bit better than native square-and-multiply. For example preparing lookup tables for small exponents so you can square several times and only then multiply, or building optimized square-addition chains for fixed exponents. This kind of problem is integral to important cryptographic algorithms, so there as been quite a bit of work on optimizing it. NP hardness is only about worst-case asymptotics, we often can produce optimal or near-optimimal solutions for instances of the problem arising in practice.Mcmillin
The text doesn't mention a being an integer, but the code does. As a consequence of that, I wonder about the accuracy of the result of the text's "very efficient" computation.Conney
S
71

If freely available C version of pow is any indication, it does not look like anything you would expect. It would not be of much help to you to find the .NET version, because the problem that you are solving (i.e. the one with integers) is orders of magnitudes simpler, and can be solved in a few lines of C# code with the exponentiation by squaring algorithm.

Silurian answered 15/1, 2012 at 14:43 Comment(2)
Thanks for your answer. The first link surprised me as I wasn't expecting such massive technical implementation of Pow() function. Although Hans Passant answer confirms that its the same in .Net world too. I think I can solve the problem at hand by making use of some of the technique listed in the squaring algorithm link. Thanks Again.Exarch
I don't believe that this code is efficient. 30 local variables just should bump all registers. I only suppose that it's ARM version, but on x86 30 local variables in the method is awesome.Idaline
A
1

Going through the answers, learned a lot about behind-the-scene calculations: I've tried some workarounds on a coding platform which has an extensive test coverage cases, and found a very effective way doing it(Solution 3):

public double MyPow(double x, int n) {
    double res = 1;
    /* Solution 1: iterative : TLE(Time Limit Exceeded)
    double res = 1;
    var len = n > 0 ? n : -n;
    for(var i = 0; i < len; ++i)
        res *= x;   
    
    return n > 0 ? res : 1 / res; 
    */
    
    /* Solution 2: recursive => stackoverflow exception
    if(x == 0) return n > 0 ? 0 : 1 / x;
    if(n == 1) return x;
    
    return n > 0 ? x * MyPow(x, n - 1) : (1/x) * MyPow(1/x, -n); 
    */
    
    //Solution 3:
    if (n == 0) return 1;
    
    var half = MyPow(x, n / 2);
    if (n % 2 == 0) 
        return half * half;
    else if (n > 0) 
        return half * half * x;
    else 
        return half * half / x;
    
    /* Solution 4: bitwise=> TLE(Time Limit Exceeded)
    var b = n > 0 ? n : -n;        
    while(true) {
        if ((b & 1) != 0) 
            res *= x;
        
        b = b >> 1;
        
        if (b == 0) break;
        
        x *= x;
    }   
    return n > 0 ? res : 1 / res; 
    */
}
Alanealanine answered 16/7, 2020 at 14:49 Comment(0)
C
0

Answer that is accepted on Leetcode:

public class Solution {
    public double MyPow(double x, int n) {
        if(n==0) return 1;
        
        long abs = Math.Abs((long)n);
        
        var result = pow(x, abs);
            
        return n > 0 ? result : 1/result;
    }
    
    double pow(double x, long n){
        if(n == 1) return x;
        
        var result = pow(x, n/2);
        result = result * result * (n%2 == 1? x : 1);
        return result;
    }
}
Conlee answered 24/4, 2022 at 15:33 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.