Thread overview
Let's kill 80bit real at CTFE
Sep 09, 2016
Stefan Koch
Sep 09, 2016
rikki cattermole
Sep 11, 2016
tsbockman
Sep 11, 2016
Manu
[OT] Re: Let's kill 80bit real at CTFE
Sep 11, 2016
Marco Leise
Sep 11, 2016
Manu
Sep 11, 2016
Stefan Koch
September 09, 2016
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!
September 09, 2016
On 09/09/2016 11:50 PM, Stefan Koch 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 only need to say one thing which this will utterly break. std.math.
Done!
September 11, 2016
On Friday, 9 September 2016 at 11:50:08 UTC, Stefan Koch wrote:
> Hi,
>
> In short 80bit real are a real pain to support cross-platform.
> emulating them in software is prohibitively slow,

Supposedly, well-optimized 128-bit software floating-point is actually a bit faster than hardware 80-bit. (Intel keeps 80-bit around primarily for backwards compatibility, and hasn't put much effort or die space into improving its hardware speed in many years.)

> 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!

This has been debated to death in the past. Reducing the maximum precision is not an acceptable option, but 80-bit could potentially be replaced by a software 128-bit implementation if that helps.

September 11, 2016
On 9 September 2016 at 21:50, Stefan Koch via Digitalmars-d <digitalmars-d@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!
September 11, 2016
Am Sun, 11 Sep 2016 15:00:12 +1000
schrieb Manu via Digitalmars-d <digitalmars-d@puremagic.com>:

> On 9 September 2016 at 21:50, Stefan Koch via Digitalmars-d <digitalmars-d@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

September 12, 2016
On 12 September 2016 at 00:31, Marco Leise via Digitalmars-d <digitalmars-d@puremagic.com> wrote:
> Am Sun, 11 Sep 2016 15:00:12 +1000
> schrieb Manu via Digitalmars-d <digitalmars-d@puremagic.com>:
>
>> On 9 September 2016 at 21:50, Stefan Koch via Digitalmars-d <digitalmars-d@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 ;)
September 11, 2016
On Sunday, 11 September 2016 at 14:52:18 UTC, Manu wrote:
>
> 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 ;)

Hey I just made another PR, to fix PowExp!
It might as well depend on phobos :)