4
\$\begingroup\$
// ms is expected to be between 0 and 2000
// fakeDouble is expected to be {0, 0}
// uint64_t isn't used because it doesn't work with /NODEFAULTLIB when compiling 32-bit
void msToDouble(uint32_t ms, uint32_t* fakeDouble) {
    if (ms == 0) {
        return;
    } else if (ms >= 2000) {
        fakeDouble[1] = 0x40000000;
        return;
    }
    
    int exponent = 0;
    
    while (ms < 1000) {
        ms *= 2;
        exponent += 1;
    }
    
    fakeDouble[1] = (1023 - exponent) << 20;
    
    for (int bitPlace = 51; ms != 1000 && bitPlace >= 0; bitPlace--) {
        ms -= 1000 * (ms > 1000);
        ms *= 2;
        fakeDouble[bitPlace >= 32] |= (ms >= 1000) << (bitPlace - (32 * (bitPlace >= 32)));
    }
}

I read about compiling with the /NODEFAULTLIB option, and I'm planning on doing that with a small program I'm working on since it doesn't use much standard library stuff. This is a function I'll need to use due to not having access to atof or similar functions.

\$\endgroup\$
7
  • 6
    \$\begingroup\$ Why uint32_t* fakeDouble? Why not an actual double? \$\endgroup\$ Dec 10, 2023 at 20:12
  • \$\begingroup\$ I just need to make the bit representation of a double. \$\endgroup\$ Dec 10, 2023 at 22:39
  • 15
    \$\begingroup\$ Is there any reason you're not doing the obvious x % 1000 / 1000.? \$\endgroup\$ Dec 11, 2023 at 8:31
  • 9
    \$\begingroup\$ @justhalf Maybe, but OP should be explicit. My first guess was that they're thinking it'll be faster if they do it manually. \$\endgroup\$ Dec 11, 2023 at 15:49
  • 5
    \$\begingroup\$ "This is a function I'll need to use due to not having access to atof or similar functions" -- I don't see what atof has to do with integer-to-fp conversion. Arithmetic and type conversions are built into the language. They do not depend on the standard library. I suggest reserving custom code for stuff that actually requires it. \$\endgroup\$ Dec 12, 2023 at 15:29

3 Answers 3

12
\$\begingroup\$

random garbage in initial word

// fakeDouble is expected to be {0, 0}

This is insanity. The code needs to implement this comment. Either assert that caller supplied 64 zero bits, or assign zeros upon entry.

        fakeDouble[1] = 0x40000000;

Relying on the [0] element to be zero is far too adventurous. As written it is indeterminate what double value this function returns. We strive to avoid Undefined Behavior.

\$\endgroup\$
3
  • \$\begingroup\$ Does this look okay? pastebin.com/1cnpMVp6 \$\endgroup\$ Dec 10, 2023 at 22:39
  • 2
    \$\begingroup\$ Assigning 64 zero bits up front is terrific, yes that looks okay. It is doing fakeDouble[0] = 0; fakeDouble[1] = 0 \$\endgroup\$
    – J_H
    Dec 10, 2023 at 23:52
  • 5
    \$\begingroup\$ Note that undefined behavior has a very specific meaning when it comes to the C programming language, usually shortened as UB. It's confusing if you use the same expression to just mean "indeterminate". \$\endgroup\$
    – pipe
    Dec 12, 2023 at 16:46
11
\$\begingroup\$

Missing header

uint32_t is not defined - we need to include <stdint.h> for this to compile.

Unit tests

The tests are completely absent.

Explanation

The rationale for returning a pair of integers instead of simply returning ms / 1000.0 as a double is missing.

Magic numbers

The code seems to assume a 64-bit IEEE double format. Either we should be using the target platform's floating-point type, in which case we should include <float.h> and use the standard macros such as DBL_RADIX and DBL_MANT_DIG, or we should define the output format ourselves and define our own constants for the layout.

There's also an assumption that 1/1000.0 is not denormal. That's more forgivable, bit it's worth a comment that this code won't work for arbitrary divisors, in case someone attempts to re-use the function in a different context.

\$\endgroup\$
11
\$\begingroup\$

Code fails every time

With ms == 0, code can return a non-zero fake double due to unassigned fakeDouble[].

Non-zero ms errors approach 1 part in 10,000.

Robust code errors would be less than 1 part in 1015.


Sample test code

#include <stdio.h>

// Return error code
int msToDouble_int_test(uint32_t ms, int print) {
  union {
      double d;
      uint32_t fakeDouble[2];
  } u;
  msToDouble(ms, u.fakeDouble);
  double ref = ms/1000.0;
  if (ref != u.d) {
    static double worst = 0.0;
    double err = (u.d - ref)/ref;
    if (ms != 0  && err > worst) {
      worst = err;
    }
    if (print || err == worst) {
      printf("%5lu --> %-24.17g %-21a, expect %-24.17g %-21a   %.3g%%\n",
        (unsigned long)ms,u.d, u.d, ref, ref, err*100.0);
    }
    return 1;
  }
  return 0;
}

