The Thomae function defined on `double`

The Thomae function $f$ is defined as $$f(x) = \left\{ \begin{array}{cl} 0 & \textrm{if $x$ is irrational.} \\ \frac{1}{n} & \textrm{if $x = \frac{m}{n}$ with $\gcd(m, n) = 1$ and $n > 0$.} \\ \end{array}\right.$$ Because all floating-point numbers are base two (and therefore rational), it therefore follows that the Thomae function evaluated at any double-precision floating-point number must be representable by a non-zero double-precision floating-point number. This makes this function definitely less-than-interesting, mathematically, because $\frac{1}{3}$ cannot be represented as a double precisely, and as the representing in binary is a repeating decimal, thomae( 1.0/3.0 ) equals 5.55112e-17, or $2^{-54}$ as the representation is 3fd5555555555555, and thus an unormalized representation would be $0.0000000000001_{16} \times 2^{-2}$, and normalizing this produces $1.0000000000000_{16} \times 2^{-2 - 52}$, as there are $52$ bits in the mantissa.

Here is an implementation, but before you look at it, consider that this is an interesting exercise in bit manipulation, for if the number is $1.0010110010111101 \times 2^n$, the result is $0.0000000000000001 \times 2^n$, and thus, the result should be a double with all zeros in the mantissa or a subnormal number with only one bit set to 1.

The three cases you must deal with are where:

The argument is a normalized double and the result is a normalized double.
The argument is a normalized double but the result is a subnormal double.
The argument is a subnormal double in which case the result must be subnormal double.

Please note, the "fractional" part of a real number are all bits to the right of the radix point, so given $101.110101$ and $0.00101110101$, the fractional parts are $0.110101$ and $0.00101110101$, respectively. The "mantissa", however, is all the bits to the right of the radix point when the number is written in scientific notation: $101.110101 = 1.01110101 \times 2^2$ and $0.00101110101 = 1.01110101 \times 2^{-3}$, and so for both of these, the mantissa is $0.01110101$. In the literature, all the significant bits, including the leading "1" are described as the "significand", so for both these numbers, the significand is $1.01110101$. The double double-precision floating-point representation stores $52$ bits in the mantissa, and therefore the significand has $53$ bits, as there is also the implied leading $1$ when the number is normalized. For subnormal numbers, the number of significant bits is less than $53$, possibly as low as $1$ bit when storing $2^{-1022 - 52}$.

Good luck, enjoy, and here is the implementation:

#include <cstdint>
#include <cmath>
#include <cassert>

double thomae( double x );

double thomae( double x ) {
    // The Thomae function is not defined for infinity or not-a-number
    if ( std::isnan( x ) || std::isinf( x ) ) {
        return NAN;
    }

    // thomae( -x ) == thomae( x )
    if ( x < 0.0 ) {
        x = -x;
    }

    // A 'double' is 8 bytes or 64 bits, with
    //    1 bit for the sign,
    //   11 bits for the exponent, and
    //   52 bits in the mantissa.
    constexpr int mantissa_bits{ 52 };

    // Break the double 'x' into its integer and fractional parts
    //  - In this case, discarding the integer part
    double integer_part;
    x = std::modf( x, &integer_part );

    // If the fractional part is zero,
    //    the double stores an integer value 'n', and thomae( n ) == 1
    if ( x == 0.0 ) {
        return 1.0;
    }

    // Ensuring that 0 < x < 1
    assert( (0.0 < x) && (x < 1.0) );

    using u64 = std::uint64_t;

    // Convert the double to a long in order to perform bit operations
    //  - Bit operations are not defined on floating-point types
    u64 n{ std::bit_cast<u64>(x) };

    // The mantissa are the last 52 bits
    //  - These are all bits to the right of the decimal point
    //    when written using scientific notation
    //                     exponent
    //      1.mmmmm...m x 2
    u64 mantissa{ n & ((u64(1) << mantissa_bits) - 1) };

    // If the mantissa is zero, then thomae( x ) == x,
    //                                                -exponent 
    //   as there is only the leading bit, so x == 1/2         ,
    //   recalling that the exponent at this point must be negative.
    if ( mantissa == 0 ) {
        return x;
    }

    // Extract the exponent from the representation
    //  - We took the absolute value, so the sign bit is '0'.
    u64 exponent{ n >> mantissa_bits };

    // Count the number of trailing zeros in the mantissa.
    //  - There must be at least one '1', so this must be less than 52.
    int trailing_0s{ std::countr_zero( mantissa ) };

    if ( exponent == 0 ) {
        // The argument was a subnormal double,
	//       so just return a '1' in the location of the last trailing '1'
	//  - We previously set the sign bit to '0'
	//     S 00000000000 0000100000001000000011000000100000000001000010000000
	//  -> 0 00000000000 0000000000000000000000000000000000000000000010000000
	//                                       trailing one -----------^
        return std::bit_cast<double>( u64(1) << trailing_0s );
    } else {
        if ( exponent <= (mantissa_bits - trailing_0s) ) {
            // The result is a subnormal double
	    //  - This is subtle, but correct... Why?
            return std::bit_cast<double>( u64(1) << (
                (mantissa_bits - 1) - ((mantissa_bits - trailing_0s) - exponent)
            ) );
        } else {
            // The result is a normalized double
            exponent -= mantissa_bits - trailing_0s;
            return std::bit_cast<double>( exponent << mantissa_bits );
        }
    }
}

The Thomae function defined on double

The Thomae function defined on `double`