uint32_t n is const, no?

similar to the above

I don't think that would add anything, really.

Sure, I can can do this. It won't change anything in practice though, ints are passed by value so the caller really don't care either way.

Sure but this is a bit OCD at this point. We are not using this at compile time anyways.

Sure, but then we need to put the exact size in there so changing the test is more burdensome. I'm not convinced this adds much value.

reative -> relative

aproximation -> approximation

aproximate -> approximate

sevral -> several

sevral -> several

require -> required

A better name would be ONE_32, as it is representing the fixed point number 1.

Using work instead of target introduces more rounding errors than using target directly and makes it more complicated porting this code, as you need division on 256 bit numbers.

It would be much easier and more precise to work directly on the bits representation of the difficulty.

As noted by Marc Lundeberg, why not use a simple linear approximation. There is no reason why an exponential should work better here. The absolute difference between a linear function and the exponential should also be just a few seconds.

The implicit cast to uint32_t here that silently adds 2^32 for negative x32 must be documented. The comment should also explain how the effect of this is taken care of by the right shift in line 138 when x32 is negative.

For my taste these bit-fiddling tricks are dangerous and will lead to subtle bugs when trying to port this code to other languages. Note also that the earliest time such a bug will trigger for the main chain is in several years when the blocks will get back to normal 10 minutes intervals.

These constants are almost the same as ((k0 + POW2_32) * LN2_32) >> 32. Does the extra precision really merit doubling the number of magic constants?

Why the >> 1 here and then >>31 in the next line. Are you afraid of overflows? The prescaling should ensure that x is less than 2^27, so there is no overflow in u1+u2.

It deserves more comments, e.g., that the halfing of u1 is counteracted by the >>31 in the next line, while for u2 there is >> 2 because the Taylor series contains 0.5x^2.

Why so complicated? Something like

// scale x so that exp(xscaled) = 2^x, i.e. xscaled= x*ln(2).
int64_t xscaled = (x*LN2_32)>>32;
// approximate exp(xscaled)-1 by the Taylor series x + x^2/2 (the right shift by 33 divides x^2 by 2).
int64_t u = xscaled + (xscaled * xscaled)>>33;
// multiply result by 1+k0.  Note that we need to return (1+k0)*exp(xscaled) - 1 = (1+k0)*(exp(xscaled)-1) +  k0
return k0 + (u * (k0 + POW2_32))>>32

should be similarly accurate, doesn't need k1, and is probably easier to understand.

Why the >> 1 here and then >>31 in the next line. Are you afraid of overflows? The prescaling should ensure that x is less than 2^27, so there is no overflow in u1+u2.

I'm afraid you are incorrect on this one, you can definitively get an oveflow for u1+u2 . This is something I double checked for.

I will play with some of the proposed simplification and see how it goes.

For my taste these bit-fiddling tricks are dangerous and will lead to subtle bugs when trying to port this code to other languages.

This is interesting. My experience is somewhat the reverse. Not that more bit fiddling is always good - obviously not - but that trying to avoid them and pretend that computer math somehow resemble "real" math is a common source of problems. I will definitively add a comment to explain the right shift business, however. Good call.

Copy over the comment from the header file? Maybe also expand a bit.

This function returns an approximation of 2^n - 1 in 32-bit fixed point notation. Thus the argument is is interpreted as the number (n/2^32) in the interval [0, 1[. The result is also in the interval [0, 1[ and scaled by 2^32.

In other words, this function returns an approximation of (2^(n / 2^32) - 1) * 2^32.

nitpicking my own code :)
It returns exactly the same value but is a tiny bit simpler.

Beware that this diff doesn't delete all old lines. Is there a way to do a proper multi-line diff in phabricator?

While your proposal was very good, @markblundeberg ended up sending me something that is both simpler and more precise, that I'm working on integrating now. The progress made here is fantastic, thanks to you and Mark.

buket => bucket

I have my own code to compute Remez-approximated polynomials. It produces a slightly more precise polynom:
2977042781 * x + 1031830288 *x**2 + 237606570 * x**3 + 44970684*x**4.
The difference is small (3.1e-9 vs. 3.8e-9), though. Main purpose of this exercise was to verify the code.

Here is the gnuplot command for plotting relative error (my code optimizes for absolute error, but the difference is small).

plot [0:.25] (1+2977042554./2**32*x + 1031834945./2**32*x**2 + 237575834./2**32*x**3 + 45037112./2**32*x**4)/exp(x*log(2))-1,(1+2977042781./2**32*x + 1031830288./2**32*x**2 + 237606570./2**32*x**3 + 44970684./2**32*x**4)/exp(x*log(2))-1

OTOH, the polynomial from Mark is precise for x=.25, which may be important for monotonicity. If I optimize with this constraint I get:
2977042597*x + 1031835569*x**2 + 237564652*x**3 + 45069083*x**4. The error is around 3.5e-9 so still very slightly better.

plot [0:.25] (1+2977042554./2**32*x + 1031834945./2**32*x**2 + 237575834./2**32*x**3 + 45037112./2**32*x**4)/exp(x*log(2))-1,(1+2977042597./2**32*x + 1031835569./2**32*x**2 + 237564652./2**32*x**3 + 45069083./2**32*x**4)/exp(x*log(2))-1

I guess my suggestion is to use these coefficients: 2977042597, 1031835569, 237564652, 45069083, unless Mark has a reason why his coefficients are better. The difference as I said is very tiny.

If you want to go fancy, you can even get to 6e-10 on the full interval [-.5,.5]
by using an approximation for R(x) = x + 2x/(2^x-1), which is easier to approximate with polynomials:

R(x) * 2^30 = 3098164009 + 124043115 * x^2 + -988834 * x^4

Then compute 2^x - 1 = 2x/(R(x)-x):

plot [-.5:.5] (1 + 2*x*2**30/(3098164009. + 124043115*x**2 + -988834 * x**4 - 2**30 * x)) / (exp(x*log(2)))-1

This requires rescaling x to a signed integer.

One of the problem with that approach is that x^4 does not fit in a natively supported type. So we' end up with a ton of bit manips, or we'd shift out a few bits, in which case we'll get back to something similar to what Mark produced.

Simplicity is also a benefit here. I'm not convinced the extra precision justify the extra complexity.