Fast trunc() algorithm

From	Bonita Montero <Bonita.Montero@gmail.com>
Newsgroups	comp.lang.c++, comp.lang.c, de.comp.lang.c
Subject	Fast trunc() algorithm
Date	2022-05-29 10:14 +0200
Organization	A noiseless patient Spider
Message-ID	<t6va0v$eno$1@dont-email.me> (permalink)

Cross-posted to 3 groups.

Show all headers | View raw

The trunc() implementation of VC++'s runtime was too slow for me,
so I've written my own:

#include <cstdint>
#include <limits>
#include <cfenv>
#include <cstring>

using namespace std;

double ftrunc( double d )
{
	static_assert(sizeof(double) == 8 && numeric_limits<double>::is_iec559, 
"double must be IEEE-754");
	// assume size_t is our CPU's native register-width
	static_assert(sizeof(size_t) == 8 || sizeof(size_t) == 4, 
"register-width must be 32 or 64 bit");
	if constexpr( sizeof(size_t) == 8 )
	// we have 64 bit registers
	{
		unsigned const
			MANTISSA_BITS = 52,
			EXP_BIAS = 0x3FF,
			INF_NAN_BASE = 0x7FF;
		uint64_t const
			EXP_MASK = (uint64_t)0x7FF << MANTISSA_BITS,
			SIGN_MASK = (uint64_t)0x800 << MANTISSA_BITS ,
			MANTISSA_MASK = 0x000FFFFFFFFFFFFFu,
			NAN_MASK = 0x0008000000000000u,
			MIN_INTEGRAL_DIGITS_EXP = (uint64_t)EXP_BIAS << MANTISSA_BITS,
			MIN_INTEGRAL_ONLY_EXP = (uint64_t)(EXP_BIAS + MANTISSA_BITS) << 
MANTISSA_BITS,
			INF_NAN_EXP = (uint64_t)INF_NAN_BASE << MANTISSA_BITS;
		int64_t const  MANTISSA_SHIFT_MASK = 0xFFF0000000000000u;
		uint64_t dx;
		memcpy( &dx, &d, 8 );
		auto retBack = [&]() -> double
		{
			memcpy( &d, &dx, 8 );
			return d;
		};
		uint64_t exp = dx & EXP_MASK;
		if( exp >= MIN_INTEGRAL_DIGITS_EXP )
			// value has integral digits
			if( exp < MIN_INTEGRAL_ONLY_EXP )
			{
				// there are fraction-digits to mask out, mask them
				unsigned shift = (unsigned)(exp >> MANTISSA_BITS) - EXP_BIAS;
				dx &= MANTISSA_SHIFT_MASK >> shift;
				return retBack();
			}
			else
				if( exp < INF_NAN_EXP )
					// value is integral
					return d;
				else
				{
					uint64_t mantissa = dx & MANTISSA_MASK;
					// infinite, NaN: return value
					if( !mantissa || mantissa & NAN_MASK )
						return retBack();
					// SNaN: raise exception on SNaN if necessary
					feraiseexcept( FE_INVALID );
					return retBack();
				}
		else
		{
			// below +/-1.0
			// return +/-0.0
			dx &= SIGN_MASK;
			return retBack();
		}
	}
	else
	// we have 32 bit registers
	{
		unsigned const
			MANTISSA_BITS = 52,
			HI_MANTISSA_BITS = 20,
			EXP_BIAS = 0x3FF,
			INF_NAN_BASE = 0x7FF;
		uint32_t const
			EXP_MASK = (uint32_t)0x7FFu << HI_MANTISSA_BITS,
			SIGN_MASK = (uint32_t)0x800u << HI_MANTISSA_BITS,
			HI_MANTISSA_MASK = 0x000FFFFFu,
			NAN_MASK = 0x00080000u,
			MIN_INTEGRAL_DIGITS_EXP = (uint32_t) EXP_BIAS << HI_MANTISSA_BITS,
			MAX_INTEGRAL32_EXP = (uint32_t)(EXP_BIAS + HI_MANTISSA_BITS) << 
HI_MANTISSA_BITS,
			MIN_INTEGRAL_ONLY_EXP = (uint32_t)(EXP_BIAS + MANTISSA_BITS) << 
HI_MANTISSA_BITS,
			INF_NAN_EXP = (uint32_t)INF_NAN_BASE << HI_MANTISSA_BITS,
			NEG_LO_MANTISSA_SHIFT_MASK = 0xFFFFFFFFu;
		int32_t const HI_MANTISSA_SHIFT_MASK = 0xFFF00000;
		uint64_t dx;
		memcpy( &dx, &d, 8 );
		uint32_t
			lo = (uint32_t)dx,
			hi = (uint32_t)(dx >> 32);
		auto retBack = [&]() -> double
		{
			dx = lo | (uint64_t)hi << 32;
			memcpy( &d, &dx, 8 );
			return d;
		};
		uint32_t exp = hi & EXP_MASK;
		if( exp >= MIN_INTEGRAL_DIGITS_EXP )
			// value has integral digits
			if( exp < MIN_INTEGRAL_ONLY_EXP )
				// there are fraction-digits to mask out
				if( exp <= MAX_INTEGRAL32_EXP )
				{
					// the fraction digits are in the upper dword, mask them and zero 
the lower dword
					unsigned shift = (unsigned)(exp >> HI_MANTISSA_BITS) - EXP_BIAS;
					hi &= HI_MANTISSA_SHIFT_MASK >> shift;
					lo = 0;
					return retBack();
				}
				else
				{
					// the fraction digits are in the lower dword, mask them
					unsigned shift = (unsigned)(exp >> HI_MANTISSA_BITS) - EXP_BIAS - 
HI_MANTISSA_BITS;
					lo &= ~(NEG_LO_MANTISSA_SHIFT_MASK >> shift);
					return retBack();
				}
			else
				if( exp < INF_NAN_EXP )
					// value is integral
					return retBack();
				else
				{
					uint32_t hiMantissa = hi & HI_MANTISSA_MASK;
					// infinite, NaN: return value
					if( !(hiMantissa | lo) || hiMantissa & NAN_MASK )
						return retBack();
					// SNaN: raise exception on SNaN if necessary
					feraiseexcept( FE_INVALID );
					return retBack();
				}
		else
		{
			// below +/-1.0
			// return +/-0.0
			hi &= SIGN_MASK;
			lo = 0;
			return retBack();
		}
	}
}

