How Math.Pow (and so on) actually works

Question

So I was googling for a long time and i found almost nothing. I found some info about possible implementation of Math.Pow from this url, but they are inaccurate, for example this code

public static double PowerA(double a, double b)
{
    int tmp = (int)(BitConverter.DoubleToInt64Bits(a) >> 32);
    int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);
    return BitConverter.Int64BitsToDouble(((long)tmp2) << 32);
}
static void Main(string[] args)
{
    double x = 12.53, y = 16.45;
    Console.WriteLine(Math.Pow(x, y));
    Console.WriteLine(PowerA(x, y));
}

provides output:

1,15158266266297E+18
8,9966384455562E+17

So inaccurate...

I was thinking that it works like a sum of series but I don't know for certain.

Vanadous

The obligatory What Every Computer Scientist Should Know About Floating-Point Arithmetic. I'm sure there are bits you didn't know about floating point in that paper. – Ergosterol

You can download code that should be equivalent to the actual code, if you want to. SSCLI. It's not been updated for a while, but Math.Pow is fairly old. – Powwow

See How is Math.Pow() implemented in .Net Framework? – Truong

I would suggest that PowerA is the one that is inaccurate. PowerA(12.0, 5.0) returns the incorrect 227008.5 whereas Math.Pow(12.0, 5.0) returns the correct 248832, which is small enough to test. Plus a whole number double raised to a whole number double should return another whole number double, and not something ending with .5! – Luba

@RickDavin I read that as what the OP was saying. The method he has doesn't correctly mirror Math.Pow, it is incorrect, and he's looking for an implementation that will correctly mirror Math.Pow. – Cadel

to: Tim S - i've readed it before asking, but there's only explanation that it's native C++ code so this is why i published this question. to: ppeterka 66 - i've readed an article based on article from your URL. So the greater part of it i know, but i'l defenitly read the original. But how can it help me? Link to disassembler manual should be more practice :) – Vanadous

The linked page does state "In my tests it usually within an error margin of 5% to 12%, in extreme cases sometimes up to 25%." @RickDavin 's results are well within defined and documented error margins. – Depurate

those are approximations not computations of pow which means that it is not accurate on whole RxR space just in some small portion of it. To compute pow you need to transform a from a^b to known base usually 2 so you compute 2^c, for more info see my answer – Duad

Answer 1

pow is usually evaluated by this formula:

x^y = exp2(y*log2(x))

Functions exp2(x),log2(x) are directly implemented in FPU. If you want to implement bignums then they can also be evaluated by basic operators with use of precomputed table of sqrt-powers like:

2^1/2, 2^1/4, 2^1/8, 2^1/16, 2^1/32 ...

to speed up the process

In case you need to handle also rooting for negative bases see this:

• real domain pow based on complex domain math

Duad

functions exp2(x),log2(x) are directly implemented in FPU - it's the most important info that i didn't know. Tnx. – Vanadous

functions exp2(x),log2(x) are directly implemented in FPU - are you sure you don't mean FP library, since FPU is a hardware entity that varies enormously with CPU architecture (not explicit in question) and may NOT NECESSARILY implement high-order functions like this. – Hadwyn

@Eight-BitGuru I work on x87 FPU's and there these functions are present natively , so yes I mean implemented in FPU not in FP lib. – Duad

@Eight-BitGuru see instructions like f2xm1,fyl2x,fyl2xp1 ... by "work on" I mean "I work with" of coarse ... the wrong translation just occur to me :) – Duad

@Duad I think you missed the core of my point - not all FPU architectures implement high-order functions of this type, so asserting that 'functions exp2(x),log2(x) are directly implemented in FPU' is making an assumption about the nature of the hardware which is not explicitly specified in the question. – Hadwyn

@Eight-BitGuru OP tags and text implies PC platform which contains x87 and the question is about pow not exp2,log2 so that is why I do not see your point. if exp2 and log2 are not presents then they can be computed by elementary functions. for example log can be done like this: first exponent is computed (on integers) then FP is converted to fixed point and binary search is applied on mantisa. this is for example sqrt implementation https://mcmap.net/q/833431/-any-faster-rms-value-calculation-in-c in this way. – Duad

@Eight-BitGuru exp2 can be computed directly by table of 2^(1/(2^i)) powers where i is integer>1 value by exploiting 2^(a0+a1+a2+...)=(2^a0)*(2^a1)*... btw mine implementation of log2 use this same table for speed up,... and there are many other ways for sure – Duad

@Spektre: I suppose log2x is base 2 logarithm but what is exp2x that is implemented in FPU? – Excerpt

@Excerpt exp2 is pow(2,x) its also implemented in FPU directly ... you know log,exp are inverse to each other ... log2(exp2(x)) = x – Duad