Sample output

    0 --> 4.9406564584124654e-323  0x1.4p-1071          , expect 0                        0x0p+0                  inf%
    1 --> 0.0009999999999999998    0x1.0624dd2f1a9fbp-10, expect 0.001                    0x1.0624dd2f1a9fcp-10   -2.17e-14%
    2 --> 0.0019999999999999996    0x1.0624dd2f1a9fbp-9 , expect 0.002                    0x1.0624dd2f1a9fcp-9    -2.17e-14%
    3 --> 0.0030000004843486728    0x1.89374fefbfffbp-9 , expect 0.0030000000000000001    0x1.89374bc6a7efap-9    1.61e-05%
    4 --> 0.0040000006408717993    0x1.0624dfefbfffbp-8 , expect 0.0040000000000000001    0x1.0624dd2f1a9fcp-8    1.6e-05%
    5 --> 0.0050000026822090106    0x1.47ae1fffffffbp-8 , expect 0.0050000000000000001    0x1.47ae147ae147bp-8    5.36e-05%
    6 --> 0.0060000009834766345    0x1.89374fffffffbp-8 , expect 0.0060000000000000001    0x1.89374bc6a7efap-8    1.64e-05%
    7 --> 0.0070000030100345568    0x1.cac08fffffffbp-8 , expect 0.0070000000000000001    0x1.cac083126e979p-8    4.3e-05%
    8 --> 0.0080000013113021764    0x1.0624dfffffffbp-7 , expect 0.0080000000000000002    0x1.0624dd2f1a9fcp-7    1.64e-05%
    9 --> 0.0090000033378600987    0x1.26e97fffffffbp-7 , expect 0.0089999999999999993    0x1.26e978d4fdf3bp-7    3.71e-05%
   10 --> 0.010000005364418021     0x1.47ae1fffffffbp-7 , expect 0.01                     0x1.47ae147ae147bp-7    5.36e-05%
   11 --> 0.011000007390975943     0x1.6872bfffffffbp-7 , expect 0.010999999999999999     0x1.6872b020c49bap-7    6.72e-05%
   22 --> 0.022000014781951887     0x1.6872bfffffffbp-6 , expect 0.021999999999999999     0x1.6872b020c49bap-6    6.72e-05%
   35 --> 0.035000026226043666     0x1.1eb85fffffffbp-5 , expect 0.035000000000000003     0x1.1eb851eb851ecp-5    7.49e-05%
   65 --> 0.065000057220458915     0x1.0a3d7fffffffbp-4 , expect 0.065000000000000002     0x1.0a3d70a3d70a4p-4    8.8e-05%
  125 --> 0.12500011920928941      0x1.00000fffffffbp-3 , expect 0.125                    0x1p-3                  9.54e-05%
  250 --> 0.25000023841857905      0x1.00000ffffffffp-2 , expect 0.25                     0x1p-2                  9.54e-05%
  500 --> 0.50000047683715809      0x1.00000ffffffffp-1 , expect 0.5                      0x1p-1                  9.54e-05%
 1000 --> 1.0000009536743162       0x1.00000ffffffffp+0 , expect 1                        0x1p+0                  9.54e-05%
 2000 --> 2.0000019073486324       0x1.00000ffffffffp+1 , expect 2                        0x1p+1                  9.54e-05%
Error count 2001

Alternate

Do the division using wide integer math to form the quotient. Below keeps track powers of 2 separately. Then form the floating point fake double.

If code adjusts the dividend / divisor prior to division and pays attention to rounding, there is no error for all uint32_t values [0 ... 0xFFFFFFFF].

#define DBL_IMPLIED_BIT (1ull << 52)

// Divide an integer by 1000 and form a floating point quotient
void msToDouble(uint32_t ms, uint32_t *fakeDouble) {
  unsigned long long dividend = ms;
  unsigned long long divisor = 1000;
  unsigned long long quotient = 0;
  unsigned long long expo = 0;
  if (dividend) {
    expo = (1024-1) + 52;  // Exponent for 1.0, scaled
    // Adjust dividend to maximize quotient precision.
    while (dividend < 0x8000000000000000u) {
      dividend <<= 1;
      expo--;
    }
    // Adjust divisor to maximize quotient precision.
    while (divisor % 2 == 0) {
      divisor >>= 1;
      expo--;
    }
    // Do the division
    quotient = dividend/divisor;
    while (quotient <  DBL_IMPLIED_BIT) {
      quotient <<= 1;
      expo--;
    }
    // Round based on last bit shifted out
    unsigned last_bit_out = 0;
    while (quotient >= DBL_IMPLIED_BIT*2) {
      last_bit_out = quotient & 1;
      quotient = quotient >> 1;
      expo++;
    }
    quotient += last_bit_out;
    // Detect carry
    if (quotient >=  DBL_IMPLIED_BIT*2) {
      quotient = quotient >> 1;
      expo++;
    }
    // Lop off implied bit
    quotient -= DBL_IMPLIED_BIT;
  }
  quotient |= expo << 52;
  fakeDouble[1] = quotient >> 32;
  fakeDouble[0] = quotient;
}
\$\endgroup\$
2
  • \$\begingroup\$ "Robust code errors would be less than 1 part in 10^15" I'm curious, where did you get that number from? \$\endgroup\$
    – Mast
    Dec 12, 2023 at 17:39
  • 1
    \$\begingroup\$ @Mast Typical double encodes most values with a 53-bit precision: 1 part in 2^53. This is akin to at least 15 decimal digits. I would expect ms/1000.0 per a robust msToDouble() to be correct to at least that many significant digits. \$\endgroup\$ Dec 12, 2023 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.