float ftrunc( float f )
{
	static_assert(sizeof(float) == 4, "sizeof(float) not equal to 
sizeof(uint32_t)");
	static_assert(numeric_limits<float>::is_iec559, "float must be IEEE-754");
	unsigned const
		MANTISSA_BITS = 23,
		EXP_BIAS = 0x7F,
		INF_NAN_BASE = 0xFF;
	uint32_t const
		EXP_MASK = (uint32_t)0xFF << MANTISSA_BITS,
		SIGN_MASK = (uint32_t)0x100 << MANTISSA_BITS,
		MANTISSA_MASK = 0x007FFFFFu,
		NAN_MASK = 0x00400000u,
		MIN_INTEGRAL_DIGITS_EXP = (uint32_t) EXP_BIAS << MANTISSA_BITS,
		MIN_INTEGRAL_ONLY_EXP = (uint32_t)(EXP_BIAS + MANTISSA_BITS) << 
MANTISSA_BITS,
		INF_NAN_EXP = (uint32_t)INF_NAN_BASE << MANTISSA_BITS;
	int32_t const MANTISSA_SHIFT_MASK = 0xFF800000u;
	uint32_t fx;
	memcpy( &fx, &f, 4 );
	auto retBack = [&]() -> float
	{
		memcpy( &f, &fx, 4 );
		return f;
	};
	uint32_t exp = fx & EXP_MASK;
	if( exp >= MIN_INTEGRAL_DIGITS_EXP )
		// value has integral digits
		if( exp < MIN_INTEGRAL_ONLY_EXP )
		{
			// there are fraction-digits to mask out, mask them
			unsigned shift = (unsigned)(exp >> MANTISSA_BITS) - EXP_BIAS;
			fx &= MANTISSA_SHIFT_MASK >> shift;
			return retBack();
		}
		else
			if( exp < INF_NAN_EXP )
				// value is integral
				return retBack();
			else
			{
				uint32_t mantissa = fx & MANTISSA_MASK;
				// infinite, NaN: return value
				if( !mantissa || mantissa & NAN_MASK )
					return retBack();
				// SNaN: raise exception on SNaN if necessary
				feraiseexcept( FE_INVALID );
				return retBack();
			}
	else
	{
		// below +/-1.0
		// return +/-0.0
		fx &= SIGN_MASK;
		return retBack();
	}
}

memcpy() is the only valid way to have writeable aliasing in C++.
memcpy() is normally an intrinsic function and in the above code
the memcpy()s between integral and floating point types are just
moves between general purpose and FPU-registers.

Back to comp.lang.c++ | Previous | Next — Next in thread | Find similar | Unroll thread

Thread

Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-29 10:14 +0200
  Re: Fast trunc() algorithm Muttley@dastardlyhq.com - 2022-05-29 08:21 +0000
    Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-29 10:29 +0200
      Re: Fast trunc() algorithm Freethinker <freethinker@mymail.com> - 2022-05-29 15:20 +0200
        Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-29 15:27 +0200
          Re: Fast trunc() algorithm Freethinker <freethinker@mymail.com> - 2022-05-29 15:37 +0200
            Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-29 15:46 +0200
    Re: Fast trunc() algorithm John McCue <jmccue@magnetar.hsd1.ma.comcast.net> - 2022-05-29 21:32 +0000
  Re: Fast trunc() algorithm Christian Gollwitzer <auriocus@gmx.de> - 2022-05-29 23:36 +0200
    Re: Fast trunc() algorithm Juha Nieminen <nospam@thanks.invalid> - 2022-05-30 06:22 +0000
      Re: Fast trunc() algorithm David Brown <david.brown@hesbynett.no> - 2022-05-30 08:59 +0200
      Re: Fast trunc() algorithm muttley@dastardlyhq.com - 2022-05-30 07:52 +0000
      Re: Fast trunc() algorithm Paavo Helde <eesnimi@osa.pri.ee> - 2022-05-30 11:58 +0300
    Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-30 11:32 +0200
  Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-30 11:34 +0200
  Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-05-31 05:41 +0200
    Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-06-01 18:07 +0200
      Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-06-01 18:09 +0200
      Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-06-02 04:43 +0200
        Re: Fast trunc() algorithm jak <nospam@please.ty> - 2022-06-02 07:44 +0200
          Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-06-02 08:14 +0200
            Re: Fast trunc() algorithm jak <nospam@please.ty> - 2022-06-02 11:14 +0200
              Re: Fast trunc() algorithm Bonita Montero <Bonita.Montero@gmail.com> - 2022-06-02 12:48 +0200

csiph-web