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

[Marco Leise via Digitalmars-d]

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