[OT] Re: Let's kill 80bit real at CTFE

[Manu via Digitalmars-d]

On 12 September 2016 at 00:31, Marco Leise via Digitalmars-d
<digitalmars-d at puremagic.com> wrote:
> Am Sun, 11 Sep 2016 15:00:12 +1000
> schrieb Manu via Digitalmars-d <digitalmars-d at puremagic.com>:
>
>> On 9 September 2016 at 21:50, Stefan Koch via Digitalmars-d
>> <digitalmars-d at puremagic.com> wrote:
>> > Hi,
>> >
>> > In short 80bit real are a real pain to support cross-platform.
>> > emulating them in software is prohibitively slow, and more importantly hard
>> > to get right.
>> > 64bit floating-point numbers are supported on more architectures and are
>> > much better supported.
>> > They are also trivial to use at ctfe.
>> > I vote for killing the 80bit handling at constant folding.
>> >
>> > Destroy!
>>
>> I just want CTFE '^^'. Please, can we have that?
>> It's impossible to CTFE any non-linear function. It's ridiculous, I
>> constantly want to generate a curve at compile time!
>
> I have experimented with a few iterative algorithms from
> around the web that are now in my module for "random stuff not
> in Phobos":
>
> /*******************************************************************************
>  *
>  * Computes the arcus tangens at compile time.
>  *
>  **************************************/
> enum ctfeAtan(real x)
> in
> {
>         assert(x == x && abs(x) != real.infinity);
> }
> body
> {
>         if (abs(x) == 0)
>                 return x;
>
>         // Reduce x to <0.5 for effective convergence of the Taylor series.
>         x /= 1 + sqrt(1 + x * x);
>         x /= 1 + sqrt(1 + x * x);
>
>         // Sum up Taylor series to compute atan().
>         immutable xSqr = -x * x;
>         real mul = x;
>         real div = 1;
>         real x_old;
>         do
>         {
>                 x_old = x;
>                 mul *= xSqr;
>                 div += 2;
>                 x += mul / div;
>         }
>         while (x !is x_old);
>
>         // Compensate for the initial reduction by multiplying by 4.
>         return 4 * x;
> }
>
>
> /*******************************************************************************
>  *
>  * Computes the arcs sinus at compile time.
>  *
>  **************************************/
> enum ctfeAsin(real x)
> in
> {
>         assert(x.isWithin(-1,+1));
> }
> body
> {
>         if (abs(x) == 0)
>                 return x;
>
>         immutable div = 1 - x * x;
>         return x / abs(x) * (div == 0 ? PI / 2 : ctfeAtan(sqrt(x * x / div)));
> }
>
>
> /*******************************************************************************
>  *
>  * Computes `x` to the power of `y` at compile-time.
>  *
>  * Params:
>  *   x = The base value.
>  *   y = The power.
>  *
>  * Source:
>  *   http://stackoverflow.com/a/7710097/4038614
>  *
>  **************************************/
> @safe @nogc pure nothrow
> real ctfePow(real x, real y)
> {
>         if (y >= 1)
>         {
>                 real temp = ctfePow( x, y / 2 );
>                 return temp * temp;
>         }
>
>         real low = 0, high = 1;
>         real sqr = sqrt( x );
>         real acc = sqr;
>         real mid = high / 2;
>
>         while (mid != y)
>         {
>                 sqr = sqrt( sqr );
>
>                 if (mid <= y)
>                 {
>                         low = mid;
>                         acc *= sqr;
>                 }
>                 else
>                 {
>                         high = mid;
>                         acc *= 1 / sqr;
>                 }
>
>                 mid = (low + high) / 2;
>         }
>
>         return acc;
> }
>
>
> /*******************************************************************************
>  *
>  * Computes the natural logarithm of `x`at compile time.
>  *
>  **************************************/
> @safe @nogc pure nothrow
> FloatReg ctfeLog(FloatReg x)
> {
>         if (x != x || x <= 0)
>                 return -FloatReg.nan;
>         else if (x == 0)
>                 return -FloatReg.infinity;
>         else if (x == FloatReg.infinity)
>                 return +FloatReg.infinity;
>
>         uint m = 0;
>         while (x <= ulong.max)
>         {
>                 x *= 2;
>                 m++;
>         }
>
>         @safe @nogc pure nothrow
>         static FloatReg agm(FloatReg x, FloatReg y)
>         in
>         {
>                 assert(x >= y);
>         }
>         body
>         {
>                 real a, g;
>                 do
>                 {
>                         a = x; g = y;
>                         x = 0.5 * (a+g);
>                         y = sqrt(a*g);
>                 }
>                 while(x != a || y != g);
>                 return x;
>         }
>
>         return PI_2 / agm(1, 4/x) - m * LN2;
> }
>
>
> Especially the `log` function seemed like a good compromise
> between execution speed and accuracy. FloatReg is "the native
> float register type", but the way CTFE works it should be
> `real` anyways. The code is mostly not by me, but from
> StackOVerflow comments and Wikipedia. Most of it is common
> knowledge to every mathematician.
>
> --
> Marco

That's cool, but surely unnecessary; the compiler should just hook
these and do it directly... They're intrinsics in every
compiler/language I've ever used! Just not DMD.
If your results are compatible, why not PR this implementation for
when ctfe in std.math?
Incidentally, this is how we used to do these operations in early
shader code, except cutting some corners for speed+imprecision ;)