Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.c++ > #91407

Re: How can I get an implementation of fmod()?

Show all headers | View raw

Am 25.08.2023 um 08:25 schrieb Megumin:
> How does remainder() in cmath works and is there any simple but fairly precise fmod that can be implemented without using double or long double?
> I'm currently trying to implement fmodf(float x,float y) myself in assembly, and seems like in IEEE754, if I try to do as follows, the precision will f-up badly if x>>y in fact if x/y is more than 2^23 this will just makes no sense at all.
> r1 = f32(x);
> r2 = f32(y);
> r3 = div_f32(r1,r2);        //precision:2ulp
> r4 = trunc(r3)// this is where things go bad if r3 is more than 2^23 it does nothing and I will get no remainder.
> r4 = multi_f32(r4,r2);
> r4 = subtract_f32(r1,r4);
> result = r4;
> 

Here's some fmod()-code I implemented myself for float and double.
It always produces results without precision loss because the result
is always smaller than the dividend.

template<floating_point Float>
static Float myFmodT( Float counter, Float denom )
{
	constexpr bool IS_DOUBLE = same_as<Float, double>;
	static_assert(IS_DOUBLE || same_as<Float, float>);
	constexpr unsigned
		BITS = IS_DOUBLE ? 64 : 32,
		MANT_BITS = IS_DOUBLE ? 52 : 23,
		EXP_BITS = IS_DOUBLE ? 11 : 8;
	static_assert(sizeof(Float) * 8 == BITS && 
numeric_limits<Float>::is_iec559);
	using bin_t = conditional_t<IS_DOUBLE, uint64_t, uint32_t>;
	constexpr int MAX_EXP = (1 << EXP_BITS) - 1;
	constexpr bin_t
		SIGN_BIT = (bin_t)1 << BITS - 1,
		EXP_MASK = (bin_t)MAX_EXP << MANT_BITS,
		IMPLCIT_BIT = (bin_t)1 << MANT_BITS,
		MANT_MASK = IMPLCIT_BIT - 1;
	static auto bin = []( Float f ) -> bin_t { return bit_cast<bin_t>( f ); };
	static auto flt = []( bin_t b ) -> Float { return bit_cast<Float>( b ); };
	bin_t
		bCounter = bin( counter ),
		bDenom = bin( denom ) & ~SIGN_BIT,
		bSign = bCounter & SIGN_BIT;
	bCounter &= ~SIGN_BIT;
	static auto maxExp = []( bin_t b ) -> bool { return b >= EXP_MASK; };
	static auto infMant = []( bin_t b ) -> bool { return !(b & MANT_MASK); };
	auto floatingResult = [&]() -> Float { return (counter * denom) / 
(counter * denom); };
	if( !bDenom ) [[unlikely]]
		// ... % +/-0.0
		return floatingResult();
	if( maxExp( bCounter ) ) [[unlikely]]
		// [Inf|NaN|SNaN] % non-zero
		return floatingResult();
	if( maxExp( bDenom ) ) [[unlikely]]
		if( infMant( bDenom ) ) [[likely]]
			// finite % Inf
			return counter;
		else
			// finite % [Inf|NaN|SNaN]
			return floatingResult();
	if( !bCounter ) [[unlikely]]
		// +/-0.0 % non-zero-finite = +/-0.0
		return counter;
	int
		counterExp = bCounter >> MANT_BITS & MAX_EXP,
		denomExp = bDenom >> MANT_BITS & MAX_EXP;
	bin_t
		counterMant = (bin_t)(bool)counterExp << MANT_BITS | bCounter & MANT_MASK,
		denomMant = (bin_t)(bool)denomExp << MANT_BITS | bDenom & MANT_MASK;
	static auto normalize = []( bin_t &mant, int &exp )
	{
		unsigned bits = countl_zero( mant ) - EXP_BITS;
		mant <<= bits;
		exp -= bits;
	};
	if( !counterExp ) [[unlikely]]
		// normalize counter
		normalize( counterMant, counterExp ),
		++counterExp;
	if( !denomExp ) [[unlikely]]
		// normalize denominator
		normalize( denomMant, denomExp ),
		++denomExp;
	int remExp = counterExp;
	bin_t remMant = counterMant;
	for( ; ; )
	{
		int below = remMant < denomMant;
		if( remExp - below < denomExp ) [[unlikely]]
			break;
		remExp -= below;
		remMant <<= below;
		if( !(remMant -= denomMant) ) [[unlikely]]
			return flt( bSign );
		normalize( remMant, remExp );
	}
	if( remExp <= 0 ) [[unlikely]]
		// denormal result
		remMant >>= -remExp + 1,
		remExp = 0;
	return flt( bSign | (bin_t)remExp << MANT_BITS | remMant & MANT_MASK );
}

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

Thread

How can I get an implementation of fmod()? Megumin <chrisramirez0227@gmail.com> - 2023-08-24 23:25 -0700
  Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 09:19 +0200
    Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 09:35 +0200
  Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 09:49 +0200
  Re: How can I get an implementation of fmod()? James Kuyper <jameskuyper@alumni.caltech.edu> - 2023-08-25 04:34 -0400
    Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 10:37 +0200
      Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 10:43 +0200
        Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-25 11:20 +0200
          Re: How can I get an implementation of fmod()? Megumin <chrisramirez0227@gmail.com> - 2023-08-30 22:52 -0700
            Re: How can I get an implementation of fmod()? Bonita Montero <Bonita.Montero@gmail.com> - 2023-08-31 11:00 +0200
            Re: How can I get an implementation of fmod()? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-08-31 21:48 +0100
            Re: How can I get an implementation of fmod()? Keith Thompson <Keith.S.Thompson+u@gmail.com> - 2023-08-31 14:37 -0700
    Re: How can I get an implementation of fmod()? James Kuyper <jameskuyper@alumni.caltech.edu> - 2023-08-25 04:44 -0400

csiph-web